变异检测

简介

当我们利用实验技术对于待检测的目标（动植物、微生物）进行了测序，通过获得的数据经过比对后。作为分析，特别是针对于DNA的基因组分析中至关重要的环节就是变异的检测。而变异检测的方法从诞生之初如同生物进化般迅速的发展。从最开始的所见即所得，到现在所得不一定所见的检测方法呈现一种急速的发展状态。下面就从几类较为常用的检测算法中介绍。

NGS数据的变异检测 目前已经有大量的方法及工具利用NGS数据进行单碱基变异的检测。这些方法不仅仅是计算工具，而他们本身又是和试验方法密不可分的。随着NGS数据越来越多，这些工具也变得越来越重要，利用度也越来越广。

一般情况下，检测的方法都基于以下5个方面：

1. 过滤read，从源头去除一些偏好及错误
2. 将read比对到参考基因组
3. 利用统计或者其他算法，预测各个位点发生变异的可能性
4. 由于不同算法有着不同的优势及劣势，所以预测出的结果一般会进行部分过滤
5. 注释

而下面将从遗传变异检测算法 和体细胞变异检测算法来进行详细介绍。

遗传变异检测算法

有很多的SNV检测算法都是用来检测基因组中的遗传变异。 这些遗传变异表示这些变异是继承自父母的变异。一般情况来讲，这些变异发生有一定的人群频率（即便有些发生的频率很低），这些变异我们会称之为单核苷酸多态（SNP）。一般来说，SNP仅仅被用来描述这类变异，但是在文章中又常常会被用来形容同义体单核苷酸突变（未造成氨基酸变异的突变）。而对于这类变异的检测主要有下面介绍的这几种工具及算法。

SAMtools+BCFtools

SAMtools是一款强大的SNV检测工具。它是一个可以从比对结果（SAM/BAM格式）中发现SNV的软件，他有许多可调节的参数。SAMtools可以通过考察覆盖度、包含参考碱基及变异碱基read数来评价并过滤SNP。他的检测原理很像是将变异read展现到人的眼前，然后使用硬性的指标来进行判断是否属于变异。这种方式就像是人眼能看到，就算检测到的一种方式。较为依赖于比对算法，并且对于一些较长的插入/缺失无法进行检测。

下面就是一个使用SAMtools配合BCFtools来进行检测的一个简单例子。

samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -bvcg - > var.raw.bcf  
bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf

UnifiedGenotyper

除了SAMtools还有一种广为应用的软件叫做GATK（Genome Analysis Toolkit），而GATK这个套件中有专门用来检测遗传突变的就是UnifiedGenotyper这个方法。他主要基于的就是贝叶斯理论的统计模型来判断位点是否为真实变异。

他利用的对于同一个位置变异及非变异的碱基数目构建贝叶斯统计模型。

$P(G|D)={\frac{P(D|G)P(G)}{P(D)}}={\frac{P(D|G)P(G)}{\sum\limits_{i=1}^{n}P(D|G_{i})P(G_{i})}}$ 上式中：

• $D$表示观察到的数据，也就是比对的reads
• $G$表示被计算的基因型
• $G{i}$表示一共有_n种基因型的第i种基因型

整个公式理解起来就是，在我们的观测到的数据 $D$ 的条件下，基因型为 $G$ 的可能性。而假设的一般情况 $G$（基因型）有 $AA$（二倍体，未发生突变）， $AB$ （杂合突变）， $BB$ （纯合突变），这三种情况。而根据条件概率我们可以将 $P(G|D)$ 转化为 $\frac{P(D|G)P(G)}{P(D)}$ ，也就是在我们假设的某一种基因型下，得到我们观察到的数据的可能性最大的那个基因型就是我们的预测结果。

HaplotypeCaller

而GATK中的UnifiedGenotyper方法只是解决了部分的预测问题，也就是对结果有一个准确性的评估，对于结果给出的贝叶斯估计的概率，但是由于比对所造成的只能检测短小的插入/或缺失并不能解决。所以在此引入了HaplotypeCaller。HaplotypeCaller的计算方式与UnifiedGenotyper相同，只是在预测的前一阶段加入了de novo。对于比对的结果会进行组装，通过组装出来的结果再重新比对，这样如果有大片段的插入/缺失就能够检测出来。

DeepVariant

而在近几年，随着深度学习的发展，已经有了像DeepVariant这样基于深度学习的检测工具出现。DeepVariant主要依据的模型是CNN（卷积神经网络），CNN通常用在图像类的识别中。而DeepVariant开发团队非常新颖的将其利用到变异检测中，通过对于标准数据的学习，训练模型，最终使得模型能够很好的进行预测。DeepVariant流程图如下图所示。

图来自于Creating a universal SNP and small indel variant caller with deep neural networks

体细胞变异检测算法

随着癌症基因组学的发展，检测体细胞变异成为描述癌症基因组的一个至关重要的环节。但是由于体细胞变异与遗传变异有些重要差异，所以并不能直接使用遗传变异检测算法检测体细胞变异检测。例如基于贝叶斯统计模型的检测算法，对于遗传突变我们可以假设其基因型是三种，野生型、杂合突变、纯合突变，而这三种基因型反映到数据上就是变异的read比例大约在0%、50%、100%。然而对于体细胞变异则变化更多，read的比例可以从0.01%到100%。而且对于体细胞变异的检测，需要达到非常高的准确性还需要对照样本。判断一个变异是否为体细胞变异不能仅仅从突变频率进行判断。更准确的应该是确定这个变异不存在其他正常细胞中，而只是存在于这些特殊样本中。下面就来介绍几种常用的体细胞变异检测算法。

MuTect

MuTect就是一种利用贝叶斯分类器来进行体细胞变异检测的工具。对于基因组的每个位点的原本碱基（参考基因组）用 r ∈{A,C,G,T} 来表示，$bi$表示覆盖这个位点的第 _i 条read在这个位置的碱基，$ei$ 表示这个位点的错误可能性（每个碱基都有碱基质量值表示他测序错误的可能性 $e_i=10^{-\frac{q_i}{10}}$ ）。MuTect用两种模型来描述数据，第一：model $M_0$，该位点没有发生变异，所看到到的不同的碱基都是由于测序错误导致；第二：model $M^m_f$ 描述 _m 变异确实存在，并且频率为 f。然后构建第二种模型的似然率，公式如下：

$L(M^m_f)=P(\{b_i\}|\{e_i\},r,m,f) =\prod^d_{i=1}P(b_i|e_i,r,m,f)$ 这公式假设了，不同read之间的测序错误是独立的（所以联合概率才能表示为单独概率的乘积）。如果所有的替换错误都是相等，那么可得到下面的公式：

$P(b_i|e_i,r,m,f)=\left\{ \begin{array}{rcl} f*{e_i/3}+(1-f)(1-e_i) & & {\text{if }b_i=r}\\ f(1-e_i)+(1-f)*e_i/3 & & {\text{if }b_i=m}\\ e_i/3 & & {\text{otherwise}} \end{array} \right.$ 变异的检测需要同时计算两种模型的似然率，然后计算 LOD score，计算公式如下：

$LOD_T(m,f)=\log_{10}{(\frac{L(M^m_f)P(m,f)}{L(M_0)(1-P(m,f))})} \geq log_{10}{\theta_{T}}$ 通过此模型MuTect达到检测变异的目的，然而要彻底获得体细胞变异，还需要有对照样本进行检测，如果同一个位点在对照中也检测出来，则说明，这个位点更可能是遗传变异而不是体细胞变异。

Varscan

而目前使用较多的除了GATK中的Mutect以外，还有Varscan。Varscan与其他算法有一个特别大的差异，就是他的输入文件并不是普通的比对结果（SAM/BAM），而是通过SAMtools处理过的pileup格式文件。而且必须同时输入待测样本及对照样本（Mutect可以使用单样本，但是会有其他需求来过滤非体细胞变异），然后通过以下几步来进行预测：

1. 待测样本及对照样本均达到覆盖度要求（默认值：至少三条read的碱基质量>=20）；
2. 根据观察到的base数据确定样本基因型（默认值：至少2条read支持变异，频率>=8%）；
3. 没有样本达到上述要求，则认为该位点为野生型； 如果变异的频率>=75%，则认为该位点为纯合变异；其他情况均认为为杂合变异；
4. 对于至少一个样本（待测样本或对照样本）检测出变异的位点进行样本之间的比较，参照如下规则：
   1. 如果基因型不一致，则利用read数来进行单尾fisher精确检验。如果 P 值显著（默认为0.1），则认为是体细胞变异或杂合缺失（前提是对照样本为杂合）否则则为遗传变异。
   2. 如果基因型一致，则变异直接认为是遗传变异。