DESeq2与limma/edgeR的padj结果差异显著是否正常？

我在分析5例健康样本和5例疾病样本的转录组数据时，发现DESeq2、limma（voom流程）和edgeR输出的log2FC结果相近，但padj（校正后P值）差异极大，想确认这种情况是否正常。

DESeq2分析代码及结果 #

dds <- estimateSizeFactors(ddsTxi)
filtered_genes <- rowSums(counts(dds) >= 10) >= (ncol(dds) / 2)
dds_filtered <- dds[filtered_genes, ]
dds_filtered <- DESeq(dds_filtered)
res <- results(dds_filtered, contrast = c("group", "disease", "healthy"), lfcThreshold = 1, cooksCutoff = F)
res_filtered <- res[res$padj <= 0.05, ]
> res_filtered.df
                  baseMean log2FoldChange     lfcSE      stat       pvalue         padj
IGHV3-30         748.81639      12.520299 1.7518336  6.576138 4.828254e-11 
IGLV3-25         787.73323      11.327897 1.4383719  7.180269 6.957447e-13 
IGHV3-49         240.35750      10.881116 1.0428510  9.475099 2.665089e-21 
IGLV2-18         240.34083      10.880756 1.4121829  6.996796 2.618818e-12 
IGKV1-6          471.52213      10.838236 0.9403616 10.462184 1.288512e-25 
IGLV3-19        3394.63949      10.677759 0.7953877 12.167348 4.639110e-34 
IGKV1D-33       1289.91551      10.524039 0.7194581 13.237795 5.308483e-40 
IGLV3-1         1748.26674      10.470238 0.5973499 15.853754 1.324292e-56 
IGLV10-54        180.39686      10.466864 1.6913951  5.597074 2.179993e-08 
IGHV3-53         467.31180      10.466634 0.8807581 10.748279 6.037750e-27 
MMP1           70806.12090      10.453865 0.6319851 14.959001 1.360532e-50 

limma（voom）分析代码及结果 #

group <- factor(targets$group)
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
v.DEGList.filtered.norm <- voom(myDGEList.filtered.norm, design, plot = TRUE)
fit <- lmFit(v.DEGList.filtered.norm, design)
contrast.matrix <- makeContrasts(infection = disease - healthy,levels=design)
fits <- contrasts.fit(fit, contrast.matrix)
ebFit <- eBayes(fits)
myTopHits <- topTable(ebFit, adjust ="BH", coef=1, number=40000, sort.by="logFC")
                   logFC      AveExpr          t      P.Value    adj.P.Val           
MMP1           10.811364  6.656104979  13.344542 1.850951e-08 1.244090e-06  
IGLC3          10.288004  3.492650878   7.518879 7.952738e-06 7.604019e-05  
IGKV1D-33      10.229349  1.116511588   9.331072 8.822073e-07 1.565046e-05  
IGKV2D-28      10.082185  0.898462482   9.523328 7.124072e-07 1.354091e-05  
IGLV3-25       10.010268  0.242702175   6.908813 1.813279e-05 1.431296e-04  
IGLV3-1        10.005929  1.748391611  11.842773 6.883173e-08 2.791843e-06  
IGKV1-5         9.925822  3.363145736  13.850135 1.224740e-08 1.007992e-06 
IGHV1-3         9.886930  0.447926148   6.261069 4.583851e-05 2.964023e-04  
IGLV3-19        9.842028  2.184963159   6.958590 1.692390e-05 1.361195e-04  
IGHV3-21        9.795755  1.417188571   9.024589 1.249408e-06 1.973411e-05  
IGLV2-8         9.773432  2.060931959  13.844952 1.229851e-08 1.007992e-06  
IGKV1D-39       9.681327  2.583496010   9.671441 6.056265e-07 1.202065e-05  
IGHV4-34        9.662238  1.369932576   6.696774 2.441521e-05 1.812382e-04  
IGLV3-21        9.641737  3.087064678   5.233231 2.246810e-04 1.073818e-03  
IGKV1-6         9.638834 -0.112334051   7.940951 4.615981e-06 5.060817e-05  
IGLV2-23        9.624652  3.125414661   9.009834 1.270807e-06 1.999166e-05  

这种padj差异是正常现象，核心原因在于三个工具的统计模型、方差估计策略以及P值校正的基础逻辑存在本质区别：

• 统计模型与方差估计差异
  DESeq2使用负二项分布模型，通过离散度估计（针对每个基因或共享离散度）建模计数数据方差，同时利用收缩估计稳定方差，尤其针对低表达基因，这使得它的原始P值通常更小，后续校正后的padj也更低；limma-voom将计数数据转换为近似正态分布的表达量，再用线性模型结合经验贝叶斯收缩方差，方差估计基于全局趋势，对极端值的处理和DESeq2不同，原始P值相对更大。

• P值校正的基数差异
  DESeq2展示的是padj <=0.05的过滤后结果，而limma输出的是topTable的前N个结果。若DESeq2过滤后保留的基因数远少于limma分析时的总测试基因数，BH校正的基数（总测试次数）不同会导致padj差异——DESeq2的校正基于过滤后的基因集，limma的校正基于所有输入voom的基因，基数越大，padj膨胀越明显。

• lfcThreshold的影响
  DESeq2的results()中设置了lfcThreshold=1，检验的是log2FC绝对值≥1的显著性；而limma默认检验log2FC≠0的显著性。这个阈值设置改变了统计检验的原假设，直接导致原始P值计算逻辑不同，进而影响padj。

总结：三者的padj差异由模型设计、方差估计、检验假设以及校正基数的不同导致，属于正常情况。若要更公平比较，需：

• 统一基因过滤标准，确保三个工具使用完全相同的基因集分析
• 统一检验的log2FC阈值（比如在limma中设置类似阈值，或移除DESeq2的lfcThreshold）
• 确认P值校正的基因总数一致

内容的提问来源于stack exchange，提问作者CrisF