3. DESeq2基本原理与特点：负二项分布模型、数据归一化、统计检验方法
好，咱们进入第三章节。讲DESeq2，我习惯先给它定个性：这工具专为小样本、高噪声的RNA-seq数据设计的。说白了，它解决的问题很实际——你手头可能就两三个生物学重复，转录组数据又特别零散，怎么才能靠谱地找出差异基因？DESeq2的思路，就是把数据的「离散性」和「均值-方差关系」都考虑进去。
3.1 负二项分布模型：为什么不是泊松？
最早一批差异表达工具用的是泊松分布。但有个问题：泊松假设均值和方差相等。真实RNA-seq数据中，方差几乎总是大于均值——这叫「过离散」。我在项目中处理过一个酵母数据集，技术重复间方差还能接受，一旦换成生物学重复，方差直接翻了几倍。泊松模型在这时就完全不适用了。
DESeq2采用负二项分布，本质上是在泊松基础上加了一个伽马分布的误差项，用来捕捉生物学变异。它的概率质量函数长这样：
P(K = k) = Γ(k + θ) / (Γ(θ) * k!) * (μ / (μ + θ))^k * (θ / (μ + θ))^θ其中：
- μ —— 均值，代表基因表达水平的期望
- θ —— 离散参数，控制方差的大小，θ越小方差越大
- k —— 观测到的计数值
为什幺要这样设计？因为对于同一个基因，不同生物学重复之间的表达量波动，服从的就是这种「厚尾」分布——大的离群值概率比正态分布更高。我记得有次分析免疫处理后的T细胞数据，有一个基因在某一重复中出现极端高表达，用泊松模型直接给出极显著的p值，换DESeq2的负二项模型一算，嘿，没那么夸张。因为模型把离散性算进去了，就不会被个别点带偏。
3.2 数据归一化：DESeq2的独特策略
归一化是所有差异分析的前提。你想想看，两个样本测序深度不同、文库大小不同，直接比原始计数毫无意义。但怎么归一化？不同工具有不同做法。
DESeq2用的是 中位数比值法（Median-of-Ratios）。这个思路我特别喜欢，因为特别直观：
- 对每个基因，计算它在样本i中的计数与所有样本中几何均数的比值
- 取该样本所有基因比值的中位数作为归一化因子
- 用这个因子缩放原始计数
数学上表示为：
s_j = median_{g ∈ G} ( k_{gj} / ( ∏_{j'=1}^m k_{gj'} )^{1/m} )其中 s_j 是样本j的归一化因子，k_{gj} 是基因g在样本j中的计数。
举个例子方便理解：假设你有三个样本，基因A的计数分别是10、100、1000。几何均数是 (10×100×1000)^(1/3) ≈ 100。那么比值就是0.1、1、10。取中位数……如果靠比值集中在1附近，说明这个样本不需要大幅缩放；如果多数基因的比值都小于1，说明这个样本文库偏小，需要向上调整。
归一化前后效果，我用一个表格总结下实际项目中常见的差异：
| 指标 | 原始计数 | DESeq2归一化后 |
|---|---|---|
| 样本间总计数差异 | 可达3-5倍 | 基本消除 |
| 低表达基因波动 | 噪声主导 | 方差稳定 |
| 高表达基因可比性 | 受文库大小影响大 | 直接可比 |
| 离群样本的影响 | 严重偏离 | 被中位数压缩 |
3.3 统计检验方法：Wald检验与LRT
归一化完事儿，数据准备好了，接下来就是检验差异。DESeq2提供了两种检验方式，我分别说说。
3.3.1 Wald检验
这是DESeq2默认的检验方式。思路很简单：先拟合负二项模型，得到每个基因在处理组和对照组间的对数倍数变化（log2 fold change）估计值及其标准误，然后用Wald统计量检验这个倍数变化是否显著不等于0。
公式不复杂：
W = β / SE(β)其中β就是log2倍变化估计值，SE是标准误。W近似服从标准正态分布，p值也就有了。
我个人做快速筛选时更喜欢Wald检验——它运行快，结果直接。一般样品量在3-5个重复时，Wald检验的表现很不错。
3.3.2 似然比检验
LRT适合更复杂的实验设计。比如你不仅有处理vs对照，还有时间点、剂量梯度等因素。LRT的思想是：比较一个全模型（包含你感兴趣的变量）和一个简化模型（去掉该变量），看加入该变量后模型拟合是否有显著改善。
统计量是：
D = -2 * (logLik_reduced - logLik_full)D服从卡方分布，自由度等于两个模型参数个数之差。
我举个例子：一个2×2因子设计（基因型×处理），你想看看有没有交互作用。全模型包含基因型、处理以及基因型×处理的交互项；简化模型只包含主效应。用LRT比较两者，就能判断交互作用是否显著。
3.4 核心流程总览
下面这张图，把DESeq2的整个分析流程串起来了。你看一遍就能把握它的骨架结构：
嗯，这张图基本就把DESeq2的骨架画出来了。从计数矩阵进来，经过归一化，再到离散参数估计，然后是模型拟合，最后根据实验设计选择合适的检验方法输出结果。每一步之间都有数据传递，缺一步都不行。
我刚开始用DESeq2那会儿，老觉得归一化这步无所谓，直接把原始计数丢进去跑。结果发现同一批数据用三个不同归一化方法跑出来的差异基因列表重叠不到60%。从那以后我才老老实实先把归一化的逻辑搞清楚。你想想看，工具帮你做了归一化，但你至少要知道它怎么做的，否则出问题了连调试方向都没有。
sizeFactors）——如果某些样品的因子特别大或特别小（比如超过3倍范围），往往提示数据质量有问题，或者存在批次效应。我通常习惯先跑一行代码看一下分布：
dds <- estimateSizeFactors(dds)
colData(dds)$sizeFactor总体来说，DESeq2的设计思路非常「生物学家友好」——它不要求你手动调一堆参数，核心的负二项模型已经内化了RNA-seq数据的多重噪声来源。无论是新手还是老手，都可以快速上手并得到可靠的结果。不过话说回来，了解它底层的数学原理，能让你在遇到诡异结果时，保持冷静，找到问题关键所在。