r - 错误概率函数

Question

我的 DNA 扩增子可能在 PCR 扩增过程中出现碱基错配。我感兴趣的是，给定每个碱基的错误率、不匹配的数量和扩增子中的碱基数量，序列包含错误的概率是多少。

我看到一篇文章[Cummings, S. M. 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 个独特序列与最相似的序列各有 3 bp 差异。由三个独立的 PCR 错误生成新序列的概率等于 P = 0.0011。

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

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

我在实现过程中可能做错了什么？我是否误解了什么？

谢谢

Answer 1

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

您想要进行的计算是:

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

即假设 30 次扩增，每个扩增的出错几率为 1.85e-5，那么在 800 个独立碱基中看到少于 3 个修饰的概率是多少？ IE。您正在计算 30 次不正确的概率。

有点统计，可能值得一动……

更多地考虑这一点，当在这里处理非常小的概率时，您将开始看到浮点不准确。 IE。当 x 的绝对值小于大约 1e-10 时，1-x(其中 x 是一个小数字)就会开始出错。此时使用对数概率是一个好主意，特别是 log1p 函数很有帮助。使用:

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

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