构建 DNA 测序的负二项分布

机器算法验证 泊松分布 负二项分布 生物信息学 序列分析
2022-01-25 21:46:27

负二项分布已成为生物信息学中计数数据(特别是来自给定实验的基因组给定区域内的预期测序读数)的流行模型。解释不一:

  • 一些人将其解释为类似于泊松分布但有一个附加参数,允许更多自由地模拟真实分布,方差不一定等于平均值
  • 一些人将其解释为泊松分布的加权混合(在泊松参数上具有伽马混合分布)

有没有办法将这些基本原理与负二项分布的传统定义相匹配,即在看到一定数量的失败之前对伯努利试验的成功次数进行建模?或者我应该认为这是一个快乐的巧合,泊松分布与伽马混合分布的加权混合具有与负二项式相同的概率质量函数?

4个回答

恕我直言,我真的认为使用负二项分布是为了方便。

所以在 RNA Seq 中有一个普遍的假设,如果你在无限次重复中对同一个基因进行无限次测量,那么真实分布将是对数正态分布。然后通过泊松过程(带有计数)对该分布进行采样,因此每个基因在重复中的真实分布读数将是泊松对数正态分布。

但是在我们使用的包中,例如 EdgeR 和 DESeq,这个分布被建模为负二项分布。这并不是因为编写它的人不知道泊松对数正态分布。

这是因为泊松对数正态分布是一件很糟糕的事情,因为它需要数值积分来进行拟合等。所以当你实际尝试使用它时,有时性能真的很差。

负二项分布具有封闭形式,因此更容易使用,并且伽马分布(基础分布)看起来很像对数正态分布,因为它有时看起来有点正态,有时有尾巴。

但是在这个例子中(如果你相信这个假设)它不可能在理论上是正确的,因为理论上正确的分布是泊松对数正态分布,这两个分布是彼此的合理近似,但并不等价。

但我仍然认为“不正确”的负二项分布通常是更好的选择,因为从经验上看,它会给出更好的结果,因为积分执行缓慢并且拟合可能表现不佳,尤其是对于长尾分布。

我浏览了几个网页并找不到解释,但我想出了一个用于整数值的解释。假设我们有两个放射源,分别以 α 和 β 的速率独立产生 αβrαβ

个 β 粒子之前,α 粒子的数量分布是什么?r

  1. 将 α 粒子视为成功,将 β 粒子视为失败。当检测到一个粒子时,它是 alpha 粒子的概率是所以,这是负二项分布αα+βNB(r,αα+β)

  2. 考虑第个 beta 粒子的时间这遵循伽马分布如果以条件,则时间之前的 alpha 粒子数遵循泊松分布个 β粒子之前的 α 粒子数量分布是 Gamma 混合泊松分布。trrΓ(r,1/β).tr=λ/αtrPois(λ).r

这就解释了为什么这些分布是相等的。

我只能提供直觉,但伽马分布本身描述了(连续的)等待时间(发生罕见事件需要多长时间)。因此,离散泊松分布的伽马分布混合会导致离散等待时间(试验直到 N 次失败)这一事实似乎并不令人惊讶。我希望有人有更正式的答案。

编辑:我总是证明负二项式分布是合理的。测序如下:实际的测序步骤只是从大型分子库(泊松)中采样读取。然而,该文库是通过 PCR 由原始样本制成的。这意味着原始分子呈指数级放大。而伽马分布描述了k个独立指数分布的随机变量的总和,即在相同数量的PCR循环中扩增k个样本分子后文库中有多少分子。

因此,负二项式模型 PCR 然后测序。

我将尝试给出一个简单的机械解释,我在考虑这个问题时发现它很有用。

假设我们在文库准备之前对基因组进行了完美的统一覆盖,并且我们观察到平均覆盖一个位点的假设测序是一个挑选原始 DNA 片段的过程,将其置于一个随机过程中,该过程经历 PCR、二次采样等,并以频率从片段中得出一个碱基,否则失败。如果排序一直持续到失败,则可以使用负二项分布对其进行建模。μpμ1ppNB(μ1pp,p)

计算这个分布的矩,我们得到期望的成功次数对于成功次数的方差,我们得到 - 片段库准备失败的速率增加了观察到的覆盖率的方差。μ1ppp1p=μσ2=μ(1p)1

虽然以上是对测序过程的稍微人为的描述,并且可以对 PCR 步骤等建立适当的生成模型,但我认为它可以让我们深入了解过度分散参数的起源直接来自负二项分布。我确实更喜欢将速率集成的泊松模型作为一般解释。(1p)1