变异检测（三）之点突变检测

点突变和插入删除突变是比较常见的变异类型。所谓点突变就是简单的单个碱基的替换，插入删除突变是一个及一个以上碱基的删除或插入，在比对上就是会引入空位。这样的突变的检测，理论上来说在比对结果中可以被识别。但是由于背景噪声的存在，经常需要应用统计学的方法进行去噪声，然后才能对真正的变异进行识别。这样的统计学方法在处理体细胞突变的检测问题中显的尤为重要。因为体细胞突变的突变频率有时往往很低，隐藏在背景噪声中，很难被发现。为此，科学家开发了很多针对性的软件解决这个问题，如Mutect，VarScan等等。近年来，也有科学家将深度学习的算法引入到突变的检测识别方法中来，主要是由于基因数据积累的越来越多，可以通过深度学习的方法将真实突变从复杂的噪声中识别出来，这样的方法也取得了比较好的效果，比如GATK4中就引入了深度学习的方法进行变异位点的过滤。

针对插入删除变异的检测识别问题，同样存在一个难点，那就是比对的问题。复杂的插入删除突变有时在比对会与原序列出现很大差异，这给比对造成很大困难，而造成检测识别的失败。针对这个问题，科学家也开发了很多针对性的算法，如重比对，局部组装等等。

点突变的检测

简单的来说，点突变的检测主要是一个分类问题，以体细胞突变的检测识别为例，这个检测识别主要是区分真实突变与背景噪声，其中背景噪声包括“比对错误”，“测序错误”，“遗传性变异”等等。不同的噪声有其一定的特征，比如可以根据测序质量及链偏向性去估计测序错误的可能性。其实，这就是一个统计学上的分类问题。既然是一个分类问题，那么就可以使用很多统计学的分类算法去解决它。下面，我举几个例子，他们是比较著名的软件中使用的统计学方法，主要会介绍他们的算法，并介绍他们选择的特征。

很多软件都使用的是贝叶斯算法。下面，我就以一个经典的体细胞变异检测的软件的算法为例，简单介绍下贝叶斯公式在变异检测中的应用。
所谓的贝叶斯公式是基于条件概率的统计学方法。贝叶斯公式中，事件Bi的概率为P(Bi)，事件Bi已发生条件下事件A的概率为P(A│Bi)，事件A发生条件下事件Bi的概率为P(Bi│A)。
MuTect软件包含了两个贝叶斯分类器，第一个是区分该位置是否发生突变，第二个是区分该变异是否为体细胞变异。

第一个贝叶斯分类器

下面，我们首先介绍第一个贝叶斯分类器。首先，我们设定是某个位置，然后开始使用这个分类器去区分这个位置是否发生突变。
r表示参考基因组：{A,T,C,G}
b表示该位点的碱基：{A,T,C,G}
e表示该位点的质量，也就是发生测序错误的可能性，它的计算方法与该碱基的测序质量Q有关。e等于10的-q/10次方。
假设了两个模型：
M0:该位点没有发生突变，所有与参考基因组不一样的碱基都是测序错误；
Mf：该位点发生了突变，且突变频率是f，M0其实等于f=0的Mf。

下面我可以计算M的似然：
LOD(M) = P(bi|ei,r,m,f)的积
其中每一个bi都是每i个reads比对到该位置碱基，ei是该碱基的质量。计算出得到P，具体计算方法如下：

如果bi是参考基因组的碱基类型，那么

P = f * ei/3 + (1-f)(1-ei)

意义: 因为该位点的突变频率是f，所以说明有f的碱基是非参考基因组的碱基类型，但bi是参考基因组的碱基类型，所以f的碱基是测序错误，而其他的1-f的碱基是没有测序错误的，其中ei的概率是测错的，那么某个碱基类型的测错的概率就是ei/3.

如果bi不是突变的碱基类型，那么

P = (1-f) * ei /3 + f(1-ei)

如果bi是其他的情况，那么：

P = ei / 3

m是突变的碱基类型，也就是有除了参考基因组碱基类型外还有三种碱基类型。需要比较两个模型的似然，才能识别该位点是否发生突变：

LOD(m,f) = log10(L(Mf)P(m,f)/L(M0)(1-P(m,f)))

如果LOD(m,f)大于log10Qt，那么就认为m是一个突变碱基。Qt在MuTect最初的文献中被设定为6.3，表示该位点突变的置信度是背景噪声的6.3倍以上。根据以上的公式，现在可以区分该位点是否存在突变。其中P(m,f)与f的值还未知，P(m,f)可以通过P(m)和P(f)计算得到，假设两者是统计学上独立的，P(f)是均一分布，P(f) = 1,P(m)是一个经验值，文献中将它设置为 3 × 10−6。
针对f的值，可以通过最大似然估计计算得到，也可以通过突变碱基数处理总测序深度得到。

第二个贝叶斯分类器

第二个贝叶斯分类器是用于区分遗传性突变和体细胞突变的。基本的模型与第一个分类器类似，其中f被设置为0.5，因为理论上遗传性突变是杂合的。

LOD(N) = log10(L(M0)P(m,f)/L(M0.5)P(germline))

就是只有大于log10Qn时才被认为是体细胞突变。
其中P(germline)和Qn阈值设定如下：

P(germline|non-dbSNP site) = 5 × 10−5 Qn|non–dbSNP site = 2.2

P(germ line|dbSNP site) = 0.095 θn|dbSNP site = 5.5