1、DESeq2 简介与背景
RNA-seq 差异表达分析,说白了就是回答一个问题:「基因 A 在对照组和处理组之间,到底有没有显著变化?」。我在刚入行时,觉得这事很简单——算个倍数、做个 t 检验不就完了?结果第一次跑真实数据,假阳性率高得吓人。后来才明白,RNA-seq 数据的「脾气」很特殊,普通统计方法根本镇不住场子。DESeq2 就是专门为这种数据设计的工具,它背后的统计学原理,是这个领域绕不开的核心。
RNA-seq 差异表达分析概述
差异表达分析的目标,是从成千上万个基因中找出那些在不同条件下表达水平有显著变化的基因。听起来像「大海捞针」,但难点不在于数量多,而在于数据的噪声结构。
RNA-seq 数据有几个明显特点:
- 计数数据——每个基因的 read 计数是非负整数,不是连续值。
- 动态范围大——高表达基因可能有几万条 read,低表达基因只有几条甚至零条。
- 过度离散——生物学重复之间的变异,往往比泊松分布预期的要大。
- 样本量小——通常每个条件只有 3~6 个生物学重复,经典统计方法很难施展。
我见过不少初学者直接拿标准化后的 FPKM 去做 t 检验,结果跑出来一堆「显著基因」,换一批样本全不重复。为什么?因为 FPKM 本身的统计分布没有被正确建模。DESeq2 从原始计数出发,用负二项分布直接建模,绕开了这些中间步骤的偏差。
| 方法 | 输入数据 | 分布假设 | 适用场景 |
|---|---|---|---|
| t 检验 / limma | log2(FPKM/RPKM) | 正态分布 | 样本量大、数据已标准化 |
| edgeR | 原始计数 | 负二项分布 | 小样本、过度离散数据 |
| DESeq2 | 原始计数 | 负二项分布 + 自适应离散度 | 小样本、离散度不均、有离群值 |
| voom + limma | 原始计数 → log2(CPM) | 均值-方差关系建模 | 中等样本量、设计复杂 |
为什么是负二项分布?
RNA-seq 计数数据最直觉的分布假设是泊松分布。但泊松分布有一个硬约束——均值等于方差。我曾在项目中检验过这个假设,拿一组真实数据的基因表达量画了均值-方差图,结果方差往往是均值的 2~8 倍。这就是典型的过度离散问题。
负二项分布(也叫 Gamma-Poisson 分布)在泊松基础上引入了一个额外的离散度参数,让方差可以独立于均值变化:
方差 = μ + α × μ²
其中 μ 是均值,α 是离散度参数。当 α = 0 时,负二项退化为泊松分布。α 越大,数据越「散」。
DESeq2 的巧妙之处在于,它不假设所有基因共享同一个 α,而是让 α 随基因的均值自适应调整——低表达基因的离散度被「压缩」向整体趋势,高表达基因则保留自身估计。这个策略我特别喜欢,因为它兼顾了统计强度和生物学异质性。
核心理解: DESeq2 用负二项分布建模计数数据,通过「适度压缩」的离散度估计,在小样本下依然能稳健地识别差异表达基因。这是它相比早期工具最大的优势。
Wald 检验:差异基因的统计裁判
模型建好了,接下来要回答「某个基因的变化是否显著」。DESeq2 提供了两种检验方式,最常用的是 Wald 检验。
Wald 检验的思路非常直接:用系数估计值除以它的标准误,得到一个 z 统计量,然后对照标准正态分布算出 p 值。在 DESeq2 的语境里,这个「系数」就是处理组相对于对照组的 log2 倍数变化。
Wald 统计量 = β / SE(β)
β: log2 倍数变化的估计值
SE(β): 该估计值的标准误
如果 β 远大于它的标准误,说明变化信号强过噪声,p 值就小。嗯,这里要注意——Wald 检验在小样本或稀疏数据下可能会偏激进。我曾经在某个低表达基因上看到 Wald 检验给出 10⁻⁶ 级别的 p 值,但可视化一看,两组数据几乎重叠。后来检查发现,这个基因的离散度被压缩得太低,导致标准误被低估了。
避坑指南: 我曾经在筛选标记基因时,过度依赖 Wald 检验的 p 值,结果选出的基因在独立验证集中表现很差。后来我养成了一个习惯:始终结合 log2FC 和 baseMean 做双重筛选,不要只看 p 值。DESeq2 的 results() 函数默认会做独立过滤,但你自己仍要设置合理的 log2FC 阈值。
除了 Wald 检验,DESeq2 还支持似然比检验(LRT),适用于比较嵌套模型——比如你想检验「时间 + 处理」与「只含时间」的模型是否有显著差异。LRT 更保守,在探索性分析中我更喜欢用它来减少假阳性。
DESeq2 分析流程总览
把上面的原理串起来,完整的 DESeq2 分析流程就像一条流水线:从原始计数出发,依次经过归一化、离散度估计、模型拟合、假设检验,最后得到差异基因列表。下面这张图是我自己梳理的流程框架,几乎每个项目我都会先对照它过一遍。
个人习惯: 每次拿到新的 RNA-seq 数据,我都会先用这个流程框架过一遍,确认每一步的输入输出是否对得上。尤其在离散度估计这一步,我会用 plotDispEsts() 画图检查压缩效果——如果离散度散点图出现明显的「阶梯状」偏离,说明数据可能有批次效应或样本污染,这时候跑出来的差异基因我是不敢直接用的。
最后聊一下实际使用时最简单的代码入口。DESeq2 的核心函数调用非常简洁,但背后跑的就是上面整个流程:
library(DESeq2)
# 构建数据对象
dds <- DESeqDataSetFromMatrix(
countData = counts_matrix,
colData = sample_info,
design = ~ condition
)
# 一键执行差异表达分析
dds <- DESeq(dds)
# 提取结果
res <- results(dds, alpha = 0.05)
你想想看,这四行代码背后做了多少事:归一化、离散度估计、负二项 GLM 拟合、Wald 检验、独立过滤、多重假设校正…… 理解每个环节的原理,才能用好这个工具。我见过太多人把 DESeq() 当黑盒跑,结果换了数据不会调参,出了问题不知道怎么排查。嗯,这正是为什么我们要花一整个章节来聊背景和原理。