r - 错误概率函数

Question

我有在 PCR 扩增过程中可能出现碱基错配的 DNA 扩增子。我的兴趣是，考虑到每个碱基的错误率、错配数和扩增子中的碱基数，序列包含错误的概率是多少。

我遇到了一篇文章[Cummings, SM et al (2010)。群体遗传分析中 PCR、克隆和测序错误的解决方案。保护遗传学，11（3），1095-1097。doi:10.1007/s10592-009-9864-6] 提出了这个公式来计算这种情况下的概率质量函数。

在此处输入图像描述

我用 R 实现了公式，如下所示

pcr.prob <- function(k,N,eps){
  v = numeric(k)
  for(i in 1:k) {
    v[i] = choose(N,k-i) * (eps^(k-i)) * (1 - eps)^(N-(k-i))
    }
  1 - sum(v)
}

从文章中，建议我们使用 30 个循环的 PCR 分析一个 800 bp 的扩增子，1.85e10-5每个循环每个碱基的错误掺入，并发现10每个3bp 与其最相似的序列不同的独特序列。由三个独立的 PCR 错误产生新序列的概率等于P = 0.0011。

但是，当我使用公式的实现时，我得到了不同的值。

pcr.prob(3,800,0.0000185)
[1] 5.323567e-07

在我的实施中我可能做错了什么？我误解了什么吗？

谢谢

score 2 · Accepted Answer

我认为他们得到了正确的数字（0.00113），但在他们的论文中解释得很糟糕。

您要做的计算是：

pbinom(3, 800, 1-(1-1.85e-5)^30, lower=FALSE)

即在 800 个独立碱基中看到少于三个修饰的概率是多少，给定 30 个扩增，每个有 1.85e-5 出错的机会。即你正在计算它不保持正确 30 次的概率。

有点统计，可能值得一试……</p>

仔细考虑这一点，当在这里使用非常小的概率时，您将开始看到浮点不准确。即a 1-xwhere xis a small number 将在绝对值x小于约1e-10 时开始出错。在这一点上，使用对数概率是一个好主意，特别是这个log1p函数很有帮助。使用：

pbinom(3, 800, 1-exp(log1p(-1.85e-5)*30), lower=FALSE)

即使错误合并率非常低，也将继续工作。

1 回答 1