"FASTQ to Paper" 시리즈 하이라이트 — 우리는 지금 어디에 서 있나

• S17에서 우리는 원시 FASTQ를 Salmon으로 유사 정렬해 유전자별 카운트를 뽑았다.
• S18에서는 featureCounts와 RSEM 정량 방식을 비교했다.
• 이 편(S19)에서 우리는 그 카운트를 실제 논문 그림 — volcano plot — 으로 완성한다.
• 다음 S20에서는 edgeR / limma-voom 대안을, S21에서는 GSEA 경로 분석을 배운다.

즉 이 편의 실습을 마치면, 원시 FASTQ가 논문 게재용 그림으로 변신하는 전 여정 중 가장 중요한 통계 파트를 손에 넣게 된다. Wet-lab 연구원이 논문을 쓸 때 반드시 만나는 벽이다.

왜 정규분포가 아니라 음이항인가

먼저 우리 데이터를 이해하자. RNA-seq 카운트 매트릭스는 다음과 같이 생겼다.

첫 번째로 눈치챌 것: 카운트다. 정수. 음수 없음. 두 번째로: 분산이 평균에 비례한다. ACTB처럼 평균이 8,000이면 분산도 크고, TP53처럼 평균이 90이면 분산이 작다. 세 번째로: 평균과 분산의 관계가 정규분포가 예측하는 것보다 훨씬 넓다 (over-dispersion).

순진한 접근: 포아송

카운트 데이터니까 자연스럽게 포아송 분포를 떠올릴 수 있다. 포아송은 Var(X) = E(X) = μ이다. 그런데 실제 RNA-seq 데이터에서 반복 샘플간 분산은 평균보다 압도적으로 크다. 특히 저발현 유전자에서 그렇다. 포아송을 그냥 써버리면 통계량이 붕괴한다.

음이항 분포로 넘어가기

음이항 분포(Negative Binomial, NB)는 포아송에 여분의 흩어짐 파라미터 α (dispersion) 를 붙인 것이다.

Var(X) = μ + α · μ²

α가 0이면 포아송으로 수렴하고, α가 크면 훨씬 넓게 흩어진다. RNA-seq 실측 데이터로 (μ, Var) 산점도를 그리면 이 관계가 깜짝 놀랄 만큼 잘 맞는다. DESeq2는 여기서 시작한다.

Dispersion Shrinkage — DESeq2의 마법

여기서 실무의 벽을 만난다. 논문 실험이 보통 3~6개 반복만 있다. 소수의 샘플로 유전자 20,000개의 α를 각각 추정하려니 추정치가 무섭게 흔들린다. 어떤 유전자는 α = 0.01, 어떤 유전자는 α = 5.0 — 대부분이 표본 오차다.

DESeq2 저자 Mike Love는 여기서 Empirical Bayes shrinkage를 꺼낸다. 아이디어는 이렇다.

1. 유전자별 원시 α_g를 추정한다 (매우 흔들림).
2. 모든 유전자의 α를 평균-분산 그림에 얹어 매끄러운 fit line을 긋는다 (이게 "trend").
3. 유전자별 최종 α_g*는 원시 추정치와 trend 값의 가중 평균으로 정한다. 표본 수가 적을수록 trend 쪽으로 강하게 끌려간다.

이게 왜 옳은가? "비슷한 발현 수준의 유전자는 비슷한 정도로 흩어진다"는 생물학적 사전 확률(prior)을 도입한 것이다. 이 shrinkage 하나로 위양성이 극적으로 줄고, 소수 반복 실험의 재현성이 급상승한다.

StatQuest로 감 잡기

이 파트는 개념이 무거우니 시각적 도움이 절대 필요하다. StatQuest의 한국어 자막 강의로 감을 잡고 오면 좋다.

• StatQuest — A Gentle Introduction to RNA-seq & DESeq2 (Joshua Starmer, 26분, 한국어 자막 완벽 지원). YouTube 설정에서 자막을 한국어로 켜면 정확한 번역이 나온다. 카드 매칭 비유가 오래 남는다.

Wald test — 어떻게 유의성을 판정하나

각 유전자별로 두 조건(예: 유방암 조직 vs 정상 조직) 간 로그 발현 변화량 β를 추정하고, Wald test로 β ≠ 0 여부를 검정한다.

z = β̂ / SE(β̂)

z가 표준정규를 따른다고 가정하고 양측 p-값을 뽑는다. DESeq2는 log fold change 자체에도 shrinkage를 걸어 저발현 유전자의 극단적 β 추정치를 안정화한다. lfcShrink() 함수가 그 역할이다.

Benjamini-Hochberg — 20,000번 검정하는 죄

여기서 우리는 CS 다중 검정 문제와 정면충돌한다. 유전자 20,000개를 동시에 검정하면, p-값 임계 0.05로 잡아도 순전히 우연히 1,000개가 "유의"하게 나온다. 이걸 그대로 volcano plot에 찍으면 논문은 폐기된다.

Benjamini-Hochberg (BH) 조정은 이렇게 작동한다.

1. 모든 유전자의 p-값을 오름차순 정렬한다: p₍₁₎ ≤ p₍₂₎ ≤ ⋯ ≤ p₍ₘ₎.
2. FDR 목표(예: 0.05)를 정한다.
3. 다음 조건을 만족하는 최대 k를 찾는다:

p₍ₖ₎ ≤ (k / m) × FDR

4. 그 k까지의 유전자를 유의로 판정한다.

DESeq2가 반환하는 padj 컬럼이 BH 조정된 p-값이다. Volcano plot의 유의 기준선은 padj < 0.05가 표준이다.

BH의 이유가 아직 안 잡히면, StatQuest의 이 짧은 강의도 한국어 자막이 완벽하다.

• StatQuest — False Discovery Rates (FDR), clearly explained (Joshua Starmer, 18분, 한국어 자막 완벽 지원). 카드 놀이 비유로 왜 20,000번 검정하면 위양성이 폭발하는지 직관적으로 보여준다.

R 30줄로 TCGA 유방암 volcano plot 완성

이제 실습이다. Colab에서 R 커널을 켜자.

1단계: TCGA 데이터 가져오기

# Bioconductor 설치 (Colab 첫 실행 시 5분)
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("DESeq2", "TCGAbiolinks", "EnhancedVolcano"))

library(DESeq2)
library(TCGAbiolinks)

# TCGA-BRCA raw count 데이터 (종양 vs 정상 조직)
query <- GDCquery(
  project    = "TCGA-BRCA",
  data.category = "Transcriptome Profiling",
  data.type  = "Gene Expression Quantification",
  workflow.type = "STAR - Counts",
  sample.type = c("Primary Tumor", "Solid Tissue Normal")
)
GDCdownload(query)
data <- GDCprepare(query)

TCGA에서 유방암 환자와 정상 조직 raw count를 받는다. 처음 다운로드는 시간이 좀 걸리지만, cBioPortal이나 UCSC Xena를 쓰면 훨씬 빠르다.

2단계: DESeq2 실행

dds <- DESeqDataSet(data, design = ~ sample_type)
dds <- DESeq(dds)  # 여기서 dispersion 추정 + Wald test + BH 조정 다 일어난다

res <- results(dds, contrast = c("sample_type", "Primary Tumor", "Solid Tissue Normal"))
res <- lfcShrink(dds, coef = "sample_type_Primary Tumor_vs_Solid Tissue Normal", type = "apeglm")
summary(res)

summary(res)가 뱉는 "up-regulated / down-regulated / outlier / low-count" 요약이 딱 논문 methods에 그대로 들어가는 문장이다.

3단계: Volcano plot

library(EnhancedVolcano)
EnhancedVolcano(res,
    lab      = rownames(res),
    x        = "log2FoldChange",
    y        = "padj",
    pCutoff  = 0.05,
    FCcutoff = 1.0,
    title    = "TCGA-BRCA Tumor vs Normal"
)

이 세 줄이면 논문에 실을 수 있는 volcano plot이 나온다. 가로축 = log2 fold change (음수는 정상에서 높음, 양수는 종양에서 높음), 세로축 = -log10(padj). 오른쪽 위 · 왼쪽 위 코너에 찍히는 점들이 종양-특이적으로 발현이 크게 바뀐 후보 유전자다.

Wet-lab에서 이 그림을 보면 익숙한 이름이 나올 것이다. BRCA1, TP53, ERBB2, MKI67 …. 다음 편(S21)에서 이 유전자 리스트로 GSEA 경로 분석을 돌리면 우리가 잘 아는 세포주기·DNA 손상 반응 경로가 정확히 튀어나온다. 그때 진짜 논문의 모습이 완성된다.

왜 edgeR/limma가 아닌 DESeq2인가

• DESeq2: 저발현 유전자에서 shrinkage 강함. 소수 반복(3~6)에 강함. 학계 표준.
• edgeR: DESeq2와 형제. Tag-wise dispersion 방식이 약간 다름. bulk에서 거의 등치.
• limma-voom: log CPM으로 변환 후 정규분포 가정. 샘플 수가 많을 때 강함(30개+).

BPD 실습은 표본 수가 적어 DESeq2로 시작하지만, 다음 편(S20)에서 세 방식이 어떻게 갈라지는지 실측한다.

CS 매핑

• 다중 가설 검정 (Multiple hypothesis testing): DryBench에서 다룬 개념. BH는 그중 가장 실용적 조절법.
• Empirical Bayes: 유전자 20,000개의 dispersion을 각각 독립 추정하지 않고, 서로에게 정보를 빌려주는 방식. CS의 "regularization"과 정확히 같은 아이디어.
• Shrinkage estimator: James-Stein 추정량의 바이오 버전. 개별 추정치를 공통 목표로 끌어당긴다.

자주 만나는 결함

• 원시 카운트가 아닌 정규화된 값(TPM, RPKM)을 넣음 → DESeq2 실패. DESeqDataSet에는 반드시 raw count.
• padj가 아닌 pvalue로 유의성 판단 → 위양성 폭발.
• 저발현 유전자(모든 샘플 카운트 < 10)를 미제거 → 무작위 스파이크가 유의로 잡힘. filterByExpr() 사용.
• Batch effect 미제어 → 종양 vs 정상이 아니라 "시퀀싱 배치 A vs B"를 보고 있을 수 있음. design = ~ batch + sample_type으로 명시.

더 깊게 파고 싶다면

본문은 BPD가 직접 재구성한 서술이다. 심화는 아래 정통 자료를 순서대로.

• Harvard STAT115 — Xiaole Shirley Liu 교수의 DESeq2 library normalization & variance stabilization (52분, 자막 완비). Size factor 계산과 VST 변환 유도가 강하다.
• StatQuest — Joshua Starmer의 두 편 (둘 다 한국어 자막 완벽 지원):
  • A Gentle Introduction to RNA-seq & DESeq2 (26분)
  • False Discovery Rates (FDR), clearly explained (18분)
• 참고 무료 웹북: Holmes & Huber, Modern Statistics for Modern Biology — 이 편의 통계 원리 전체가 huber.embl.de/msmb에서 무료로 열려 있다. Chapter 8 "High-throughput count data" 정독을 권한다.
• 원 논문: Love, Huber, Anders (2014), Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biology 15:550. Shrinkage 유도의 원천.

다음 편 S20에서 edgeR / limma-voom과 실측 비교하고, S21에서 이 결과로 GSEA를 돌린다. 우리는 이제 "FASTQ to Paper" 여정의 정확히 절반을 지났다.