GSEA·fgsea 완벽 이해 — enrichment score부터 MSigDB 실습까지 (한국어)

volcano plot에 500개 점이 찍혔습니다. 그래서 무슨 일이 일어난 건가요

S19와 S20을 거치며 우리는 유의 유전자 리스트를 손에 넣었습니다. 그런데 리뷰어도, 공동연구자도, 심지어 여러분 자신도 묻습니다. "이 500개 유전자가 대체 무슨 이야기를 하는 거죠?" 유전자 이름 500개를 나열하는 논문은 없습니다. 논문은 이렇게 씁니다. "종양에서 세포주기 경로가 활성화되고 면역 반응 경로가 억제되었다."

개별 유전자에서 **경로(pathway)**로 시야를 올리는 것. 이게 "FASTQ to Paper" 여정의 마지막 조각이며, 이 편의 목표입니다. 두 가지 통계 철학 — ORA와 GSEA — 를 구분하는 것이 핵심입니다.

첫 번째 발상: 그냥 세어보자 (ORA)

가장 직관적인 방법은 과대표현 분석(Over-Representation Analysis)입니다. "내 유의 유전자 500개 중 세포주기 유전자가 몇 개인가? 우연보다 많은가?"를 초기하 분포로 검정합니다.

P(X \ge k) = \sum_{i=k}^{n} \frac{\binom{K}{i}\binom{M-K}{n-i}}{\binom{M}{n}}

전체 유전자 $M$ 개, 그중 세포주기 경로 유전자 $K$ 개, 내 유의 리스트 $n$ 개, 그 안의 세포주기 유전자 $k$ 개. Fisher exact test와 같은 계산입니다.

ORA는 간명하지만 치명적 약점이 있습니다. 임계선(padj < 0.05)을 넘긴 유전자만 세고, 나머지는 버립니다. 경로 전체가 조금씩(각각 1.3배씩) 일관되게 올라갔지만 아무도 임계선을 못 넘긴 경우, ORA는 그 경로를 완전히 놓칩니다. 생물학에서 가장 흥미로운 신호가 바로 이런 "약하지만 일관된 이동"입니다.

두 번째 발상: 순위 전체를 보자 (GSEA)

Gene Set Enrichment Analysis는 임계선을 없앱니다. 유전자 20,000개 전부를 fold change(또는 signed statistic) 기준으로 한 줄로 세운 뒤, 특정 경로의 유전자들이 그 줄의 위쪽에 몰려 있는지를 봅니다.

핵심 도구는 running enrichment score입니다. 정렬된 리스트를 위에서 아래로 훑으며,

경로에 속한 유전자를 만나면 점수를 올리고(가중치는 그 유전자의 통계량 크기에 비례),
경로 밖 유전자를 만나면 점수를 내립니다.

ES = \max_i \left( \sum_{g_j \in S, j \le i} \frac{|r_j|^p}{N_R} - \sum_{g_j \notin S, j \le i} \frac{1}{N - N_S} \right)

경로 유전자가 리스트 최상단에 몰려 있으면 누적합이 초반에 크게 치솟고, ES(그 최댓값)가 커집니다. 무작위로 흩어져 있으면 오르내림이 상쇄되어 ES가 0 근처에 머뭅니다. 이 곡선은 Kolmogorov-Smirnov 통계량과 정확히 같은 구조입니다 — 두 분포의 최대 편차를 보는 것.

손으로 running score 따라가기

유전자 8개를 순위대로 세웠고, 관심 경로 집합 $S = \{g_1, g_2, g_5\}$ (가중치 단순화를 위해 hit는 +0.33, miss는 −0.20이라 가정)라고 해봅시다.

순위	유전자	경로원소?	증감	누적합
1	g1	● hit	+0.33	0.33
2	g2	● hit	+0.33	0.66
3	g3	miss	−0.20	0.46
4	g4	miss	−0.20	0.26
5	g5	● hit	+0.33	0.59
6	g6	miss	−0.20	0.39
7	g7	miss	−0.20	0.19
8	g8	miss	−0.20	−0.01

누적합의 최댓값 ES = 0.66(순위 2에서). 경로 유전자가 상단에 몰렸으니 초반에 봉우리가 섰습니다. 이 봉우리 높이가 통계량이고, 봉우리 위치는 "어느 유전자까지가 핵심 기여자인가"(leading edge)를 알려줍니다.

유의성 — permutation

ES가 우연인지 판정하려면 귀무분포가 필요합니다. GSEA는 유전자 라벨(또는 표현형 라벨)을 무작위로 섞어 ES를 수천 번 다시 계산하고, 실제 ES가 그 분포의 어디에 위치하는지로 p-값을 뽑습니다. 경로마다 크기가 다르므로 ES를 정규화한 NES(Normalized ES)로 경로 간 비교를 합니다.

여기서 옛 GSEA의 고통이 있었습니다. permutation 수천 번 × 경로 수천 개 = 몇 시간. fgsea(fast GSEA)는 다단계 몬테카를로(adaptive multilevel)로 이 계산을 수백 배 가속해, 극히 작은 p-값도 정확히 추정합니다. 오늘날 실무 표준입니다.

fgsea로 MSigDB Hallmark 돌리기

이제 실습입니다. S19의 DESeq2 결과를 그대로 받아 경로 분석까지 이어봅시다.

1단계: 순위 벡터 만들기

if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("fgsea", "msigdbr"))
library(fgsea); library(msigdbr)

# res = S19의 DESeq2 결과 (log2FoldChange, stat 컬럼 포함)
# signed statistic으로 순위 벡터 생성 (NA 제거, 유전자명 = SYMBOL)
ranks <- res$stat
names(ranks) <- rownames(res)
ranks <- sort(ranks[!is.na(ranks)], decreasing = TRUE)

순위 기준으로 stat(Wald 통계량)을 쓰는 것이 fold change만 쓰는 것보다 안정적입니다. 크기와 신뢰도를 함께 반영하기 때문입니다.

2단계: 경로 세트 로드 + 실행

# MSigDB Hallmark 50개 경로 (사람)
h <- msigdbr(species = "Homo sapiens", category = "H")
pathways <- split(h$gene_symbol, h$gs_name)

set.seed(42)
fgseaRes <- fgsea(pathways = pathways, stats = ranks,
                  minSize = 15, maxSize = 500)
head(fgseaRes[order(padj)], 10)

3단계: 대표 경로 그림

topUp <- fgseaRes[NES > 0][order(padj)][1, pathway]
plotEnrichment(pathways[[topUp]], ranks) +
  ggplot2::labs(title = topUp)

plotEnrichment이 그려주는 곡선이 바로 위에서 손으로 따라간 running score입니다. 봉우리가 왼쪽(고발현 쪽)에 서면 그 경로는 종양에서 활성화된 것입니다. 유방암 데이터라면 HALLMARK_E2F_TARGETS, HALLMARK_G2M_CHECKPOINT 같은 세포주기 경로가 상단에 뜹니다. 우리가 S19에서 본 BRCA1·MKI67이 바로 그 경로의 구성원입니다. 이야기가 닫힙니다.

CS 매핑

Kolmogorov-Smirnov 통계량: ES는 경로 유전자의 누적분포와 균등분포의 최대 편차입니다. KS 검정의 바이오 응용 그 자체입니다.
running-sum / 최대 부분합: enrichment score 곡선은 누적합의 최댓값을 찾는 것으로, Kadane 알고리즘 계열의 최대 부분배열 문제와 형제입니다.
permutation test: 해석적 귀무분포 대신 라벨 셔플로 분포를 만드는 것은 비모수 통계의 핵심이며, fgsea의 multilevel은 rare-event 추정을 위한 몬테카를로 가속입니다.

자주 만나는 결함

유전자 ID 불일치 — 순위 벡터는 SYMBOL인데 경로 세트는 ENTREZ면 매칭이 0입니다. minSize 경고를 무시하지 맙시다.
fold change만으로 순위 — 저발현 유전자의 극단 LFC가 상단을 오염시킵니다. stat(shrunk LFC/SE) 사용을 권합니다.
ORA와 GSEA를 혼동 — "유의 유전자만 넣는" clusterProfiler enrichKEGG는 ORA입니다. 임계선 없는 전체 순위 분석은 GSEA/fgsea입니다. 목적에 맞게 고릅시다.
다중검정 미조정 — 경로 수천 개를 검정하므로 padj(BH)로 판단합니다. pval로 보고하면 위양성 경로가 쏟아집니다.

더 깊게 파고 싶다면

본문은 BPD가 직접 재구성한 서술입니다. 심화는 아래 정통 자료를 활용합시다.

Harvard STAT115 — Xiaole Shirley Liu 교수의 Pathway and Gene Set Enrichment Analysis (자막 완비, 자동번역 우수). ORA와 GSEA의 통계적 구분이 명료합니다.
원 논문(GSEA): Subramanian et al. (2005), Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles, PNAS 102:15545. running score의 원천입니다.
원 논문(fgsea): Korotkevich, Sukhov, Sergushichev (2019), Fast gene set enrichment analysis, bioRxiv 060012. multilevel 가속 유도입니다.
데이터베이스: MSigDB (gsea-msigdb.org) — Hallmark(H), C2 canonical pathways(KEGG/Reactome), C5(GO)를 무료 제공합니다. 라이선스는 연구용 CC 계열이니 상업 사용 시 조건 확인이 필요합니다.

이로써 원시 FASTQ가 논문 그림과 경로 이야기로 완성되는 "FASTQ to Paper" 여정이 닫힙니다. 다음 편 S22부터는 완전히 다른 무대 — 참조 게놈 없이 밑바닥부터 서열을 조립하는 de novo 어셈블리 — 로 넘어갑니다.

BioPlayground

GSEA와 fgsea: 유전자 하나가 아니라 경로 전체가 움직였다고 말하는 법