4. 模型假设：负二项分布与均值-方差关系
说到DESeq2和edgeR的结果为什么经常不一样，我得先聊一个核心问题——它们俩对数据的基本假设。说实话，很多人跑完差异分析，看到两个工具给的结果不同，就开始怀疑代码写错了。其实啊，很可能是忽略了这两个工具的模型起点就不完全一样。
我习惯把这个问题归结为三个字：方差假设。你想想看，RNA-seq数据的本质是什么？是计数。对计数数据建模，最自然的想法是泊松分布。但为什么实际没人直接用泊松呢？
核心结论：RNA-seq数据的方差通常大于均值，违反了泊松分布“方差=均值”的假设。DESeq2和edgeR都采用负二项分布来建模，引入了一个额外的离散度参数来处理过度离散。
为什么泊松分布不够用？
先来个简单的例子。假设我们测了一个基因，它在两个条件下的平均表达量都是100。按照泊松分布，方差也应该是100，标准差就是10。也就是说，大部分样本的计数会在90到110之间波动。
但实际RNA-seq数据呢？我做过一个项目，某个看家基因的平均表达量是5000，可6个重复样本的计数却是：4200、5100、4900、5600、4800、5300。你算算方差，远大于均值。这种方差过大的现象，统计上叫过度离散。
| 分布类型 | 均值-方差关系 | 适用场景 |
|---|---|---|
| 泊松分布 | 方差 = 均值 | 技术重复、测序深度极高且生物学变异极小 |
| 负二项分布 | 方差 = 均值 + 离散度 × 均值² | 生物学重复、真实RNA-seq数据 |
表格里这个公式很关键。你看，负二项分布的方差由两部分组成：一部分是泊松部分的均值（测序噪声），另一部分是一个二次项，由离散度参数控制。这个二次项就是处理过度离散的关键。
均值-方差关系的具体形式
DESeq2和edgeR都用负二项分布，但它们在如何估计和处理离散度上存在差异。这就导致了最终结果的不同。
先看基础公式。假设基因g在样本j中的真实表达量为μgj，那么它的计数Kgj服从负二项分布：
K_gj ~ NB(μ_gj, α_g)
Var(K_gj) = μ_gj + α_g * μ_gj²这里的α_g就是基因特异的离散度参数。当α_g = 0时，负二项就退化成了泊松。α_g越大，方差偏离均值越厉害。
我在项目中遇到过一种情况：某个基因的表达量很低，平均只有5左右，但它的离散度特别大。结果DESeq2和edgeR对这个基因的p值判定就完全相反。为什么？因为两个工具对离散度的先验假设和收缩方式不同。
个人习惯：我每次跑差异分析前，都会先画一个均值-方差散点图。如果看到方差随均值呈二次增长的趋势，那就放心用负二项。如果方差和均值差不多是一条斜线，那可能数据本身就有问题（比如重复样本太相似）。
DESeq2 vs edgeR：离散度估计的差异
来个表格，清晰对比一下两个工具在模型假设上的关键差异：
| 维度 | DESeq2 | edgeR |
|---|---|---|
| 离散度先验 | 拟似然先验，默认所有基因共享一个全局趋势 | 经验贝叶斯，基于加权条件似然 |
| 收缩方式 | 将基因特异离散度向全局趋势收缩 | 将基因特异离散度向全局共同离散度和趋势离散度收缩 |
| 低表达基因处理 | 强收缩，低表达基因离散度被拉向全局趋势 | 中等收缩，保留更多基因特异信息 |
| 均值-方差曲线 | 参数化拟合（由所有基因拟合一条光滑曲线） | 非参数化或半参数化拟合 |
看这个表格，你就能明白为什么结果会有差异了。特别是低表达基因，DESeq2的收缩更狠，结果更保守；edgeR相对更放松，容易筛出更多基因。我自己做模拟数据对比时，这个差异非常明显。
什么时候两个工具的结果会一致？
这里分享一个经验：当生物学重复足够多（比如每组6个以上）时，两个工具的结果高度一致。因为重复多了，离散度估计就稳定了，先验假设的影响就小了。
我曾经踩过的坑：有个项目每组只有2个生物学重复，我同时跑了DESeq2和edgeR。结果筛选出的差异基因，只有30%是重叠的。当时我还以为是数据有问题，后来发现就是离散度估计不稳定导致的。
均值-方差关系的可视化理解
直接看公式和表格可能有点抽象。我来画一张图，把整个逻辑串起来。
这张图展示了从原始计数数据出发，到最终模型假设建立的完整流程：
这张图把整个逻辑链条串起来了。从原始计数出发，我们发现泊松分布太理想化，实际数据方差更大。所以必须引入负二项分布，用离散度参数α来描述均值-方差关系。而DESeq2和edgeR的差异，归根结底就是对α的估计策略不同。
实践中的关键参数：离散度
在实际操作中，有几个数值值得关注：
- α < 0.01：离散度很低，数据接近泊松，DESeq2和edgeR结果几乎一致。
- α 在 0.01 ~ 0.1：中度离散，两个工具开始在低表达基因上出现差异。
- α > 0.1：高离散度，结果分歧显著，需要谨慎解释。
我个人习惯在跑差异分析前，先跑一个简单的离散度估计：
# R代码示例：查看基因的均值-方差关系
library(DESeq2)
# 假设dds是已构建的DESeqDataSet
plotDispEsts(dds, main="Mean-Dispersion Relationship")
# 如果看到点（基因）围绕红色曲线紧密分布，说明模型拟合良好
# 如果点很分散，可能是数据质量有问题避坑指南：我曾经处理过一批肿瘤样本数据，离散度特别大（α接近0.5）。当时两个工具给出的差异基因列表，重叠率连40%都不到。后来我重新检查了样本，发现有3个样本的library size异常大，去除了它们之后，离散度降到了0.05，两个工具的结果也基本一致了。
总结一下核心要点
- RNA-seq计数数据的方差 > 均值，必须用负二项分布处理过度离散。
- 负二项分布的方差公式为：μ + α·μ²，α是离散度参数。
- DESeq2和edgeR在离散度α的估计方式上存在差异，导致结果不同。
- 样本量越少、离散度越大，两个工具的结果差异越明显。
- 检查均值-方差关系，是判断模型是否合适的第一步。
嗯，到这里，模型假设这部分就讲清楚了。你理解了负二项分布和均值-方差关系，就等于抓住了DESeq2和edgeR差异的根源。下一块内容，我们具体看看这两个工具是怎么处理那些异常值的。