DESeq2中手动计算系数与自动对比的差异原因探究

你在DESeq2中构建了包含replicate作为阻塞因素的模型：

dds <- DESeqDataSetFromMatrix(countData = CPEB4_featureCounts_3utr_matrix,
                              colData = CPEB4_sample_list,
                              design = ~ replicate + sample_name)
dds <- DESeq(dds)

对应的元数据如下：

sample_name replicate
0195_2022       INPUT         4
0196_2022         IgG         4
0197_2022       CPEB4         4
0198_2022       INPUT         5
0199_2022         IgG         5
0200_2022       CPEB4         5
2125_2021       INPUT         1
2126_2021         IgG         1
2127_2021       CPEB4         1
2235_2021       INPUT         2
2237_2021       CPEB4         2
2238_2021       INPUT         3
2239_2021         IgG         3
2240_2021       CPEB4         3

你尝试了两种方式提取CPEB4 - IgG的对比：

方式1：自动对比 #

CPEB4vsIgG <- results(dds, contrast=c("sample_name","CPEB4","IgG"))

得到的差异基因统计结果：

summary(CPEB4vsIgG)
# out of 17300 with nonzero total read count
# adjusted p-value < 0.1
# LFC > 0 (up)       : 598, 3.5%
# LFC < 0 (down)     : 30, 0.17%
# outliers [1]       : 0, 0%
# low counts [2]     : 7637, 44%
# (mean count < 41)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results

方式2：手动计算系数均值 #

mod_mat <- model.matrix(design(dds), colData(dds))
CPEB4 <- colMeans(mod_mat[dds$sample_name == "CPEB4", ])
IgG <- colMeans(mod_mat[dds$sample_name == "IgG", ])
CPEB4vsIgG_2 <- results(dds,  contrast = (CPEB4 - IgG))

得到的差异基因统计结果：

summary(CPEB4vsIgG_2)
# out of 17300 with nonzero total read count
# adjusted p-value < 0.1
# LFC > 0 (up)       : 672, 3.9%
# LFC < 0 (down)     : 81, 0.47%
# outliers [1]       : 0, 0%
# low counts [2]     : 7637, 44%
# (mean count < 41)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results

你检查了两组系数：

CPEB4
#      (Intercept)       replicate2       replicate3       replicate4       replicate5   sample_nameIgG 
#              1.0              0.2              0.2              0.2              0.2              0.0 
# sample_nameINPUT 
#              0.0 

IgG
#      (Intercept)       replicate2       replicate3       replicate4       replicate5   sample_nameIgG 
#             1.00             0.00             0.25             0.25             0.25             1.00 
# sample_nameINPUT 
#             0.00 

你疑惑：为什么两种方法结果不同？且不纳入replicate时结果一致？

这个差异的核心在于两种对比方式的实际含义完全不同，而你的实验设计（replicate2缺少IgG样本）放大了这个差别：

1. 自动对比的实际含义 #

当你使用contrast=c("sample_name","CPEB4","IgG")时，DESeq2会直接构建一个精准的对比向量：

• 这个向量中，所有replicate相关的系数权重都是0，仅保留sample_nameCPEB4（系数+1）和sample_nameIgG（系数-1）的差值。
• 它的实际意义是：在控制了replicate效应的前提下，CPEB4相对于IgG的平均差异——相当于在每个有配对的replicate（1、3、4、5）中计算CPEB4 vs IgG的差异，再合并这些差异，完全排除了replicate本身的效应干扰。

2. 手动均值法的问题 #

你手动计算CPEB4组和IgG组的模型矩阵均值再做差，这里的问题出在样本分布不均衡：

• 你的CPEB4组有5个样本（覆盖了replicate1-5），而IgG组只有4个样本（缺少replicate2的样本）。
• 这导致计算IgG组的模型矩阵均值时，replicate2的系数均值是0（没有IgG样本在replicate2），而CPEB4组的replicate2系数均值是0.2（5个样本中有1个属于replicate2）。
• 最终你的手动对比向量中，除了sample_name的系数差（-1 for IgG, +1 for CPEB4），还混入了replicate2的系数差（0.2 - 0 = 0.2），以及其他replicate的微小系数差（比如replicate3的0.2 - 0.25 = -0.05）。

简单说，你的手动对比实际上是CPEB4 vs IgG + replicate2的效应差异，这和自动对比想要的“纯CPEB4 vs IgG差异”完全不是一回事，自然会得到不同的差异基因结果。

3. 为什么去掉replicate后结果一致？ #

当模型不纳入replicate时，模型矩阵只有(Intercept)、sample_nameIgG、sample_nameCPEB4三个系数：

• CPEB4组的均值向量是：(Intercept)=1, sample_nameIgG=0, sample_nameCPEB4=1
• IgG组的均值向量是：(Intercept)=1, sample_nameIgG=1, sample_nameCPEB4=0
• 手动计算的差值向量是(0, -1, 1)，和自动对比的向量完全一致，所以结果自然相同。

如果你想得到“控制replicate后CPEB4 vs IgG”的对比，一定要用DESeq2的自动对比方式，手动均值法在样本分布不均衡时会引入额外的无关效应，导致结果偏差。

内容的提问来源于stack exchange，提问作者Ilario De Toma