edgeR로 RNA-seq 차등발현(DEG) 분석하기 — DESeq2와 무엇이 다른가 (R 실전)

TL;DR — RNA-seq 차등발현(differentially expressed gene, DEG) 분석의 양대 표준이 edgeR과 DESeq2입니다. 둘 다 음이항(negative binomial, NB) 모델 + 경험적 베이즈(empirical Bayes)라는 같은 뼈대를 쓰지만, 정규화·검정 방식이 다릅니다. edgeR은 TMM (trimmed mean of M-values) 정규화 + 준가능도(quasi-likelihood, QL) F-검정을 씁니다.

이 글은 공개 데이터 pasilla (초파리)로 복붙하면 그대로 도는 edgeR 코드를 따라갑니다 — DGEList → filterByExpr → calcNormFactors(TMM) → estimateDisp → glmQLFit / glmQLFTest → topTags. 마지막에 edgeR vs DESeq2를 정규화·모델·검정·선택 기준으로 정리합니다. (코드는 edgeR 4.x 공식 User's Guide 기준입니다.)

🔗 관련 글  ·  같이 보면 좋은 R 분석: DESeq2로 차등발현 분석  ·  p값과 FDR(다중검정보정)  ·  clusterProfiler 농축분석

그림 1. counts → DGEList → TMM → estimateDisp(NB) → glmQLFTest → topTags(FDR)로 이어지는 edgeR 파이프라인.

0. 왜 edgeR인가 — 1분 직관

RNA-seq 카운트는 정수이고, 평균이 크면 분산도 커지는(과대산포, overdispersion) 성질이 있습니다. 그래서 정규분포 t-검정이 아니라 음이항(NB) 분포로 모델링합니다 — 이 지점은 edgeR과 DESeq2가 똑같습니다.

• 문제는 샘플 수가 적다는 것입니다(보통 그룹당 2~3개). 유전자마다 분산을 따로 추정하면 너무 불안정합니다.
• edgeR의 해법은 경험적 베이즈(empirical Bayes) — 수만 개 유전자의 분산 정보를 빌려 와 각 유전자의 분산 추정을 안정화합니다(common → trended → tagwise 분산).
• 여기에 준가능도(QL) F-검정을 얹어, 분산 추정의 불확실성까지 검정에 반영합니다. 이게 edgeR의 현재 권장 경로입니다.

📌 한 줄 요약 — NB로 분포를 잡고, 경험적 베이즈로 분산을 빌려 오고, QL F-검정으로 보수적으로 본다.

1. 준비물 — 설치

Bioconductor 패키지 edgeR과 예제 데이터 pasilla를 씁니다. volcano용으로 ggplot2도 설치합니다.

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("edgeR", "pasilla"))   # 최초 1회
install.packages("ggplot2")

library(edgeR)
library(pasilla)
library(ggplot2)

2. 예제 데이터 — pasilla 카운트 불러오기

pasilla는 초파리 세포에서 pasilla 유전자를 RNAi로 억제한 RNA-seq로, 처리(treated) 3개 + 비처리(untreated) 4개 = 7 샘플입니다(DESeq2·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)으로
group

[1] untreated untreated untreated untreated treated treated treated
Levels: untreated treated

levels의 첫 값이 기준(reference)입니다. untreated를 앞에 두면 logFC가 "비처리 대비 처리" 방향이 됩니다 — 부호 해석이 직관적이죠.

3. DGEList — edgeR의 기본 자료구조

edgeR의 모든 함수는 DGEList 객체를 입력으로 받습니다. 카운트 행렬과 그룹만 넣으면 됩니다(DESeq2의 DESeqDataSet에 해당).

y <- DGEList(counts = cts, group = group)
y

An object of class "DGEList"
$counts          : 14599 genes × 7 samples
$samples         : group, lib.size, norm.factors

y$samples에 라이브러리 크기(lib.size)가 자동 계산돼 들어갑니다. norm.factors는 아직 전부 1 — 다음 단계에서 TMM으로 채웁니다.

4. filterByExpr — 저발현 유전자 거르기

거의 발현되지 않는 유전자는 검정력만 깎고 다중검정 보정 부담만 키웁니다. edgeR은 filterByExpr()로 "그룹 크기를 고려해 충분히 발현된" 유전자만 남깁니다(수동 CPM 컷오프보다 권장).

keep <- filterByExpr(y, group = group)         # 그룹 설계를 반영한 자동 필터
y <- y[keep, , keep.lib.sizes = FALSE]         # 필터 후 라이브러리 크기 재계산
table(keep)

keep
FALSE  TRUE
 6500  8099

keep.lib.sizes = FALSE가 중요합니다 — 유전자를 버렸으니 라이브러리 크기를 다시 합산해야 정규화가 정확합니다. 빠뜨리면 옛 lib.size가 남습니다.

5. calcNormFactors — TMM 정규화

샘플마다 시퀀싱 깊이도 다르고, 소수의 초고발현 유전자가 구성 편향(composition bias)을 일으킵니다. edgeR의 기본 정규화가 TMM (trimmed mean of M-values) — 두 샘플의 로그비(M값)를 양극단을 잘라낸 뒤 평균 내, 그 샘플의 스케일 인자를 구합니다.

y <- calcNormFactors(y)        # method = "TMM"가 기본값
y$samples

            group    lib.size norm.factors
untreated1  untreated  ...        1.012
treated1    treated    ...        0.951   # 1보다 작으면 상향 보정
...

• TMM은 카운트를 바꾸지 않습니다. norm.factors라는 곱셈 보정값을 만들어, 모델이 effective library size로 쓰게 합니다.
• edgeR 4.x에서는 같은 함수가 normLibSizes()라는 이름으로도 제공됩니다(동작 동일). 기존 코드의 calcNormFactors()도 그대로 작동합니다.

⚠️ 핵심 — TMM은 "대부분의 유전자는 차등발현하지 않는다"는 가정 위에서 동작합니다. 전사체 절반이 통째로 바뀌는 극단적 설계에서는 spike-in 등 다른 정규화를 고려하세요.

6. design — 실험 설계 행렬

검정에 앞서 설계 행렬(design matrix)을 만듭니다. 2그룹 비교는 ~ group 한 줄이면 됩니다. 배치(batch)가 있으면 ~ batch + group처럼 보정 변수를 먼저 넣습니다.

design <- model.matrix(~ group)    # 1열=절편, 2열=grouptreated(우리가 볼 효과)
colnames(design)

[1] "(Intercept)"   "grouptreated"

관심 효과는 2번째 열 grouptreated입니다 — 뒤에서 coef = 2로 지목합니다. 셀 단위 평균이 필요하면 ~ 0 + group 후 makeContrasts()로 대비(contrast)를 직접 짤 수도 있습니다.

7. estimateDisp — 산포(dispersion) 추정

NB 모델의 핵심 모수가 산포(dispersion) — "평균 대비 분산이 얼마나 더 퍼지는가"입니다. estimateDisp() 한 번이 common·trended·tagwise 세 가지를 모두 추정합니다(경험적 베이즈로 유전자별 추정을 공통 추세 쪽으로 수축).

y <- estimateDisp(y, design)       # common → trended → tagwise 한 번에
sqrt(y$common.dispersion)          # BCV(생물학적 변동계수) = sqrt(분산)

[1] 0.31   # BCV ≈ 0.31 → 사람 데이터(0.4 안팎)보다 낮은, 깔끔한 셀라인

• BCV (biological coefficient of variation) = sqrt(dispersion). 세포주 0.1~0.2, 사람 임상 0.4 이상이 흔합니다.
• edgeR 4.x에서는 glmQLFit()이 산포를 직접 추정할 수 있어 이 단계를 건너뛸 수도 있습니다. 다만 plotBCV() 진단을 보려면 estimateDisp()를 명시적으로 두는 편이 학습에 좋습니다.

8. glmQLFit + glmQLFTest — 준가능도 F-검정 (권장)

이제 유전자마다 NB 일반화선형모형(GLM)을 적합하고, QL F-검정으로 차등발현을 검정합니다. 이게 edgeR의 현재 1순위 권장 경로입니다 — Wald/LRT보다 위양성(false positive)을 더 잘 통제합니다.

fit <- glmQLFit(y, design)             # 준가능도 NB GLM 적합
qlf <- glmQLFTest(fit, coef = 2)       # 2번째 열(grouptreated) 검정
topTags(qlf, n = 5)                    # 상위 5개 DEG

        logFC   logCPM      F       PValue          FDR
FBgn..  -4.05    5.49    234.6   1.5e-12       3.1e-09
FBgn..   2.71    7.83    180.2   4.4e-11       4.5e-08
...

각 열의 뜻 — logFC (log2 발현 변화, 음수=처리군에서 감소), logCPM (평균 발현 수준), F (QL F-통계량), PValue (원시 p값), FDR (BH 다중검정보정된 q값).

🔁 대안 경로 두 가지
• classic exactTest — 2그룹·공변량 없을 때만. et <- exactTest(y) (설계 행렬 불필요). 가장 단순하지만 배치 보정 불가.
• glmLRT — glmFit() 후 glmLRT(fit, coef = 2). 우도비(likelihood ratio) 검정. 표본이 많을 때 무난하나, 소표본에선 QL F-검정이 더 안전합니다.

9. 결과 정리 — 유의 DEG 세기 & 표 추출

decideTests()로 FDR 컷오프 기준 상향/하향 DEG 수를 세고, topTags()로 전체 표를 데이터프레임으로 받습니다.

summary(decideTests(qlf))              # FDR < 0.05 기준 유의 유전자 수
res <- topTags(qlf, n = Inf)$table     # 전체 결과를 데이터프레임으로
res$gene <- rownames(res)

sum(res$FDR < 0.05)                    # 유의 DEG 총 개수
head(res[res$FDR < 0.05, ])            # 유의 유전자만 미리보기

       grouptreated
Down            512
NotSig         7100
Up              487

decideTests()의 기본은 adjust.method = "BH", p.value = 0.05입니다. 컷오프를 바꾸려면 인자를 직접 넘기세요(예: p.value = 0.01).

10. 시각화 — MA plot & volcano

MA plot (평균-차이)은 plotMD() 한 줄, volcano plot은 topTags 결과로 ggplot2를 그립니다.

# (a) MA plot — x=logCPM, y=logFC, 유의 유전자 자동 강조
plotMD(qlf, main = "edgeR MA plot (pasilla)")
abline(h = c(-1, 1), col = "grey40", lty = 2)   # |logFC|=1 보조선

# (b) volcano plot — x=logFC, y=-log10(PValue)
res$sig <- ifelse(res$FDR < 0.05 & abs(res$logFC) > 1, "DEG", "n.s.")
ggplot(res, aes(logFC, -log10(PValue), 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()

그림 2. MA plot은 logCPM 대비 logFC, volcano는 logFC 대비 −log10 P — 유의(FDR<0.05) 유전자를 색으로 구분.

• MA plot — 발현 변화가 평균 발현에 따라 어떻게 분포하는지. 저발현 구간에서 점이 양옆으로 퍼지면 정상입니다.
• volcano plot — 위로 갈수록 유의, 좌우로 갈수록 효과 큼. 오른쪽 위·왼쪽 위 모서리가 신뢰할 만한 DEG입니다. (volcano 해석은 Volcano plot 글 참고)
• 산포 진단도 한 줄로 — plotBCV(y)는 발현 수준별 BCV를, plotQLDisp(fit)은 QL 산포를 보여 줍니다.

💡 FDR/padj 컷오프 — FDR < 0.05는 "유의로 부른 유전자의 평균 5%가 위양성"이라는 뜻입니다. logFC 컷오프(예: |logFC| > 1, 즉 2배)를 같이 걸면 통계적 유의 + 생물학적 크기를 동시에 봅니다. 자세한 원리는 p값과 FDR 글에.

11. edgeR vs DESeq2 — 정규화·모델·검정·언제 무엇을

같은 데이터를 넣으면 보통 상위 DEG 목록이 80~90% 겹칩니다(둘 다 NB + 경험적 베이즈라서). 차이는 세부 설계에 있습니다.

항목	edgeR	DESeq2
정규화	TMM (trimmed mean of M-values)	median-of-ratios (RLE)
분포 모델	음이항(NB) GLM	음이항(NB) GLM
산포 추정	경험적 베이즈, common/trended/tagwise	경험적 베이즈, 추세 적합 후 수축
logFC 처리	원시 logFC (필요 시 glmTreat)	logFC도 수축(lfcShrink)
기본 검정	QL F-검정(glmQLFTest)	Wald 검정(LRT 옵션)
저카운트·소표본	v4에서 강화(연속 NB·편향보정)	독립 필터링·이상치 보정 내장
입력 객체	DGEList	DESeqDataSet
속도(대규모)	v4에서 빠름(estimateDisp 생략 가능)	중간

언제 무엇을?

• 둘 중 아무거나 OK — 표준 2그룹·반복 있는 설계라면 결과가 거의 같습니다. 팀·논문 관례를 따르세요.
• edgeR이 유리 — 샘플이 매우 적거나 카운트가 낮을 때(QL F-검정의 보수성), 복잡한 대비(contrast)를 세밀히 짤 때, 대규모 데이터 속도.
• DESeq2가 유리 — logFC 자체를 순위·시각화에 쓸 때(lfcShrink로 저카운트 과장 억제), 이상치 자동 처리가 편할 때.

🧭 두 도구 결과가 크게 다르면 그건 도구 탓이 아니라 필터링·설계·정규화 가정을 다시 볼 신호입니다.

그림 3. 정규화(TMM vs median-of-ratios)·모델(NB)·검정(QL F vs Wald)·언제 쓰나를 한 장으로.

12. 자주 나는 에러 & 해결

• Error in glmQLFit ... design — model.matrix(~ group)을 빠뜨렸거나 그룹과 행렬의 샘플 순서가 어긋난 경우. colnames(cts)와 group 순서를 맞추세요.
• 유의 유전자가 0개 — 필터링이 과했거나(filterByExpr 전에 임의로 더 잘라냄), estimateDisp/calcNormFactors를 건너뛴 경우. 순서대로 다시 실행하세요.
• could not find function "glmQLFit" — library(edgeR)를 안 했거나 구버전. packageVersion("edgeR")로 4.x인지 확인.
• exactTest에 design을 넣으려 함 — exactTest는 설계 행렬을 받지 않습니다(2그룹 전용). 공변량이 있으면 glmQLFit 경로로 가세요.
• logFC 부호가 반대 — factor의 levels 순서 문제. 기준(reference)을 levels 맨 앞에 두세요.

13. 확장 — limma-voom (제3의 선택지)

edgeR·DESeq2 외에 limma-voom이라는 세 번째 표준이 있습니다. 카운트를 그대로 NB로 보는 대신, voom()이 분산-평균 관계를 추정해 관측 가중치(precision weight)로 바꾸고, 그 뒤 limma의 선형모형 + 경험적 베이즈(원래 마이크로어레이용)를 적용합니다.

# 같은 y(DGEList, TMM까지 끝난 상태)에서 바로 이어진다
v   <- voom(y, design, plot = TRUE)     # 카운트 → 로그CPM + 가중치
fit <- lmFit(v, design)                 # 가중 선형모형
fit <- eBayes(fit)                      # 경험적 베이즈 조절
topTable(fit, coef = 2)                 # 결과 표(logFC, t, P.Value, adj.P.Val)

• 장점 — 샘플이 많을 때 매우 빠르고 안정적이며, limma의 풍부한 대비·블로킹(duplicateCorrelation) 기능을 그대로 씁니다.
• 선택 기준 — 대략 소표본·저카운트 → edgeR/DESeq2, 표본 많고 복잡한 설계 → limma-voom이 무난합니다. 셋 다 결과는 대체로 수렴합니다.

정리 & 다음 학습

• ✅ edgeR 핵심 경로 — DGEList → filterByExpr → calcNormFactors(TMM) → estimateDisp → glmQLFit/glmQLFTest → topTags.
• ✅ 정규화는 TMM (카운트를 안 바꾸고 norm.factors만 만든다), 검정은 QL F-검정(소표본 위양성 통제).
• ✅ 결과 표는 logFC·logCPM·PValue·FDR — FDR < 0.05 + |logFC| > 1로 DEG 선정.
• ✅ edgeR vs DESeq2 — TMM vs median-of-ratios, QL F vs Wald. 표준 설계에선 결과 80~90% 일치.
• ✅ 세 번째 선택지 limma-voom — 표본 많은 복잡한 설계에 강함.

🎓 다음 학습
• (비교) DESeq2로 차등발현(DEG) 분석 — 같은 분석을 DESeq2로
• (다음 단계) clusterProfiler 농축분석 — DEG 목록 → GO/KEGG 경로 해석
• (참고) Volcano plot 그리기 · RNA-seq 정규화 · PCA로 샘플 QC

References

1. Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1), 139–140.
2. Chen, Y., Chen, L., Lun, A. T. L., Baldoni, P., Smyth, G. K. (2024/2025). edgeR v4: powerful differential analysis of sequencing data with expanded functionality and improved support for small counts and larger datasets. (glmQLFit 산포 직접 추정·연속 NB·normLibSizes)
3. edgeR User's Guide (Bioconductor, edgeR 4.x). DGEList·filterByExpr·calcNormFactors(TMM)·estimateDisp·glmQLFit/glmQLFTest·exactTest·topTags·plotMD/plotBCV.
4. Robinson, M. D., Oshlack, A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data (TMM). Genome Biology, 11, R25.
5. Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 (median-of-ratios·Wald). Genome Biology, 15, 550.
6. 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.
7. 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