GSEA란 Gene Set Enrichment Analysis의 줄임말로 특정 유전자들의 집합이 있으면 이 유전자들이 어떠한 특성을 가지는 지 알아보는 분석 방법입니다. RNA 시퀀싱의 downstream analysis로 많이 쓰입니다.
간단하게 GSEA를 하는 방법은, 내가 원하는 Gene들을 p-value로든 Fold change로든 어떠한 점수를 기준으로 1위부터 꼴등까지 줄을 세웁니다. 그 후에 그 순위를 가지고 어떠어떠한 Gene Set에는 높은 순위의 유전자가 많이 들어있더라, 낮은 순위의 유전자가 많이 들어있더라 등을 확인합니다. 여기서 Gene set은 특정 pathway에 속하는 유전자들일 수도 있고, 공통된 기능을 가진 유전자 그룹일 수도 있고, 특정 질병에 노출되면 발현되는 유전자 리스트일 수도 있습니다. 사용자가 정하기 나름이죠.
Broad Institute와 UC San Diego에서 이러한 유전자 세트, Gene set을 데이터베이스화 해서 정리해놓은 사이트가 있는데 그게 MsigDB입니다. 인간과 쥐의 유전자를 그룹별로 잘 정리해놓아서, GSEA에 유용하게 쓸 수 있습니다.
https://www.gsea-msigdb.org/gsea/msigdb/

GSEA는 보통 그룹 간의 DEG (Differentially Expressed Gene), 즉 유의미하게 다르게 발현된 유전자들을 가지고 실행할 수가 있습니다. 이 때에 실행하는 GSEA는 내가 뽑은 DEG들이 어떠한 성격을 지니고 있는지 알아보는 것이겠죠. 예를 들어서, 두 환자 그룹에 한 그룹에는 비타민D를 1개월간 복용시키고 다른 그룹에는 플라시보를 복용시켜서 유전자 발현의 차이를 보는 경우에 두 그룹 간의 DEG를 구해 GSEA를 돌려 본다면 비타민D로 인하여 유의미하게 변한 유전자가 대체로 어떠한 그룹에 속하며 어떤 기능을 가지고 있는 지 알아볼 수 있습니다. 이 때, 이 DEG가 10개-20개 정도면 뭐 하나하나 살펴볼 수 있겠지만, 100개에서 500개 이렇게 나온다면 하나하나 보는 건 너무 많은 시간과 노력이 들겠죠. 그러면 GSEA가 이러한 시간을 아주 줄여줄 수 있습니다.
비슷하지만 다른 방법의 GSEA로는, DEG만 넣는 게 아니라 모든 유전자의 발현 정보를 넣어서 GSEA를 실행해볼 수도 있습니다. 이 때에는 입력 값이 아주 크겠지만, 위의 예제로 본다면 비타민D가 전반적으로 어떠한 영향을 미치고 어떤 pathway가 변화하는 지 종합적으로 알 수가 있습니다.
위의 두 방법 중에 어느 것이 옳고 어느 것이 틀리고 그러지는 않습니다. 다만, 내가 알아보려고 하는 질문에 따라서 전자 혹은 후자를 택하기도 합니다.
이 글에서는 Single cell RNA-Seq 분석을 할 때에 그 Downstream analysis로 Seurat 객체를 가지고 GSEA를 하는 방법을 간단히 기록하려고 합니다.

일단 먼저 Seurat을 불러오고 내가 이용하려는 데이터를 불러옵니다. 또한, 우리가 하려는 GSEA는 모든 유전자를 줄세워서 입력하는 위에 적은 후자의 방법을 쓰기 때문에, 유전자를 그룹 간에 비교한 값이 필요합니다. Seurat으로는 2만개 이상의 유전자를 모두 비교하는 데 엄청나게 오래 걸리기 때문에, presto 패키지에 있는 Wilcoxon rank sum test를 불러와서 사용하려고 합니다.

library(Seurat)
library(devtools)
install_github("immunogenomics/presto")   #인스톨을 최초로 한 다음에는 이 줄은 생략합니다
library(presto)

#또 나중에 필요한 라이브러리 들을 여럿 불러와 보겠습니다.
library(msigdbr)
library(fgsea)
library(tidyr)
library(dplyr)
library(ggplot2)
library(tibble)
library(tidyverse)
library(data.table)

그런 다음, Wilcoxon rank sum test를 모든 유전자 대상으로 수행해줍니다. 제 Seurat 데이터는 myData3에 저장되어 있으며 저는 그 중 celltype.condition에서 특정 세포 (CD14+ Monocytes)의 질병 그룹(cALD)과 건강한 비교군 (Control) 간의 차이를 보고 싶습니다.

DefaultAssay(myData3)<-"RNA"
Idents(myData3)<-"celltype.condition"
myData.genes<-wilcoxauc(myData3, "celltype.condition", 
				groups_use = c("CD14+ Mono_cALD","CD14+ Mono_Control"))

이제 앞서 말한 mSigDB에서 Gene Set들을 불러와야 합니다. 저는 가장 기본으로 유전자를 50개의 질병 및 생명 현상으로 나눈 Hallmark gene sets를 먼저 보려고 합니다. 그리고 이 불러온 Gene Set을 GSEA에 필요한 정보만 뽑아서 정리해 줍니다.

#msigdbr_show_species() 
#위의 코드는 모든 가능한 생물 종을 표시해 줍니다.

#"H"로  Hallmark gene set을 불러 옵니다.
m_df<- msigdbr(species = "Homo sapiens", category = "H") 
#불러온 Gene Set을 정리해 줍니다.
fgsea_sets<- m_df %>% split(x = .$gene_symbol, f = .$gs_name)

#불러온 Gene Set을 살펴봅니다.
head(m_df)
dplyr::count(myData.genes, group)

그런 다음에 제가 원하는 그룹을 뽑아서 logFC와 AUC score를 빼낸 뒤 정리해봅니다. 저는 CD14+ Mono_cALD 그룹에서 어떠한 pathway들이 더 혹은 덜 발현되는지 알아보고 싶습니다.

myData.genes %>%
  dplyr::filter(group == "CD14+ Mono_Control") %>%
  arrange(desc(logFC), desc(auc)) %>%
  head(n = 10)

#그리고 뽑은 정보들을 auc에 따라 1등부터 꼴등까지 줄을 세워봅니다.
CD14Mono.genes<- myData.genes %>%
  dplyr::filter(group == "CD14+ Mono_Control") %>%
  arrange(desc(auc)) %>% 
  dplyr::select(feature, auc)

ranks<- deframe(CD14Mono.genes)
head(ranks)

이제 실제로 GSEA를 실행합니다. 저희는 이번에 fgsea라는 패키지를 사용할 겁니다. 바이오컨덕터로 간편하게 설치 가능합니다.
https://bioconductor.org/packages/release/bioc/html/fgsea.html

fgsea는 인풋으로
1) 미리 지정한 Gene Set 이 필요하구요 (m_df로 불러들여서 fgsea_sets로 포맷을 맞춤)
2) 유전자의 리스트와 각각 유전자의 순위 혹은 순위를 지정할 수 있는 점수가 필요합니다 (myData3에서 auc를 계산하여 myData.genes에 저장, 그런 다음 CD14Mono.genes에 순서대로 저장)

fgseaRes<- fgsea(fgsea_sets, stats = ranks, nperm = 1000)
fgseaResTidy <- fgseaRes %>%
  as_tibble() %>%
  arrange(desc(NES))

#결과를 출력하기 좋게 정리합니다
fgseaResTidy %>% 
  dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% 
  arrange(padj) %>% 
  head()

이제 마지막 순서입니다. ggplot을 이용하여 결과를 나타내어 봅니다.

ggplot(fgseaResTidy %>% filter(padj < 0.05) %>% head(n= 50), aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill= NES<0)) +
  coord_flip() +
  labs(x="Pathway", y="Normalized Enrichment Score",
       title="Hallmark pathways NES from GSEA") + 
  theme_minimal()

다음과 같은 결과가 나옵니다.

각각 Pathway나 Gene Set에 Enrichment Plot을 보려면 다음 코드가 유용합니다.

plotEnrichment(fgsea_sets[["HALLMARK_TNFA_SIGNALING_VIA_NFKB"]],
               ranks) + labs(title="HALLMARK_TNFA_SIGNALING_VIA_NFKB")

 
참 쉽죠??
끝!!!

이 글은 사이언스지의 기사인 "New animal study raises concerns about high-does gene therapy"를 번역한 글입니다. 원문은 여기서 찾아보실 수 있습니다.

치명적인 신경 근육 질환을 가지고 태어난 아기에게 극적인 효과를 가져온 최근의 유전자 치료 시험이 비슷한 방법을 이용한 다른 질병의 치료에도 사용될 수 있는 희망을 불러 왔습니다. 그러나 새로운 동물 실험에 따르면, 이러한 치료에 사용되는 유전자를 포함한 바이러스의 고용량 투여가 인간 임상 실험에서 제시된 것처럼 항상 안전하지만은 않을 수도 있다고 합니다.

어제 몇몇 유전자 치료 회사의 주식에 급락을 불러일으킨 이 새로운 연구에서, 연구자들은 보통 해가 없어서 신경 근육 질환을 치료하기 위한 치료 유전자를 세포로 전달해주는 데에 쓰이는 아데노 관련 바이러스 9(adeno-associated virus 9; AAV9)를 어린 원숭이와 돼지에게 주입했습니다. 며칠 내로, 일부의 동물들은 심각한 간과 신경 손상을 일으켰습니다.

이 연구 결과는 10대 소년이 유전자 치료 벡터로 인한 면역 반응으로 사망한 1999년 임상 시험을 주도한 펜실베니아 대학의 유전자 치료 원구원 제임스 윌슨 (James Wilson) 교수의 연구실에서 주도했기 때문에 더욱 주목을 받았습니다. 이러한 최근의 동물 실험 결과를 바탕으로 Wilson은 고용량 AAV9 투여를 이용한 유전자 치료를 계획하고 있는 연구자들에게 비슷한 부작용을 찾아보아야 한다고 경고하고 있습니다.

그러한 경고에도 불구하고, 유전자 치료 분야의 그와 다른 사람들은 이것이 유전자 운반자 (gene carrier)에 대한 진행 중인 임상 실험이 중단되어야 한다는 것을 의미하지 않는다고 말합니다. 지난 1월 29일 새로운 연구를 온라인으로 발표한 ​Human Gene Therapy​의 편집장이자 워체스터 메사추세츠 의과 대학의 유전자 치료 연구자인 테렌스 플롯 (Terence Flotte)는 "이번 연구는 잠재적으로 큰 의미가 있다"고 말하면서, 자신의 논평에서 "지나치게 반응해서는 안된다"라고 썼습니다.

DNA를 전달하기 위해 유전자 치료에 쓰이는 벡터 중에서, AAV9는 신경 조직을 통해 퍼지는 데 특히 탁월한 것으로 입증 되었고, 고용량 투여시 혈액-뇌 장벽을 통과하여 척수 신경 및 뇌에 도달 할 수도 있습니다. 지난 11월에 발표된 AAV9를 이용한 최초의 인체 실험에서, 일반적으로 2세 이전에 사망하는 신경퇴행성 질환인 심각한 형태의 척수성 근 위축증을 앓고 있는 15명의 아기가 그들이 결핍된 유전자인 ​SMN​을 삽입한 바이러스를 주입받았습니다. 대부분의 아기들은 이제 앉을 수 있고, 그 중 두 명은 현재 걸어다닙니다.

그러나 Wilson 교수의 연구에서 SMN​이 삽입된 비슷한 AAV9을 정맥 내로 대용량 주사한 3 마리의 어린 히말라야 원숭이들은 간 독성 징후를 보여, 그 중 하나는 안락사를 받아야 했습니다. 동일한 치료를 받은 3 마리의 새끼 돼지들은 운동 신경에 손상이 있었고, 더 이상 걸을 수가 없어, 또한 안락사를 받아야 했습니다. 이러한 반응들은 이전의 AAV 연구에서 보여진 부작용과는 달리, AAV9의 단백질 외막이나 AAV9가 보유하고 있는 유전자에 대한 면역 반응과는 관련이 없는 것으로 보인다고 Wilson 교수의 그룹이 보고했습니다.

이 연구실은 다른 유전자를 지니는 또다른 AAV9 변종을 한 원숭이에 주입한 결과 비슷한 부작용을 보인다는 것을 발표할 계획으로, 그 결과가 광범위하게 적용 될 것이라고 주장합니다. Wilson 교수와 그 동료 연구팀은 적어도 5개의 계획되어 있는 혹은 진행 중인 임상 시험에서 유전자를 근육이나 신경세포에 전달하기 위해 AAV9 혹은 그 변종의 대용량 투여를 이용한다고 언급했습니다. 동물에서의 이러한 연구를 포함한 다른 안전성 연구들도 "유사한 독성에 대한 주의 깊은 모니터링을 포함시켜야 한다"고 그들은 주장합니다. 그의 발견이 이번달 초에 발표될 즈음에, Wilson 교수는 Duchenne 근이영양증 (Duchenne muscular dystrophy)를 치료하기 위해 정맥주사 AAV9 유전자 치료법을 사용할 계획인 Solid Biosciences 이사회에서 최근 사임하였습니다. (이 치료법에 관한 임상 실험은 미국 식품의약국[FDA]에 의해 정밀 조사 중이며, 현재는 고용량의 AAV9를 사용하는 것을 사실 상 중단하였습니다.)

Flotte는 이 새로운 결과들이 인체 실험에 적용되지 않을 수도 있는 몇 가지 이유에 대해 다음과 같이 설명합니다. Wilson 교수 팀이 주사한 AAV9은 임상 시험에서 사용된 것과 동일한 타입이 아니며, 그 연구는 동물에게 독성을 유발할 지도 모르는 인간 버전의 ​SMN​ 유전자를 사용했습니다. 또한, 연구자들은 치료가 오염되었다는 가능성을 배제하지 않았다고 말했습니다. Flotte와 공동 저자들은 그들의 논평에서 AAV9 유전자 치료 임상 시험이 잠재적으로 생명을 구하기 때문에 계속되어야 한다고 기술하고 있습니다. "이걸 연구해 봅시다, 하지만 그걸 무시하지는 맙시다."라고 Flotte는 덧붙였습니다. FDA 대변인은 이 연구에 대해 진술을 거부했습니다.

Wilson 교수는 또한 "이 시점에서 무엇이든 주제에서 벗어나고 싶지 않다"고 말했습니다. 그는 새로운 실험이 AAV9 안전성에 대한 이전의 연구와 크게 다른 점에 동의했습니다. 또한, 실험실과 벡터 생산 시설은 AAV9 용량을 측정하기 위해 서로 다른 방법을 사용하기 때문에 여러 연구에 걸쳐서 용량을 비교하기가 어렵습니다. Wilson 교수는 "단지 이것이 (새로운 유전자 치료를 개발하는 연구자들에게 제공되는) 유용한 정보가 되기를 희망한다"고 말했습니다.

채플 힐에있는 노스 캐롤라이나 대학 (University of North Carolina)의 유전자 치료 연구원인 Jude Samulski는 AAV9를 이용한 영장류 연구에서 이러한 독성 영향을 한번도 본 적이 없으며 앞으로 많은 실험실에서 Wilson 교수의 벡터를 요청하고 그의 실험을 반복 할 것이라고 예측했습니다. "이것이 진정으로 우려된다면 재현이 가능해야합니다."라고 그는 말합니다.