Fork me on GitHub

变异检测(三)之点突变检测

引言

点突变和插入删除突变是比较常见的变异类型。所谓点突变就是简单的单个碱基的替换,插入删除突变是一个及一个以上碱基的删除或插入,在比对上就是会引入空位。这样的突变的检测,理论上来说在比对结果中可以被识别。但是由于背景噪声的存在,经常需要应用统计学的方法进行去噪声,然后才能对真正的变异进行识别。这样的统计学方法在处理体细胞突变的检测问题中显的尤为重要。因为体细胞突变的突变频率有时往往很低,隐藏在背景噪声中,很难被发现。为此,科学家开发了很多针对性的软件解决这个问题,如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