第3章：负二项模型基础——edgeR中count数据的分布假设、均值与方差关系

3.1 为什么是负二项？

做RNA-seq差异表达分析，第一步就得琢磨清楚：count数据到底服从什么分布？

我刚开始接触这个领域时，看到别人直接用t检验比较两组count均值，心里就觉得不对劲。count数据是非负整数，而且有明显的“均值小、方差大”的现象——说白了，就是离散度远远超过了泊松分布能解释的范围。

举个实际例子。我曾在某个水稻转录组项目里，看到一个基因在对照组5个样本的count值分别是：3、4、2、6、5。均值是4，方差大约是2.5。泊松分布要求方差等于均值，这里差了快一倍，还没到极端情况。

换成另一个基因，对照组：2、8、3、15、4。均值6.4，方差呢？算出来大概30.8——是均值的近5倍！这种“过离散”在真实数据里比比皆是。

所以edgeR选择负二项分布作为基础模型，我非常认同。它引入了一个额外的离散参数，专门用来捕捉这种方差膨胀。

3.2 数学表达：从泊松到负二项

负二项分布有两种常见参数化方式。edgeR用的是“均值-离散度”参数化，我最喜欢这种，因为参数意义直观。

记第i个基因在第j个样本中的count值为Y_ij。负二项分布假设：

Y_{ij} ~ NB(μ_{ij}, φ_i)
其中：
  μ_{ij} = 样本j的文库大小因子 × 基因i的真实表达水平
  φ_i = 基因i的离散度参数（dispersion parameter）
  方差 Var(Y_{ij}) = μ_{ij} + φ_i × μ_{ij}^2

注意这个方差表达式，它包含了两个部分：

• μ：泊松部分（抽样误差）
• φ × μ²：过离散部分（生物学变异）

我个人习惯把φ理解成“每单位均值平方带来的额外方差比例”。φ越大，过离散越严重。

3.3 均值-方差关系：edgeR的灵魂曲线

方差函数 V(μ) = μ + φμ² 是一条抛物线。我过去带学生时，经常让他们画这个曲线：横轴是均值μ，纵轴是方差V。泊松分布就是一条45度线（V=μ），而负二项模型的曲线会随着μ增大，快速甩开泊松线。

下面这张图展示了不同φ取值下的均值-方差关系：

# 在R里你可以这么画：
mu <- seq(0, 50, length.out = 100)
phi <- c(0.01, 0.05, 0.1, 0.2)
plot(mu, mu, type = "l", lwd = 2, ylim = c(0, 100),
     xlab = "均值 μ", ylab = "方差 V(μ)")
for (p in phi) {
  lines(mu, mu + p * mu^2, col = rainbow(4)[which(phi == p)], lwd = 2)
}
legend("topleft", legend = paste("φ =", phi), 
       col = rainbow(4), lwd = 2)

为什么要强调这个关系？因为在edgeR的quasi-likelihood框架里，方差函数直接决定了准似然方程的权重。你想想看，高通量数据里，弱表达基因（μ很小）和强表达基因（μ很大）的方差可以差好几个数量级，如果不处理好这个关系，后续的似然比检验根本没谱。

3.4 离散度参数：φ到底在描述什么？

φ是负二项模型新增的“自由度”。它涵盖了两类变异：

• 技术重复之间的变异：同一样本在相同条件下测序，count波动——这部分通常很小，因为测序技术已经比较成熟
• 生物学个体之间的变异：同一处理组不同个体间的真实表达差异——这才是大头

edgeR在估计φ时，默认假设所有基因共享一个公共离散度（common dispersion），然后根据数据调整。后来有了tagwise dispersion（gene-wise dispersion），再后来就是quasi-likelihood的稳健估计。

离散度类型	edgeR函数	适用场景
公共离散度	estimateCommonDisp()	少量样本、基因数少的情况
测序深度趋势离散度	estimateTrendedDisp()	中等样本量，表达量跨度大
基因特异离散度	estimateTagwiseDisp()	样本量充足（≥3个重复每条件）
准似然离散度	estimateQLDisp()	复杂实验设计，需要稳健性

我个人建议：除非你的样本量极少（比如每组只有2个重复），否则不要只用common dispersion。它太粗糙，会低估高表达基因的方差，同时高估低表达基因的方差。换成tagwise或QL都会好很多。

3.5 本章知识结构图

下面这张SVG图总结了负二项模型在edgeR中的完整逻辑链：

3.6 代码实战：用edgeR感受负二项模型

理论说了半天，不如来段代码直接感受。假设我们有一个简单的DGEList对象：

library(edgeR)
# 模拟一个简化情况：10个基因，3个样本
set.seed(2025)
counts <- matrix(rnbinom(30, mu = 10, size = 5), nrow = 10, ncol = 3)
group <- factor(c("A", "A", "B"))
dge <- DGEList(counts = counts, group = group)

# 查看负二项假设下隐式假定的均值-方差关系
# 计算每个基因的均值与方差
mu_hat <- rowMeans(counts)
var_hat <- apply(counts, 1, var)

# 用edgeR的common dispersion估计
dge <- estimateCommonDisp(dge)
phi_common <- dge$common.dispersion
cat("公共离散度 phi =", phi_common, "\n")

# 预测方差（负二项模型）
var_nb <- mu_hat + phi_common * mu_hat^2

# 对比样本方差与模型预测方差
print(cbind(均值 = mu_hat, 样本方差 = var_hat, 负二项预测方差 = var_nb))

运行后你会发现，样本方差和负二项预测方差大致吻合，但不会完全相等。这是因为我们只有3个样本，每个基因的方差估计本身就有噪音。这也是为什么edgeR后来引入quasi-likelihood的原因之一——它不要求精确的分布假设，只要求方差结构写对了就行。