edgeR vs limma-voom vs DESeq2 — 언제 무엇을 쓰나, Empirical Bayes 완벽 이해 (한국어)

DESeq2를 배웠는데 왜 또 다른 도구를 배우나요

S19에서 우리는 DESeq2로 volcano plot을 완성했습니다. 그런데 논문을 투고하면 리뷰어가 이런 코멘트를 던지는 일이 잦습니다. "edgeR로도 재현되는지 확인해달라." 혹은 반대로, 샘플이 수십 개인 대규모 코호트를 받았는데 DESeq2가 너무 느려 밤을 새우게 됩니다.

세 도구 — DESeq2 · edgeR · limma-voom — 는 경쟁자가 아니라 같은 문제에 대한 세 개의 다른 통계 철학입니다. 이 편을 마치면 "어떤 상황에서 무엇을 꺼내야 하는가"를 손에 넣게 됩니다. 그리고 세 도구 모두의 심장에 Empirical Bayes라는 동일한 아이디어가 박혀 있다는 것을 확인하게 됩니다. S19에서 한 번, 이 편에서 두 번 더. 세 번 반복되면 그건 우연이 아니라 원리입니다.

edgeR — DESeq2의 형제

edgeR도 DESeq2와 똑같이 카운트를 음이항 분포로 모델링합니다. 갈라지는 곳은 두 군데입니다.

첫째, 정규화 방식입니다. DESeq2는 median-of-ratios size factor를 쓰고, edgeR은 TMM(Trimmed Mean of M-values)을 씁니다. TMM의 발상은 이렇습니다. 두 샘플의 로그 발현비(M값)와 평균 발현(A값)을 구한 뒤, 양극단을 잘라내고(trim) 남은 중앙부의 가중 평균으로 라이브러리 스케일을 잡습니다.

\log_2(\text{TMM factor}) = \frac{\sum_{g \in G^*} w_g M_g}{\sum_{g \in G^*} w_g}

여기서 $M_g = \log_2(y_{gk}/N_k) - \log_2(y_{gr}/N_r)$ 이고, $G^*$ 는 극단값을 잘라낸 유전자 집합입니다. 핵심 가정은 "대다수 유전자는 발현이 변하지 않는다"입니다. 소수의 초고발현 유전자(예: 글로빈)가 전체 라이브러리를 잡아먹어 다른 유전자를 인위적으로 낮게 보이게 만드는 composition bias를 이 trimming이 막아줍니다.

둘째, dispersion 추정입니다. edgeR은 common → trended → tagwise 세 단계로 α를 좁혀갑니다. tagwise 단계에서 유전자별 dispersion을 공통 추세 쪽으로 끌어당기는데, 이게 바로 DESeq2의 shrinkage와 같은 Empirical Bayes입니다. 이름만 다릅니다.

limma-voom — 완전히 다른 발상

limma는 원래 마이크로어레이(연속값) 시대의 도구였습니다. 정규분포와 선형모델로 20년간 다듬어진 통계 엔진입니다. RNA-seq 카운트는 정규분포가 아닌데, voom이 이 간극을 메웁니다.

voom의 아이디어는 놀랍도록 실용적입니다. 카운트를 log-CPM으로 변환하면 평균이 낮은 유전자일수록 분산이 큽니다(mean-variance trend). voom은 이 추세를 매끄러운 곡선으로 추정한 뒤, 각 관측치에 precision weight를 부여합니다.

w_{gi} = \frac{1}{s(\hat{\mu}_{gi})^2}

즉 "분산이 클 것으로 예측되는 관측치는 신뢰도를 낮춰라". 저발현 유전자의 불안정한 관측치가 회귀에 덜 기여하게 됩니다. 이 weight를 넣은 가중 최소제곱으로 선형모델을 적합하면, 카운트 데이터를 정규분포 세계의 강력한 도구로 그대로 끌고 올 수 있습니다.

moderated t — limma의 Empirical Bayes

여기서 세 번째 Empirical Bayes가 등장합니다. 일반 t-검정은 유전자별 표준오차 $s_g$ 를 각각 씁니다. 반복이 3개뿐이면 $s_g$ 가 요동칩니다. limma는 유전자별 분산 $s_g^2$ 를 전체 유전자의 사전 분산 $s_0^2$ 쪽으로 끌어당긴 값으로 t 통계량을 만듭니다.

\tilde{t}_g = \frac{\hat{\beta}_g}{\tilde{s}_g \sqrt{v_g}}, \qquad \tilde{s}_g^2 = \frac{d_0 s_0^2 + d_g s_g^2}{d_0 + d_g}

$\tilde{s}_g^2$ 는 사전 분산과 관측 분산의 가중 평균입니다. 반복이 적어 $d_g$ 가 작으면 $s_0^2$ 쪽으로 강하게 끌립니다. S19의 dispersion shrinkage와 완전히 같은 골격입니다. 대상이 dispersion이냐 residual variance냐만 다릅니다.

손으로 감 잡기

사전 자유도 $d_0 = 4$ , 사전 분산 $s_0^2 = 0.10$ 이라 하고, 어떤 유전자의 관측 분산이 $s_g^2 = 0.90$ (반복 3개라 $d_g = 2$ )라고 해봅시다.

\tilde{s}_g^2 = \frac{4 \times 0.10 + 2 \times 0.90}{4 + 2} = \frac{0.40 + 1.80}{6} = 0.367

관측 분산 0.90이 0.367로 끌려 내려왔습니다. 이 유전자는 "우연히 분산이 크게 나온" 케이스였을 가능성이 높고, moderated 분산이 그 과장을 눌러줍니다. 위양성이 줄어드는 원리가 여기 있습니다.

세 도구 언제 무엇을 쓰나

세 방식의 실무 선택 기준을 정리하면 다음과 같습니다.

도구	모델	강한 상황	속도
DESeq2	음이항 + LFC shrink	소수 반복(3～6), 저발현 강조	느림
edgeR	음이항 + QL F-test	소수 반복, 복잡한 설계식	중간
limma-voom	정규 + precision weight	샘플 많음(수십～수백), 대규모 코호트	빠름

경험칙 하나만 기억합시다. 반복이 적으면 DESeq2/edgeR, 샘플이 많으면 limma-voom. 그리고 세 도구가 상위 유의 유전자에서 크게 어긋나면, 그건 도구 문제가 아니라 데이터에 배치 효과나 이상치가 숨어 있다는 신호입니다.

Colab R로 세 방식 직접 비교

이제 실습입니다. airway 데이터셋(천식 약물 처리 vs 대조, 8샘플)으로 셋을 한자리에서 돌려봅시다.

1단계: 데이터와 패키지

if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("airway", "edgeR", "limma", "DESeq2"))

library(airway); data(airway)
counts <- assay(airway)                 # 유전자 × 샘플 raw count
group  <- airway$dex                    # trt vs untrt
counts <- counts[rowSums(counts) > 10, ] # 저발현 제거

2단계: edgeR

library(edgeR)
y <- DGEList(counts = counts, group = group)
y <- calcNormFactors(y)                 # TMM 정규화
design <- model.matrix(~ group)
y <- estimateDisp(y, design)            # common→trended→tagwise
fit <- glmQLFit(y, design)              # quasi-likelihood
qlf <- glmQLFTest(fit, coef = 2)
edger_res <- topTags(qlf, n = Inf)$table

3단계: limma-voom

library(limma)
v <- voom(y, design, plot = TRUE)       # precision weight 계산 (곡선 확인)
vfit <- lmFit(v, design)
vfit <- eBayes(vfit)                    # moderated t
voom_res <- topTable(vfit, coef = 2, number = Inf)

voom(plot = TRUE)가 그려주는 mean-variance 곡선을 꼭 눈으로 확인합시다. 이 곡선이 평평하면 이미 정규화가 잘 된 것이고, 심하게 휘어 있으면 voom weight가 큰 역할을 하고 있다는 뜻입니다.

4단계: 세 유의 유전자 집합 겹쳐 보기

sig <- function(tab, p, lfc) rownames(tab)[tab[[p]] < 0.05 & abs(tab[[lfc]]) > 1]
e <- sig(edger_res, "FDR", "logFC")
v_ <- sig(voom_res, "adj.P.Val", "logFC")
length(intersect(e, v_)) / length(union(e, v_))  # Jaccard 유사도

보통 0.8 이상이 나옵니다. 세 도구가 대체로 같은 결론을 주되, 경계선 근처 유전자에서만 갈립니다. 그 경계 유전자가 바로 후속 검증(qPCR)의 우선순위입니다.

CS 매핑

Empirical Bayes = regularization: 개별 추정치를 공통 목표로 끌어당기는 것은 CS의 정규화(L2, ridge)와 같은 골격입니다. S19의 dispersion, 이 편의 residual variance — 대상만 다를 뿐 정확히 같은 수학입니다.
precision weight = 가중 최소제곱(WLS): voom은 관측치별 신뢰도를 역분산 가중으로 넣습니다. 통계·수치해석에서 이질 분산(heteroscedasticity)을 다루는 표준 기법 그대로입니다.
TMM trimming = robust statistics: 극단값을 잘라내고 중앙부만 쓰는 것은 trimmed mean, median 같은 강건 통계의 발상입니다.

자주 만나는 결함

정규화 팩터를 안 걸고 바로 검정 → composition bias로 결과가 통째로 왜곡됩니다. calcNormFactors()는 선택이 아닙니다.
limma-voom에 이미 log/CPM 된 값을 다시 넣음 → voom은 raw count(DGEList)를 기대합니다. 이중 변환은 곡선을 망칩니다.
샘플 3개인데 limma-voom 고집 → 정규 근사가 약해집니다. 소수 반복은 DESeq2/edgeR이 안전합니다.
세 도구 결과가 크게 다른데 "도구 탓" → 대개 숨은 배치 효과입니다. design = ~ batch + group으로 먼저 잡읍시다.

더 깊게 파고 싶다면

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

Harvard STAT115 — Xiaole Shirley Liu 교수의 Differential Expression: edgeR & limma (자막 완비, 자동번역 우수). 세 도구의 모델 차이를 통계적으로 정리합니다.
원 논문(voom): Law, Chen, Shi, Smyth (2014), voom: precision weights unlock linear model analysis tools for RNA-seq read counts, Genome Biology 15:R29. precision weight 유도의 원천입니다.
원 논문(edgeR): Robinson, McCarthy, Smyth (2010), edgeR: a Bioconductor package for differential expression analysis of digital gene expression data, Bioinformatics 26:139.
원 논문(moderated t): Smyth (2004), Linear models and empirical Bayes methods for assessing differential expression in microarray experiments, SAGMB 3:Article 3. Empirical Bayes의 고전입니다.
공식 문서: limma User's Guide (bioconductor.org/packages/limma) — RNA-seq 챕터가 실무 레시피의 정석입니다.

세 도구를 손에 넣었으니, 다음 편 S21에서 이 유의 유전자 리스트를 GSEA에 넣어 "개별 유전자"에서 "생물학적 경로"로 시야를 한 단계 넓힙니다. 그때 "FASTQ to Paper"의 마지막 조각이 맞춰집니다.

BioPlayground

edgeR와 limma-voom: DESeq2 말고 다른 길, 그리고 Empirical Bayes가 세 번 반복되는 이유