3. DESeq2 介绍：负二项分布模型、设计公式与优势

好，我们来聊聊 DESeq2。这个工具在 RNA-seq 差异表达分析里，可以说是绕不开的一站。它的核心思路，说白了就是给每个基因建立一个统计模型，然后用这个模型去判断表达量的变化到底靠不靠谱。

我记得最早接触 DESeq2 是在一次处理小鼠发育时序数据的时候。样本少，批次乱，用别的工具跑出来一堆假阳性，换 DESeq2 之后结果就稳多了。从那以后，我个人习惯在数据分析初期就先用它快速摸底。

3.1 负二项分布模型：为什么不是泊松？

很多人刚开始会问：测序数据是计数数据，为什么不用泊松分布？

道理很简单。泊松分布假设均值等于方差，但 RNA-seq 数据往往出现 过离散 现象——方差远大于均值。不信你随便挑一个基因，看看它在重复样本里的 read counts，波动通常很大。要是硬用泊松，那显著差异的基因能多到让你怀疑人生。

DESeq2 采用 负二项分布（也叫 Gamma-Poisson 分布）来解决这个问题。它引入了一个额外的离散参数，让模型能同时描述均值 μ 和方差 σ² 之间的关系：

σ² = μ + α * μ²

这里的 α 就是离散参数。当 α = 0 时退化为泊松分布；α 越大，过离散越严重。

关键理解：负二项分布给每个基因配了两个参数——一个控制平均水平（均值），一个控制波动幅度（离散度）。这样，表达量高的基因允许更大的方差，表达量低的基因方差被压制，避免了低表达基因因为一次随机波动就被标成差异。

遇到过一个坑：有个合作的实验室拿到的数据文库质量很差，结果一批低表达基因的离散度被严重高估。我当时建议他们用 estimationType="parametric" 配合稳健的离群基因过滤，才把结果拉回来。

3.2 设计公式：把实验结构告诉模型

DESeq2 的设计公式（design formula）是它的灵魂。你用它告诉模型：我的实验里有哪些因素在影响表达量？ 比如最简单的二分组比较：

design = ~ condition

这里的 condition 就是你实验条件的分组变量（比如 treatment vs control）。

但真实实验往往更复杂——有批次效应、有配对设计、有连续协变量。拿我做过的一个时间序列实验来说，每组有 3 个生物学重复，但样本是分两批测的：

design = ~ batch + condition

这样模型会把 batch 的影响剥离掉，只看 condition 带来的差异。

你还可以加入交互项，比如想看在两种细胞系中药物处理的效应是否不同：

design = ~ cell_line + treatment + cell_line:treatment

我个人建议，在设计公式时遵循一个原则：把所有已知的、可能引入系统偏差的变量都放进去，但别过度拟合。说白了，样本量少的时候加太多参数，反而把信号给吃了。

小技巧：如果你的实验有配对关系（比如同一个病人的癌 vs 癌旁），可以把个体 ID 作为因子放进设计公式：design = ~ patient + condition。DESeq2 会自动把它当作配对模型来处理。

3.3 优势：稳健、灵敏、自带「刹车」

相比其他工具，DESeq2 有几点让人用得很舒服：

  • 稳健的离散度估计：它利用所有基因的信息来「借用强度」（shrinkage），把单个基因的离散度向全局趋势拉。低表达基因不会被单一样本的异常值带偏。
  • 内置的离群值处理：当某个样本在某基因上的 counts 异常高或低时，DESeq2 会用 Cook's distance 检测并自动替换或标记。我曾遇到一个样本在某个基因上有 500 counts 而其他重复都是 10 左右——DESeq2 自动把它标记为离群值，没让这个点影响结果。
  • 自带基因筛选过滤器：默认会丢弃那些平均 counts 太低的基因，减少了多重检验校正的负担。你不需要手动设 CPM 阈值。
  • 支持多种检验：除了两两比较，还支持 Wald 检验和 LRT（似然比检验）。后者适用于多水平或时间序列的全局检验。
特性 DESeq2 edgeR
分布假设 负二项（离散度由所有基因共享信息估算） 负二项（经验贝叶斯估计离散度）
设计公式 支持复杂因子、交互项、连续协变量 支持类似公式（通过 glm 框架）
离群值处理 内置 Cook's distance 自动管理 需要手动标记或使用稳健方法
低表达基因过滤 自动基于平均标准化 counts 手动（如 CPM + 组内非零样本数）
推荐样本量 每组 3 个以上生物学重复 每组 2 个以上即可

注意：样本量太少（如每组 2 个重复）时，离散度估计会非常不稳定。虽然 DESeq2 提供了 fitType="mean" 或 fitType="parametric" 的选项，但我个人建议至少每组 3 个重复，否则我宁愿用 edgeR 的 exact test 配合偏移参数。

3.4 核心工作流程

下面这张图是我自己总结的 DESeq2 分析流程，从输入到结果输出，核心就 5 步：

DESeq2 分析核心流程 1. 构建 count 矩阵 + 设计公式 2. 估计大小因子（标准化） 自动担架、不必手动做 TPM 3. 估计离散度（收缩） 所有基因共享信息 4. 拟合模型 + 假设检验 Wald 检验（两两比较） LRT 检验（多组比较） 5. 提取结果表 + 校正 p 值

整个流程你写成 R 代码其实很简洁。我最常用的最小工作流是这样：

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = meta,
                              design = ~ batch + condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treated", "control"))
resOrdered <- res[order(res$padj), ]

第三行那个 DESeq() 函数实际上依次做了标准化、离散度估计和模型拟合。嗯，一步搞定，但内部执行顺序就是图上那 1→5 步。

3.5 我常用的参数调优

说几个在实际项目里踩过的点：

  • minReplicatesForReplace：默认是 7。如果你每组只有 3-4 个重复，建议改成 minReplicatesForReplace=Inf，关闭自动替换离群值功能，否则可能替换掉本属正常的生物变异。
  • betaPrior：默认 TRUE 会给效应项加先验，使系数收缩，减少假阳性。但在处理多因子交互设计时，我遇到过冲突——有时把交互项也收缩掉了。如果交互项是你关注的重点，可以设置 betaPrior=FALSE。
  • parallel：大数据集（比如几万个基因、几十个样本）时，用 DESeq(dds, parallel=TRUE) 可以并行跑离散度估计，速度提升明显。我最多跑过 60 个样本、6 万基因，并行后从 20 分钟缩到 5 分钟。

经验之谈：别急着调参数。先用默认参数跑一遍，看结果里的离散度图（plotDispEsts(dds)）。如果拟合曲线偏离散点过远，再考虑换 fitType。我以前一上来就改 fitType="mean"，结果低估了高表达基因的方差，导致几十个假阳性……后来老老实实看图再决策。

3.6 什么时候该用 DESeq2

总结一下我个人的选择逻辑：

  • 每组 ≥3 个重复，且实验设计涉及多个因子或协变量 → 首选 DESeq2
  • 样本量很少（每组 2 个）或需要 exact test → 选 edgeR
  • 数据中有明显离群样本 → DESeq2 的自动处理更省心
  • 需要分析多组时间序列或剂量响应 → DESeq2 的 LRT 好用

当然，这两者现在都支持 count 数据格式，设计公式也很类似。我有时候会两个工具都跑一遍，看看结果的一致性。如果某个基因在两个工具里都是显著且方向一致，我就更有信心拿去验证。

最后提一句：负二项模型不是万能的。如果你的数据里零膨胀（大量基因在某些样本里 count=0）非常严重，DESeq2 默认会直接把那些基因过滤掉。遇到这种情况，可以考虑改用 edgeR 的稳健方法或干脆用专门处理 zero-inflated 模型的工具（比如 zinbwave）。不过这是后话了。

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