一个看似奇怪的事实
做RNAseq分析的时候，我们通常会说：表达量计数（read counts）服从负二项分布（Negative Binomial distribution，NB）。这个结论出现在几乎所有主流差异表达分析工具的文档里——DESeq2、edgeR、limma-voom——好像这是一个不言自明的前提。
但为什么？
我查过很多资料，大多数解释要么太浅（“因为数据过度离散，过dispersed”），要么直接跳到公式层面（“NB的概率质量函数是……"），要么扯到什么"等待时间"上去（“NB描述的是等待r次成功所需的试验次数”）——最后一种解释尤其让人困惑，因为RNAseq计数和"等待时间"八竿子打不着。
所以我想把这件事彻底说清楚。用一个最简单的例子，一步步推出为什么NB是自然的选择，而那个经典的等待时间定义，其实和我们的推导毫无关系。
从一个直观的例子说起
假设我们要研究某个基因在细胞A和细胞B中的表达量。
细胞内部发生了什么？ 一个基因要被转录，首先要有RNA聚合酶结合到启动子上，然后开始合成。聚合酶会随机地到达和离开启动子区域，所以转录本身是一个随机过程——在任何给定的时间窗口内，合成出的mRNA分子数是随机的。
如果我们把观察时间固定（比如观察一秒），转录产生的mRNA数量大致服从Poisson分布。理由是：mRNA的产生可以看作一系列独立的"成功"事件（一次转录起始），在一个固定时间窗口内，Poisson是描述这类计数过程的自然模型。
但问题来了——Poisson有一个很强的假设：均值等于方差。
对于单细胞来说，均值和方差的关系确实大致如此。但当我们比较一群细胞（比如组织样本中的上万个细胞）时，情况完全不同。即使是同一个基因，不同细胞的表达量也可能差异巨大——有些细胞可能完全不表达，有些高表达。这种细胞间的异质性（biological variability）会导致整个群体的方差远大于均值。
这就是所谓的"过度离散”（overdispersion）。Poisson处理不了这个问题。
Gamma-Poisson层级模型
解决之道来自一个层级思想。
我们先承认：不同细胞的真实表达率（$\lambda$）本身是不同的。这种差异可以用一个分布来描述。选什么分布？在统计中，Gamma分布是一个自然的选择——因为它数学上处理起来方便，而且足够灵活，可以描述各种形状的随机变量。
于是我们构建一个两层模型：
第一层（细胞间变异）： $\lambda \sim \text{Gamma}(\alpha, \beta)$
这里 $\lambda$ 代表每个细胞的真实表达率（单位时间内的mRNA平均产量）。$\alpha$ 和 $\beta$ 是Gamma分布的参数。均值 $\mu = \alpha/\beta$，方差 $= \alpha/\beta^2$。
第二层（细胞内随机性）： 给定 $\lambda$，计数 $X$ 服从Poisson分布：$X | \lambda \sim \text{Poisson}(\lambda)$
这一层描述的是：即使两个细胞有完全相同的 $\lambda$，由于转录过程本身的随机性，实际观察到的计数也会有波动。
现在，把 $\lambda$ 积掉（marginalize out），求 $X$ 的边缘分布：
$$P(X = k) = \int_0^\infty P(X=k|\lambda) \cdot P(\lambda) \, d\lambda$$这个积分的结果是：
$$P(X = k) = \frac{\Gamma(k+\alpha)}{k! \cdot \Gamma(\alpha)} \left(\frac{\beta}{1+\beta}\right)^\alpha \left(\frac{1}{1+\beta}\right)^k$$而这——恰好就是负二项分布的参数化形式。
通常NB写作NB(r, p)，这里的等价关系是：$r = \alpha$，$p = \beta/(1+\beta)$。或者用均值-方差参数化：均值 $\mu = r(1-p)/p$，方差 $= \mu + \mu^2/r$。
关键洞察： 推断得到的方差是 $\mu + \mu^2/r$，而不是 $\mu$。这意味着方差大于均值——模型天然就包含了过度离散。而且 $r$ 越小，离散程度越大。
这就是为什么负二项分布在RNAseq分析中如此自然：它同时捕获了两层变异——细胞间的异质性（Gamma层）和细胞内的随机性（Poisson层）。
为什么不是Poisson？
如果硬要用Poisson模型，唯一的办法是假设所有细胞的 $\lambda$ 都相同。但这是不现实的。真实数据中，不同实验条件、不同组织部位、不同细胞周期阶段，基因表达率都会有差异。
用Poisson拟合同一批细胞的数据，方差会被严重低估，导致假阳性率飙升——你本来没有差异，会被误判为有差异。
负二项分布通过引入额外的参数r（有时称为"dispersion parameter"），给了模型足够的灵活性来适应这种异质性。
和"等待时间"没有任何关系
经典的负二项分布教科书会告诉你：NB描述的是"进行直到第r次成功为止所需的试验次数"。比如抛硬币，抛到第3次正面朝上时总共抛了多少次。
这个定义在工业质量控制、排队论等领域确实有用。但——它和我们刚才的Gamma-Poisson推导毫无关系。
这两条路只是恰好得到同一个数学形式而已。
在概率论中，很多不同的随机过程确实会收敛到相同的分布——这叫"分布的等价性"。NB就是一个例子：从"等待时间"角度可以得到NB，从Gamma-Poisson层级模型也可以得到NB。但因果链完全不同。
在RNAseq的语境下，用等待时间来理解NB只会造成困惑（“我们在等待什么？"），而用Gamma-Poisson框架则直观得多：不同细胞有不同的表达率（Gamma），每个细胞内部的计数是随机的（Poisson），两者叠加产生NB。
所以下次你看到有人说"NB的 dispersion parameter r 代表自由度"或者扯什么"等待时间”，请保持警惕。那是另一个视角，和我们的生物学应用没有直接联系。
写在最后
统计模型的意义取决于它在说什么故事。
负二项分布在数学上有多种等价的定义/推导，但它们不是等价的直觉。对于RNAseq，Gamma-Poisson的故事最贴合实际——我们确实在测量一群异质性的细胞，细胞间确实有变异，而Poisson确实描述不了这种变异。
另一个有意思的点：当我们把dispersion parameter r估计出来之后，实际上是在推断这个群体有多异质。如果r很大，方差趋近于均值，NB退化为Poisson——说明细胞间比较均匀。如果r很小，方差远大于均值——说明异质性很强。这在单细胞数据里尤其有意义。
所以下次你跑DESeq2或者edgeR的时候，别忘了底层在讲什么故事。
附录：完整推导
第一步：明确我们要算什么。
我们有两层模型：
- 第一层：$\lambda$ 服从 Gamma 分布，即 $\lambda \sim \text{Gamma}(\alpha, \beta)$
- 第二层：给定 $\lambda$，计数 $X$ 服从 Poisson 分布，即 $X | \lambda \sim \text{Poisson}(\lambda)$
我们要找的是 $X$ 的边缘分布，也就是把 $\lambda$ “积掉”：
$$P(X = k) = \int_0^\infty P(X=k|\lambda) \cdot P(\lambda) \, d\lambda$$第二步：写出两个分布的具体形式。
Poisson 分布（给定 $\lambda$ 时 $X$ 的条件分布）：
$$P(X = k | \lambda) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, ...$$Gamma 分布（$\lambda$ 的分布）：
$$P(\lambda) = \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta\lambda}, \quad \lambda > 0$$这里 $\Gamma(\alpha)$ 是 Gamma 函数，它是阶乘的推广（当 $\alpha$ 是正整数时，$\Gamma(\alpha) = (\alpha-1)!$）。
第三步：代入积分式。
把上面两个式子代入边缘分布的公式：
$$P(X = k) = \int_0^\infty \frac{\lambda^k e^{-\lambda}}{k!} \cdot \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta\lambda} \, d\lambda$$第四步：整理被积函数。
把常数项提到积分外面，把含 $\lambda$ 的项合并：
$$P(X = k) = \frac{\beta^\alpha}{k! \cdot \Gamma(\alpha)} \int_0^\infty \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda$$注意：
- $e^{-\lambda} \cdot e^{-\beta\lambda} = e^{-\lambda(1+\beta)}$
- $\lambda^k \cdot \lambda^{\alpha-1} = \lambda^{k+\alpha-1}$
第五步：关键的换元技巧。
观察积分里面的形式：$\lambda^{\text{某个幂次}} \cdot e^{-\lambda \cdot \text{某个常数}}$。这看起来很像 Gamma 分布的密度函数！
回忆 Gamma 分布的密度：$\frac{\beta^{\alpha}}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta\lambda}$，其中 $\alpha$ 和 $\beta$ 是形状和速率参数
我们的积分缺少前面的归一化常数。但我们可以构造它！
令 $\alpha' = k + \alpha$，$\beta' = 1 + \beta$。那么：
$$\int_0^\infty \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda = \frac{\Gamma(k+\alpha)}{(1+\beta)^{k+\alpha}} \cdot \underbrace{\int_0^\infty \frac{(1+\beta)^{k+\alpha}}{\Gamma(k+\alpha)} \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda}_{= 1}$$积分里面的部分正好是 $\text{Gamma}(k+\alpha, 1+\beta)$ 分布的密度函数，积分结果等于 1。
所以：
$$\int_0^\infty \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda = \frac{\Gamma(k+\alpha)}{(1+\beta)^{k+\alpha}}$$第六步：代回原式。
$$P(X = k) = \frac{\beta^\alpha}{k! \cdot \Gamma(\alpha)} \cdot \frac{\Gamma(k+\alpha)}{(1+\beta)^{k+\alpha}}$$整理一下：
$$P(X = k) = \frac{\Gamma(k+\alpha)}{k! \cdot \Gamma(\alpha)} \cdot \frac{\beta^\alpha}{(1+\beta)^{k+\alpha}}$$$$= \frac{\Gamma(k+\alpha)}{k! \cdot \Gamma(\alpha)} \cdot \left(\frac{\beta}{1+\beta}\right)^\alpha \cdot \left(\frac{1}{1+\beta}\right)^k$$第七步：与负二项分布对比。
标准负二项分布 NB(r, p) 的概率质量函数是：
$$P(Y = k) = \frac{\Gamma(k+r)}{k! \cdot \Gamma(r)} (1-p)^r p^k, \quad k = 0, 1, 2, ...$$对比发现：
- $r = \alpha$（形状参数）
- $p = 1/(1+\beta)$（成功概率）
- $1-p = \beta/(1+\beta)$
完全匹配！
第八步：验证均值和方差（可选但推荐）。
均值：
利用全期望公式 $E[X] = E[E[X|\lambda]] = E[\lambda] = \alpha/\beta$
方差：
利用全方差公式
$$Var(X) = E[Var(X|\lambda)] + Var(E[X|\lambda]) = E[\lambda] + Var(\lambda) = \alpha/\beta + \alpha/\beta^2 = \alpha(1+\beta)/\beta^2$$用负二项分布的参数表示：
- 均值 $\mu = r(1-p)/p$
- 方差 $= \mu + \mu^2/r$
从我们的推导：
- $\mu = \alpha/\beta$
- $Var = \alpha/\beta + \alpha/\beta^2 = \mu + \mu^2/\alpha$
所以 $r = \alpha$，即 Gamma 分布的形状参数就是 NB 的 dispersion 参数。
整个推导的核心思路是：
- 写出层级模型的两个分布
- 对隐藏变量 $\lambda$ 积分（边缘化）
- 利用 Gamma 函数的积分性质完成计算
- 整理结果，发现恰好是负二项分布的形式
这个推导展示了为什么负二项分布能自然地描述 RNAseq 数据：它本质上是 Poisson 分布的"加权平均"，权重由 Gamma 分布给出。当细胞间异质性越大（Gamma 分布的方差越大），最终分布的 overdispersion 就越严重。