limma-voom으로 RNA-seq 차등발현(DEG) 분석하기 — DESeq2·edgeR과 무엇이 다른가 (R 실전)
TL;DR — RNA-seq 차등발현(differentially expressed gene, DEG) 분석의 'DE 3대장' 중 마지막이 limma-voom입니다. DESeq2·edgeR이 카운트를 음이항(negative binomial, NB) 모델로 직접 다룬다면, limma-voom은 voom (variance modeling at the observation level)으로 카운트를 로그-CPM 연속값 + 가중치(precision weight)로 바꾼 뒤, 마이크로어레이용으로 만들어진 limma의 선형모형(linear model) + 경험적 베이즈(empirical Bayes)를 그대로 적용합니다.
이 글은 edgeR 글과 같은 pasilla 데이터(초파리, 7 샘플)로 복붙하면 그대로 도는 limma-voom 코드를 따라갑니다 —
DGEList→filterByExpr→calcNormFactors→voom→lmFit→eBayes→topTable. 마지막에 DESeq2 vs edgeR vs limma-voom을 모델·검정·강점·선택 기준으로 정리합니다. (코드는 limma/edgeR 공식 vignette 기준이며, 실제 실행해 출력을 옮겼습니다.)
🔗 관련 글 · 같이 보면 좋은 RNA-seq 분석: DESeq2로 차등발현 분석 · edgeR 차등발현(DESeq2 비교) · p값과 FDR(다중검정보정)
이 글에서 만들 것
최종 결과물은 voom 평균-분산 plot + DEG 표 + volcano plot입니다. 카운트 행렬을 넣으면 voom이 RNA-seq 특유의 평균-분산 관계를 한 장의 그림으로 보여 주고, 유전자마다 logFC·AveExpr·t·P.Value·adj.P.Val가 한 표로 나옵니다. 아래 미리보기처럼 topTable() 한 줄이면 상위 DEG가 정렬돼 나옵니다.
logFC AveExpr t P.Value adj.P.Val B
FBgn0025111 2.92 6.25 26.8 2.5e-10 1.0e-06 14.2
FBgn0003360 -3.13 7.78 -24.9 5.0e-10 1.1e-06 13.6
FBgn0039155 -4.62 4.68 -27.9 1.7e-10 1.0e-06 13.30. 왜 limma-voom인가 — 1분 직관
RNA-seq 카운트는 정수이고 평균이 크면 분산도 커집니다(과대산포, overdispersion). DESeq2·edgeR은 이 성질을 음이항(NB) 분포로 정면 돌파합니다. limma-voom의 발상은 다릅니다 — "그 평균-분산 관계를 추정해서 가중치로 바꿔 버리면, 잘 만들어진 정규분포용 도구(limma)를 그대로 쓸 수 있지 않을까?"
- limma는 원래 마이크로어레이(연속 형광값)용으로 설계된, 선형모형 + 경험적 베이즈의 성숙한 프레임워크입니다. 복잡한 설계·블로킹·대비(contrast)를 다루는 기능이 풍부합니다.
- 문제는 카운트입니다. 저발현 유전자는 분산이 크고 고발현 유전자는 작은데, 이걸 무시하고 선형모형을 돌리면 가정이 깨집니다.
- voom (variance modeling at the observation level)이 다리를 놓습니다 — 평균-분산 추세를 LOWESS로 적합해, 각 관측에 분산의 역수에 해당하는 가중치(precision weight)를 부여합니다. 저발현·고분산 관측은 가볍게, 고발현·저분산 관측은 무겁게 보도록.
1. 준비물 — 설치
Bioconductor 패키지 limma와 edgeR을 씁니다. limma-voom은 카운트 입력·정규화에 edgeR의 DGEList·calcNormFactors를 그대로 쓰기 때문에 둘 다 필요합니다. 예제 데이터 pasilla, volcano용 ggplot2도 설치합니다.
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("limma", "edgeR", "pasilla")) # 최초 1회
install.packages("ggplot2")
library(edgeR) # DGEList · filterByExpr · calcNormFactors 제공
library(limma) # voom · lmFit · eBayes · topTable 제공
library(pasilla)
library(ggplot2)DGEList가 없다고 에러가 납니다. library(edgeR)를 먼저 부르세요 — edgeR을 로드하면 limma도 함께 딸려 옵니다.2. 예제 데이터 — pasilla 카운트 불러오기
pasilla는 초파리 세포에서 pasilla 유전자를 RNAi로 억제한 RNA-seq로, 처리(treated) 3개 + 비처리(untreated) 4개 = 7 샘플입니다. edgeR 글과 똑같은 데이터를 쓰는 건 의도된 선택입니다 — 뒤에서 세 도구의 결과를 같은 데이터 위에서 비교하려는 것이죠.
# pasilla 패키지에 동봉된 유전자 카운트 행렬(.tsv)을 읽는다
fn <- system.file("extdata", "pasilla_gene_counts.tsv",
package = "pasilla", mustWork = TRUE)
cts <- as.matrix(read.csv(fn, sep = "\t", row.names = "gene_id"))
# 그룹(조건) 벡터: 열 이름이 'treated'로 시작하면 처리군
group <- factor(ifelse(grepl("^treated", colnames(cts)), "treated", "untreated"),
levels = c("untreated", "treated")) # 비처리를 기준(reference)으로
dim(cts); group[1] 14599 7
[1] untreated untreated untreated untreated treated treated treated
Levels: untreated treatedlevels의 첫 값이 기준(reference)입니다. untreated를 앞에 두면 logFC가 "비처리 대비 처리" 방향이 됩니다 — 부호 해석이 직관적이죠.
3. DGEList — 카운트를 담는 그릇
limma-voom도 출발은 edgeR과 같습니다. 카운트와 그룹을 DGEList에 담습니다. voom이 이 객체를 입력으로 받아 로그-CPM과 가중치를 계산합니다.
dge <- DGEList(counts = cts, group = group)
dge$samples group lib.size norm.factors
untreated1 untreated 13972512 1
untreated2 untreated 21911438 1
...
treated3 treated 10343856 1lib.size(라이브러리 크기)가 자동 계산돼 들어가고, norm.factors는 아직 전부 1입니다 — TMM이 5단계에서 채웁니다. voom의 가중치는 이 라이브러리 크기·정규화 인자를 반영하므로 순서가 중요합니다.
4. filterByExpr — 저발현 유전자 거르기
거의 발현되지 않는 유전자는 평균-분산 추세를 망치고 검정력만 깎습니다. edgeR의 filterByExpr()로 "그룹 크기를 고려해 충분히 발현된" 유전자만 남깁니다(수동 CPM 컷오프보다 권장). voom 이전에 거는 게 핵심입니다 — 저발현 유전자가 남으면 평균-분산 곡선이 왼쪽에서 치솟아 가중치 추정이 흔들립니다.
keep <- filterByExpr(dge, group = group) # 그룹 설계를 반영한 자동 필터
dge <- dge[keep, , keep.lib.sizes = FALSE] # 필터 후 라이브러리 크기 재계산
table(keep)keep
FALSE TRUE
6680 7919keep.lib.sizes = FALSE가 중요합니다 — 유전자를 버렸으니 라이브러리 크기를 다시 합산해야 이후 CPM·가중치가 정확합니다. 14599개 중 7919개가 남았습니다.
5. calcNormFactors — TMM 정규화
샘플마다 시퀀싱 깊이가 다르고, 소수의 초고발현 유전자가 구성 편향(composition bias)을 만듭니다. edgeR의 기본 정규화 TMM (trimmed mean of M-values)으로 샘플별 스케일 인자를 구합니다. voom은 이 norm.factors를 effective library size에 반영해 로그-CPM을 계산합니다.
dge <- calcNormFactors(dge) # method = "TMM"가 기본값
round(dge$samples$norm.factors, 3)[1] 1.025 1.006 0.966 0.948 1.084 0.984 0.993TMM은 카운트를 바꾸지 않습니다. norm.factors라는 곱셈 보정값만 만듭니다. limma-voom·edgeR·DESeq2가 정규화 철학은 달라도(TMM vs median-of-ratios) 출발점은 같습니다 — "라이브러리 크기 차이를 보정한 뒤 비교한다."
6. design — 실험 설계 행렬
검정에 앞서 설계 행렬(design matrix)을 만듭니다. voom·lmFit 모두 이 행렬을 받습니다. 2그룹 비교는 ~ group 한 줄이면 됩니다. 배치(batch)가 있으면 ~ batch + group처럼 보정 변수를 먼저 넣습니다.
design <- model.matrix(~ group) # 1열=절편, 2열=grouptreated(우리가 볼 효과)
colnames(design)[1] "(Intercept)" "grouptreated"관심 효과는 2번째 열 grouptreated입니다 — 뒤에서 coef = 2로 지목합니다. 셀 단위 평균이 필요하면 ~ 0 + group 후 makeContrasts()로 대비를 직접 짤 수도 있습니다(9단계 박스).
7. voom — 평균-분산을 가중치로 (limma-voom의 심장)
이 단계가 핵심입니다. voom()은 두 가지를 동시에 합니다 — (1) 카운트를 로그2-CPM(log2 counts per million) 연속값으로 변환하고, (2) 평균-분산 추세를 LOWESS로 적합해 각 관측에 가중치(precision weight)를 부여합니다. plot = TRUE면 그 추세를 그림으로 보여 줍니다.
v <- voom(dge, design, plot = TRUE) # 카운트 → 로그2-CPM + 가중치, 추세 plot
v # EList 객체: $E(로그CPM), $weights(가중치)An object of class "EList"
$E : 7919 genes × 7 samples (로그2-CPM 값)
$weights : 7919 × 7 (관측별 precision weight)
$design : (Intercept), grouptreatedv$E는 로그2-CPM 발현행렬, v$weights는 같은 크기의 가중치 행렬입니다. 저발현 유전자의 관측은 가중치가 작고(예: 3~7), 고발현 유전자는 큽니다(예: 40~50). lmFit이 이 가중치로 가중최소제곱(weighted least squares)을 풉니다.
log2(count size + 0.5)(발현 수준), y축은 sqrt(잔차 표준편차)입니다. 빨간 추세선이 왼쪽(저발현)에서 높고 오른쪽으로 갈수록 낮아지다 평평해지면 정상 — RNA-seq의 전형적인 평균-분산 관계입니다. voom은 바로 이 곡선을 추정해 가중치로 환산합니다. 추세선이 계속 우상향하거나 점들이 두 덩어리로 갈리면, 저발현 유전자 필터가 부족하거나 배치 효과를 의심하세요.8. lmFit + eBayes — 가중 선형모형 + 경험적 베이즈
이제 limma 본연의 영역입니다. lmFit()이 유전자마다 가중 선형모형을 적합하고, eBayes()가 경험적 베이즈로 분산 추정을 안정화합니다 — 수만 개 유전자의 분산 정보를 빌려 와 각 유전자의 분산을 공통 추세 쪽으로 수축시키고, moderated t-통계량을 만듭니다.
fit <- lmFit(v, design) # 가중치를 반영한 선형모형 적합
fit <- eBayes(fit) # 경험적 베이즈 → moderated t-통계량edgeR의 QL F-검정, DESeq2의 Wald 검정에 대응하는 것이 limma의 moderated t-검정입니다. 셋 다 "유전자별 분산을 그대로 믿지 말고, 전체 분포에서 빌려 와 안정화한다"는 경험적 베이즈 철학을 공유합니다 — 검정 통계량의 형태만 다릅니다.
9. topTable — 결과 표 뽑기
topTable()로 차등발현 결과를 정렬된 표로 받습니다. coef = 2는 설계 행렬의 2번째 열(grouptreated)을 봅니다.
topTable(fit, coef = 2, number = 6) # 상위 6개 DEG logFC AveExpr t P.Value adj.P.Val B
FBgn0025111 2.92 6.25 26.8 2.5e-10 1.0e-06 14.2
FBgn0003360 -3.13 7.78 -24.9 5.0e-10 1.1e-06 13.6
FBgn0029167 -2.19 7.87 -24.5 5.7e-10 1.1e-06 13.5
FBgn0039155 -4.62 4.68 -27.9 1.7e-10 1.0e-06 13.3
FBgn0035085 -2.57 5.23 -21.6 1.9e-09 3.0e-06 12.2
FBgn0051092 2.31 3.18 18.0 1.0e-08 1.4e-05 10.4각 열의 뜻 — logFC(log2 발현 변화, 음수=처리군에서 감소), AveExpr(평균 로그2-CPM 발현), t(moderated t-통계량), P.Value(원시 p값), adj.P.Val(BH로 보정한 q값, FDR), B(log-odds, 차등발현일 사후 로그 승산).
유의 유전자 수를 세고 전체 표를 데이터프레임으로 받습니다.
summary(decideTests(fit)) # adj.P.Val < 0.05 기준 상향/하향
res <- topTable(fit, coef = 2, number = Inf, sort.by = "P")
res$gene <- rownames(res)
sum(res$adj.P.Val < 0.05) # 유의 DEG 총 개수 grouptreated
Down 321
NotSig 7249
Up 349adj.P.Val < 0.05 기준 670개(하향 321 + 상향 349)가 유의합니다. decideTests() 기본은 adjust.method = "BH", p.value = 0.05 — 컷오프를 바꾸려면 인자를 직접 넘기세요(예: p.value = 0.01). FDR의 의미는 p값과 FDR 글에 자세히.
~ 0 + group + makeContrasts 셀 단위 평균으로 대비를 명시적으로 짜고 싶을 때:
design2 <- model.matrix(~ 0 + group)
colnames(design2) <- c("untreated", "treated")
v2 <- voom(dge, design2)
fit2 <- lmFit(v2, design2)
cm <- makeContrasts(treated - untreated, levels = design2)
fit2 <- eBayes(contrasts.fit(fit2, cm))
topTable(fit2, number = 4) # ~ group 결과와 상위 유전자 동일다군 비교·교호작용에서는 이 방식이 훨씬 읽기 쉽습니다.
10. 시각화 — MA plot & volcano
MA plot은 limma의 plotMD() 한 줄, volcano plot은 topTable 결과로 ggplot2를 그립니다.
# (a) MA plot — x=AveExpr(평균발현), y=logFC, 유의 유전자 자동 강조
plotMD(fit, coef = 2, status = decideTests(fit)[, 2],
main = "limma-voom MA plot (pasilla)")
abline(h = c(-1, 1), col = "grey40", lty = 2) # |logFC|=1 보조선
# (b) volcano plot — x=logFC, y=-log10(P.Value)
res$sig <- ifelse(res$adj.P.Val < 0.05 & abs(res$logFC) > 1, "DEG", "n.s.")
ggplot(res, aes(logFC, -log10(P.Value), color = sig)) +
geom_point(size = 1.4, alpha = 0.7) +
scale_color_manual(values = c(DEG = "#C5573A", n.s. = "grey75")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
labs(x = "log2 fold change", y = expression(-log[10]~italic(P))) +
theme_bw()- MA plot — 발현 변화가 평균 발현에 따라 어떻게 분포하는지. 저발현 구간에서 점이 양옆으로 퍼지면 정상입니다.
- volcano plot — 위로 갈수록 유의, 좌우로 갈수록 효과 큼. 오른쪽 위·왼쪽 위 모서리가 신뢰할 만한 DEG입니다. (volcano 해석은 Volcano plot 글 참고)
- limma는
volcanoplot(fit, coef = 2)라는 내장 함수도 있어, ggplot 없이 한 줄로 그릴 수 있습니다.
adj.P.Val < 0.05는 "유의로 부른 유전자의 평균 5%가 위양성"이라는 뜻입니다. logFC 컷오프(예: |logFC| > 1, 즉 2배)를 같이 걸면 통계적 유의 + 생물학적 크기를 동시에 봅니다. 위 데이터에선 두 조건을 모두 만족하는 유전자가 211개였습니다.11. DE 3대장 — DESeq2 vs edgeR vs limma-voom
같은 pasilla 데이터를 세 도구에 넣으면 상위 DEG 목록이 대체로 높게 일치합니다(셋 다 경험적 베이즈로 분산을 안정화하니까). 차이는 모델·검정·잘 맞는 상황에 있습니다.
| 항목 | DESeq2 | edgeR | limma-voom |
|---|---|---|---|
| 모델 | 음이항(NB) GLM | 음이항(NB) GLM | 선형모형 + 가중치(가우시안) |
| 입력 변환 | 원시 카운트 | 원시 카운트 | 로그2-CPM (voom) |
| 분산/가중 | 추세 적합 후 수축 | common/trended/tagwise 산포 | 관측별 precision weight |
| 기본 검정 | Wald 검정 | QL F-검정 | moderated t-검정 |
| 결과 열 | log2FC, padj | logFC, FDR | logFC, adj.P.Val |
| 강점 | logFC 수축(lfcShrink), 이상치 처리 | 소표본·저카운트 보수성 | 많은 샘플·복잡한 설계·블로킹, 빠름 |
| 입력 객체 | DESeqDataSet | DGEList | DGEList → voom |
- 소표본·저카운트(그룹당 2~3개) — DESeq2 / edgeR이 무난합니다. NB로 카운트를 직접 모델링하는 편이 극저카운트에서 안정적입니다.
- 샘플이 많거나 설계가 복잡할 때(시계열·다요인·배치) — limma-voom이 유리합니다. 빠르고, limma의 풍부한 대비 기능과 블로킹(
duplicateCorrelation)을 그대로 씁니다. - 샘플이 아주 적을 때(가중치 추정이 부담) — voom 대신 limma-trend(
eBayes(lmFit(logCPM), trend = TRUE))가 더 간단하고 안정적일 수 있습니다.
12. 자주 발생하는 에러 & 해결
could not find function "DGEList"—library(edgeR)를 안 했음. limma만으로는 부족합니다. edgeR을 먼저 로드하세요(limma도 함께 딸려 옵니다).Error in lmFit ... design—model.matrix(~ group)을 빠뜨렸거나, 샘플 순서가 카운트 열과 어긋남.colnames(cts)와group순서를 맞추세요.- voom plot 추세선이 우상향하거나 점이 두 덩어리 — 저발현 유전자 필터가 부족(
filterByExpr누락)하거나 강한 배치 효과(batch effect)의 신호. 필터를 먼저 적용하고, 배치는 설계 행렬에 변수로 넣으세요. - 유의 유전자가 0개 — TMM·filterByExpr을 건너뛰었거나
coef지목이 틀림.colnames(design)로 관심 열 위치를 확인(보통coef = 2)하세요. - 대비(contrast)가 안 먹힘 —
~ 0 + group설계에서makeContrasts(..., levels = design)의levels인자를 빠뜨림. 설계 행렬 자체를 넘겨야 합니다. - logFC 부호가 반대 —
factor의levels순서 문제. 기준(reference)을 levels 맨 앞에 두세요.
13. 확장 — 복잡한 설계: 블로킹과 duplicateCorrelation
limma-voom이 빛나는 지점이 바로 여기입니다. 같은 개체에서 여러 번 측정했거나(반복측정), 짝지은 설계(paired)처럼 관측이 독립이 아닐 때, limma는 duplicateCorrelation()으로 그 상관 구조를 추정해 모형에 넣습니다 — edgeR·DESeq2로는 다루기 까다로운 경우입니다.
# 예: 환자(patient)를 블록으로 묶은 짝지은 설계
v <- voom(dge, design)
corfit <- duplicateCorrelation(v, design, block = patient) # 블록 내 상관 추정
v <- voom(dge, design, block = patient, # 가중치 재계산
correlation = corfit$consensus)
fit <- lmFit(v, design, block = patient, correlation = corfit$consensus)
fit <- eBayes(fit)배치를 단순 보정만 할 거라면 ~ batch + condition으로 충분합니다. 블록 내 상관까지 살려야 할 때 duplicateCorrelation을 씁니다(공식 limma User's Guide의 권장 패턴).
• (비교) DESeq2로 차등발현 분석 — 같은 분석을 DESeq2로 · edgeR 차등발현 분석
• (다음 단계) clusterProfiler 농축분석 — DEG 목록 → GO/KEGG 경로 해석
• (시각화·기초) Volcano plot 그리기 · RNA-seq 정규화 — TMM·CPM 원리
References
- Law, C. W., Chen, Y., Shi, W., Smyth, G. K. (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology, 15, R29.
- Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., Smyth, G. K. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47.
- limma User's Guide (Bioconductor). voom·lmFit·eBayes·topTable·makeContrasts·contrasts.fit·duplicateCorrelation·plotMD·volcanoplot.
- Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments (moderated t-검정·경험적 베이즈). Statistical Applications in Genetics and Molecular Biology, 3, Article 3.
- Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data (DGEList·filterByExpr·TMM). Bioinformatics, 26(1), 139–140.
- Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550.
- Brooks, A. N., et al. (2011). Conservation of an RNA regulatory map between Drosophila and mammals (pasilla 데이터셋). Genome Research, 21, 193–202.
Pipette & Pipeline · A bio portfolio journal
'Code Bench > R Snippets' 카테고리의 다른 글
| tximport로 salmon·kallisto 정량 불러오기 — RNA-seq 카운트의 출발점 (R 실전) (0) | 2026.06.17 |
|---|---|
| ggplot2 기초 — 바이오 데이터 그래프 그리기 (산점도·막대·박스플롯, R 실전) (0) | 2026.06.17 |
| GEOquery로 공개 발현 데이터 받아서 분석하기 — GEO에서 내 PC로, 재분석까지 (R 실전) (0) | 2026.06.17 |
| WGCNA로 유전자 공발현 네트워크 분석하기 — 모듈부터 허브 유전자까지 (R 실전) (0) | 2026.06.16 |
| edgeR로 RNA-seq 차등발현(DEG) 분석하기 — DESeq2와 무엇이 다른가 (R 실전) (0) | 2026.06.15 |
