均值-方差关系与dispersion参数深度解析

大家好。今天我们聊一个在RNA-seq差异表达分析里，说不上最难、但绝对是最绕不开的概念——均值-方差关系。我个人习惯把这个关系看作是整个统计建模的地基，地基不稳，上面盖的房子再漂亮也得塌。edgeR的quasi-likelihood方法之所以强大，很大程度上就是因为它把这套关系处理得比别人更细致。

很多做数据分析的同学都问过我：“为什么RNA-seq数据不能用t-test直接比较？”这其实就是在问均值-方差关系为什么重要。说实话，我当年刚接触高通量测序数据时也有过这个疑问。但如果你理解了下面我要讲的技术变异和生物学变异的区别，还有它们如何通过dispersion（离散参数）统一到同一个模型里，那整个edgeR的设计逻辑就一目了然了。

技术变异：泊松分布的局限

先说说技术变异。这是什么？就是你把同一个RNA样本上机测序两次，得到的count数不一样。这是由测序本身的随机采样过程造成的，我们通常用泊松分布来刻画它。

对于泊松分布，一个很关键的性质是：均值等于方差。如果一个基因在三个生物学重复样本中的表达量均值是100，那按照泊松分布预测，这三个样本的计数应该围绕100波动，方差也是100，标准差就是10。听起来挺规整的吧？

我当时做的一个细菌RNA-seq项目，初期用的就是标准泊松模型。结果在差异分析时发现很多基因的重复样本间波动特别大，根本不能用泊松拟合。怎么会这样？

实际上，RNA-seq数据的真实方差几乎总是大于泊松分布的方差。我给大家看一组典型值：想象一个基因，它在样品A、B、C中的count数是120、85、95。均值是100，而样本方差算出来将近300。泊松说你的方差应该是100，但实际是300。这个现象有个专门的名字，叫overdispersion（过离散）。

特征	泊松分布	RNA-seq实际数据
均值与方差关系	均值 = 方差	方差 > 均值
适用范围	技术重复	生物学重复
overdispersion	无	常见且明显

生物学变异：为什么数据会“过离散”

生物学变异才是造成过离散的“真凶”。说白了很简单：你从3只不同的小鼠身上取肝脏组织cuff称重、提RNA、建库、测序。即使这三只小鼠都是同一品系、同龄、同一饮食，它们的基因表达水平本身就不完全一样。这是生物体固有的差异性，不是实验操作失误。

所以RNA-seq数据的总方差，是两部分加在一起的：

• 技术变异：测序采样带来的随机误差
• 生物学变异：个体之间真实的表达差异

这两部分通过一个简单的加法关系融合在一起。在edgeR和limma中，我们用一个叫做负二项分布（negative binomial）的模型来描述这种组合方差。负二项分布在泊松的基础上加了一个额外的参数——dispersion。这个参数就是用来捕捉那个“超出泊松预期”的多余方差。

我个人觉得，用一个公式来展示会更清楚。假设某个基因在样本中的真实期望表达量为μ，那么其方差可以用以下公式表示：

Var(Y) = μ + φ × μ²

这里，φ就是dispersion值。如果φ=0，方差等于均值，模型退化为泊松分布。φ越大，意味着生物学变异相对于技术变异越显著。从我接触过的上百个数据集来看，典型的RNA-seq数据中，φ的值通常在0.01到0.2之间。

dispersion参数：建模的核心

现在重点来了。dispersion参数φ究竟是干什么的？我给大家做一个类比。我们把RNA-seq测序看作是一个测量的过程：技术变异类似测量仪器本身的噪声，而生物学变异像是被测量物体本身的抖动。如果说基因表达量是物体位置，那么dispersion就是“物体抖动幅度”的度量。

在edgeR里，dispersion估计是整个分析流程中最关键的步骤。稍微处理不好，后续所有的p值和FDR就都跟着跑偏了。

edgeR把dispersion分成了几个层级来考虑：

1. Common dispersion：所有基因共用一个离散度。适合样本量少、无法精细建模的情况。
2. Trended dispersion：根据基因的平均表达量拟合出的dispersion趋势线。表达量低的基因，其dispersion往往更大。
3. Tagwise dispersion：每个基因独立估计的dispersion。这种方式最灵活，但也最易受个别基因的表达波动影响。

在实际分析中，我通常推荐使用tagwise dispersion，因为它能更准确地反映每个基因自身的变异模式。但也要注意，如果你的样本量太少（比如每个组只有2个样本），则单个基因的dispersion估计会非常不可靠。这时候，用common dispersion或者趋势线来“收缩”估计，反而更稳定。

我曾经有过一次教训。在做一项人体肿瘤样本的研究时，每组只有3个样本，我毫不犹豫地用了tagwise dispersion。结果分析出来的差异基因列表全被一些低表达的、高变异基因主导了。后来我吸取教训，改用estimateDisp函数默认的趋势收缩策略，效果立竿见影。

从theory到practice：负二项方差公式

现在我想带大家看看负二项分布的方差公式是如何在实际中与我们打交道的。首先，技术变异部分由μ（即泊松的均值）贡献。生物学变异部分由φ × μ²贡献。所以当基因表达水平高时，生物学变异占据主导地位；当基因表达水平低时，技术变异不容忽视。

这个特性直接影响差异表达分析的判定。因为高表达基因的方差大，即使组间差异稍微明显一点，也可能被判定为不显著；而低表达基因的方差小，微小的差异就可能被判定为显著。

为了避免在这种权衡时发生偏差，edgeR利用quasi-likelihood方法对dispersion的不确定性做了进一步量化。这就是我们的课程一直强调的QLF检验的理论背景。

在edgeR中设定dispersion参数

说了这么多理论，我们来看看在edgeR代码里，dispersion到底是如何被触碰的。

首先，你需要通过estimateDisp函数计算dispersion。

y <- DGEList(counts = data, group = group)
y <- calcNormFactors(y)
y <- estimateDisp(y, design)

estimateDisp会同时计算common、trended和tagwise dispersion。如果你希望在差异表达时使用某个特定的dispersion，比如只使用common dispersion，可以直接在后面的glmQLFit中设置：

fit <- glmQLFit(y, design, dispersion = y$common.dispersion)

这么做的缺点是，它完全忽略了基因间的生物学变异差异性。当然，如果样本量太少，这也算是迫不得已的选择。

我个人习惯的做法是：先用plotBCV函数画出均值与BCV（biological coefficient of variation）的关系图。

plotBCV(y)

这个图能帮助我判断数据集的dispersion分布是否合理。green线是trended dispersion，蓝色的点tagwise dispersion。如果你看到tagwise dispersion的点大部分紧贴着trended线，说明我们的收缩策略做得不错。如果点很散乱，甚至出现很多在trended线之上的极端点，就需要仔细检查是不是有样品出了问题。

图形化理解均值-方差关系

现在我画一个框架图，把这个过程可视化一下，方便大家记忆。这张图展示了从测序数据到dispersion估计的整个流程。

这张图能帮你快速理清脉络。从技术变异和生物学变异出发，经过负二项分布建模，最终落脚到三种不同粒度的dispersion估计上。我建议大家把这个图当成后面的学习地图。

小结

均值-方差关系不是什么新概念，但在高通量测序数据分析中却处处是“坑”。理解它，你就能明白为什么差异表达分析不能简单地用t检验；理解它，你就能在分析中选择合适的dispersion建模策略；更重要的，理解它，你就能在汇报结果时，给合作者解释清楚哪些基因的差异是真的，哪些可能是回归假阳性。

如果想进一步深入下去，可以跟踪一下不同的dispersion估计方法对差异基因数量的影响。把这部分搞懂了，后面的quasi-likelihood拟合和QLF检验就好办了。