DESeq2 简介：设计理念与批次校正之道

DESeq2 这名字，做转录组的朋友肯定不陌生。我入行那会儿，市面上还没有这么顺手的工具。

说白了，它的核心任务就一个：从基因表达计数矩阵里，找出哪些基因在不同条件下真的在差异表达。你想想看，一次 RNA-seq 实验下来，基因那么多，样本还可能来自不同批次、不同测序 lane。怎么把真正的生物学差异和那些乱七八糟的技术噪声分开？DESeq2 就是干这个的。

DESeq2 的设计理念

DESeq2 的设计哲学，我总结为三点：稳健、灵活、不做过度假设。它基于负二项分布模型来处理计数数据——为什么不用普通的泊松分布？因为真实的生物学重复之间，方差往往比均值大，这就是“过度离散”。泊松模型会低估方差，导致假阳性。

嗯，这里要注意，DESeq2 不像某些方法那样，上来就把数据 log 转换。它直接对原始整数计数建模，把测序深度和基因长度的影响通过“size factor”和“normalization factor”吸收掉。我最早用 DESeq 的时候，总觉得这步有点绕，后来才明白——这才是不失真的大前提。

DESeq2 里的 shrinkage 估计也是神来之笔。对于表达量很低或者离散度不稳定的基因，它会向全局趋势“收缩”，避免那些不靠谱的极端值出来捣乱。我在项目中遇到过，有些基因只有几个 counts，如果不 shrinkage，p 值会忽高忽低，校正完一个显著的都找不到。用了 shrinkage，结果一下就稳定了。

从原始计数到差异基因——核心流程

这个流程其实不复杂，我来带你捋一遍：

  • 1. 构建 DESeqDataSet — 把 count 矩阵、样本信息和实验设计 formula 组合成一个专门的数据对象。这一步就像搭好火锅的底料，后面所有操作都基于它。
  • 2. 估计 size factor — 用来校正不同样本间的测序深度差异。我习惯用默认的中位数比值法，当然也可以用自己定义的 factor。
  • 3. 估计基因离散度 — 每个基因单独估计初始离散度，然后拟合一条趋势线，再把每个基因的离散度向趋势线做 shrinkage。
  • 4. 拟合模型并做 Wald 检验或 LRT — 根据你的设计公式，拟合 NB 模型，然后比较 full model 和 reduced model，或者直接做 Wald 检验。
  • 5. 多重假设检验校正 — 用 BH 方法控制 FDR，选出真正的差异基因。

你看，上面这五步，缺一不可。我曾经见过有人跳过离散度 shrinkage 这一步，结果跑出来的差异基因列表里，全是那些低表达基因，完全没法用。所以，该走的流程，一步都不要省。

核心代码段（我习惯的写法）：

dds <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, alpha = 0.05)

批次效应校正：把“实验室差异”写进模型

批次效应，是所有高通量实验绕不开的坎。同一批样本分两批跑测序，你今天跑的和下周一跑的，结果可能差了十万八千里。

DESeq2 校正批次效应，最标准的方法就是把批次作为一个变量，直接放到设计公式里。举个例子，如果你的实验有两个条件（treated vs control），同时有三个测序批次（batch1, batch2, batch3），你的 design 应该写成这样：

design(dds) <- ~ batch + condition

这样一来，DESeq2 在建模时会把批次当作协变量，把批次效应“吸收”掉，只关心 condition 的效应。我特别推荐这种做法——因为它是直接在模型层面做校正，不会干扰原始 counts 的整体结构。

但有个坑，我必须提醒你：

我曾经踩过的坑：

如果某个条件只出现在某一个批次里，比如所有 treated 样本都在 batch1，所有 control 样本都在 batch2。这种情况下，batch 和 condition 完全混杂。你如果还坚持在 design 里同时放 batch 和 condition，模型会告诉你“不能估算”——实际上，这是实验设计的问题，不是 DESeq2 的问题。正确的做法是重新设计实验，或者在分析时放弃批次校正，仅保留 condition。

有时，你拿到了别人已经标准化好的数据，或者想做一个更直观的批次校正可视化。这时可以用 limma::removeBatchEffect() 或 DESeq2 提供的 varianceStabilizingTransformation 加 assay 操作。我个人的经验是，对于 PCA 或热图展示，用 vst 或 rlog 转完再移除批次效应，效果非常干净。

小技巧：

如果你不确定设计公式该怎么写，不妨先跑一个 plotPCA(rlog(dds))，看看样本的聚类情况——同一批次的样本是不是扎堆了？如果是，那就说明批次效应确实存在，必须放进模型里面去。

不同批次设计下的公式写法对照

不同的实验结构，对应的公式写法也不同。我整理了一个对照表，你直接用：

实验场景 设计公式 说明
单因素，无批次 ~ condition 最简单的情形
单因素 + 批次 ~ batch + condition 批次放在前面，作为协变量
多因素 + 批次 ~ batch + factor1 + factor2 + factor1:factor2 批次+主效应+交互效应
批次与条件混杂 ~ condition 此时无法校正，设计缺陷，只能放弃批次变量

说实话，批次效应校正这件事，最关键的其实不是代码怎么写，而是你得先搞清楚——自己的实验到底有没有批次效应？它在什么位置？能不能跟条件变量分开处理？

我个人的工作流是：拿到新数据后，第一件事就是做 PCA，先看看数据概貌，再决定要不要放 batch 变量。如果 PCA 图上，第一主成分完全被批次划分，那就别犹豫了，把 batch 加进 design 里，准没错。如果批次的影响只在第二或第三主成分体现，那也得加，只是影响没那么大而已。

总之，DESeq2 的批次校正方法很成熟，但用的时候一定要留个心眼——设计公式写对了，一切都顺；写错了，结果就会骗人。

DESeq2 核心思想与批次校正流程 原始 counts 估计 size factor 离散度估计 + shrinkage 拟合 NB 模型 设计公式加入 batch Wald 检验 / LRT 校正批次效应后 差异基因列表 要点：批次变量需与条件变量可分离；设计公式中的顺序影响不大，但批次要放在前面作为协变量 避坑：若批次与条件完全混杂，模型报错“无法估算”，此时只能放弃批次校正

这张图展示了我心目中 DESeq2 的完整工作流——从原始 counts 开始，经过 size factor 校正、离散度估计与收缩、模型拟合，再到加入批次变量后的统计检验，最后输出干净的结果。你写代码的时候，就照着这个流程走，基本不会出错。

← 上一章 📖 返回目录 下一章 →