点突变的假设检验

在点突变检测中，可以用假设检验的方法来确定某个位点是真实突变，还是仅仅是测序错误。例如某个位点覆盖有100X，它有10个位点是突变碱基，测序错误率为0.01，为此做的检测检验：

  • H0：该位点突变碱基频率小于等于0.01（突变碱基源于测序错误）
  • H1：该位点突变碱基频率大于0.01

使用R进行显著性计算：

  1. binom.test(10, 100, 0.01, alternative = "greater")
  3. ###### output ######
  4.                         Exact binomial test
  6. data:  10 and 100
  7. number of successes = 10, number of trials = 100, p-value =
  8. 7.632e-08
  9. alternative hypothesis: true probability of success is greater than 0.01
  10. 95 percent confidence interval:
  11.  0.05526324 1.00000000
  12. sample estimates:
  13. probability of success 
  14.                    0.1

这里假设数据符合二项分布，结果p值<0.001，所以认为这个位点是发生了突变的。当然实际情况并没有那么简单，所以才有各种各样的软件开发出来解决这个问题，如：gatk，varscan，strelka，freebayes等。

检测需要的深度

相对应的有一个问题：“如果需要能够检测到频率为0.02，需要测多少的深度？”这通常用统计学上的功效分析来解决。首先需要明白几个概念。

  1. I型错误：当H0为真却拒绝了H0
  2. II型错误：当H0为假却没有决绝H0

功效分析通常关注四个量：

  • 样本大小：观测的数目
  • 显著性水平（alpha）：即I型错误的概率
  • 功效（power）：1 - II型错误的概率
  • 效应值：与希望检测的原假设之间的最小偏差

这四个量知道了其中三个可计算出第四个。很明显，我们要知道的深度对应的是样本大小，这是个未知量。另外三个：显著性水平一般设为0.05，功效可以是0.8-0.99，效应值是检测下限0.02。功效分析依赖于假设检验中使用的统计方法，这里我们假设突变符合参数为突变频率的二项分布，使用伯努利随机过程来建模。使用R进行的功效分析如下：

  1. # sample size: ?
  2. # alpha: 0.05
  3. # power: 0.98
  4. # effect size: 0.02
  5. library(binom)
  6. cloglog.sample.size(p.alt = 0.02, p = 0.01, power = 0.98, alpha = 0.05)
  8. ###### output ######
  9.     p.null p.alt delta alpha power    n phi
  10. 1   0.01  0.02  0.01  0.05  0.98 2352   1

求得的深度是2352x。这意味着测序深度为2352X的时候，大于等于0.02频率的突变几乎都可以检测出来。

讨论

  1. 如果测序错误率为0.01，就不可能检测到频率0.01的突变
  2. 如果真实突变频率是0.1，观测到的突变频率不可能精确的等于真实突变频率，计算95%的置信区间：
    1. 100X 对应（0.05-0.176）
    2. 200X 对应（0.062-0.150）
    3. 500X 对应(0.075, 0.130)

深度越高，频率越准确

  1. 应该通过实际检测来证实理论结果的正确性