在各类基因组高通量测序的研究软件当中，GATK（Genome analysis toolkit）可谓是当红大明星，尤其是GATK的haplotypecaller组件，现在已是重测序变异挖掘的最优选择。虽然相比于其他用于变异识别的软件来说haplotypecaller所用的时间相对来说更长一些，但并不妨碍它始终占据最佳变异识别软件的宝座。
GATK haplotypecaller“吃的是bam，挤出的是vcf”，将序列比对文件转化成标准变异格式文件，可以说非常省心了。并且vcf文件中包含的内容非常详尽，包括位点的位置、变异前后的碱基、位点的基因型、变异的质量值等等，我们可以方便地利用它进行后续的各种分析。不过有的时候，这些信息也会让我们有小小的迷惑。对于一个个体的变异信息来说，大多数位点要么只有一种reads支持（MAF=0，纯合突变/无突变），要么两种reads大致各占一半（MAF=0.5，杂合突变）。但是有时候，我们在仔细研究vcf文件的AD值的时候，会发现有一些位点明明有两种reads，其中一种reads的占比却比较低（MAF值比较低但不为0）。并且更加奇妙的是，可能MAF值差不多的两个位点，haplotypecaller给出的基因型却不一样。例如下图中的两个位点，一个MAF值约为0.048，被判定为0/0（无突变），另一个MAF值约为0.064，被判定为0/1（杂合突变）。好奇的你有没有想过，这到底是为什么？haplotypecaller是怎么计算出这些位点的基因型的？
两个位点的GT和AD信息（来源于真实数据）
其实haplotypecaller的工作原理正像它的名字一样，是通过识别单倍型（haplotype）来确定各个位点的基因型的。其主要工作流程包含4个主要的步骤：识别高变区、组装单倍型、计算似然值、分配基因型。下面我们结合GATK官方给出的流程图来一探究竟。
haplotypecaller主要针对高变区进行单倍型识别，那么第一步就是要识别基因组中发生活跃变化的区域。这一步的运算通过滑窗方法进行，计算一个区域内的碱基错配、插入缺失以及soft-clip带来的熵值，并将超过阈值的区域提取出来进行后续计算。在得到高变区之后，haplotypecaller会在这个区域执行一个局部的重比对，基于图论的方式来组装出最可能存在的单倍型，并使用Smith-Walterman（史密斯-沃特曼）算法将单倍型重新比对到参考基因组上，这一步最终得到了候选的单倍型以及可能的突变位点。（从图解中我们可以看出来，可能性很低的单倍型会被过滤掉，否则后续的计算量会成指数增长）第三步，haplotypecaller会将reads和单倍型两两配对，使用PairHMM（配对隐马尔科夫模型）来给配对结果打分。在这一步中，碱基测序质量是一个非常重要的参数。这一步会输出一个矩阵，包含每个read对应每种单倍型的似然值。之后程序会将矩阵进行转化，将单倍型的似然值转化为等位基因的似然值。得到了每个read对应每个等位基因的似然值之后，程序会据此计算出每种基因型对应的概率，并将最可能的基因型分配给这个位点。具体的计算方法在下图中有详细的介绍，我们可以通过示例看到单倍型的似然值是如何通过运算变换成基因型的概率的。现在，位点是否突变以及对应的基因型就已经计算出来了。那么运算是不是就可以结束了呢？当然不是，我们还需要对基因型的质量进行评价。我们会在vcf文件中看到突变位点的PL和GQ两个信息，这两个信息通常会用于对突变基因型的质量进行评价。PL指基因型对应的经过标准化的Phred-scaled probability，每一个可能的基因型都有一个对应的PL值，PL值为0的是最优基因型。而GQ值取的是第二低的PL值与99之间较低的数字，可能说起来比较拗口，但是仔细一想，我们就可以发现这个值可以用来判断该位点是其他基因型的可能性有多大。GQ值越高，位点是其他基因型的可能性就越小，即给出基因型的质量越高。综上所述，我们可以看出来，GATK haplotypecaller进行变异识别的算法过程确实非常缜密，不仅用到了reads对位点的支持率，还考虑了周围区域组成的单倍型，以及碱基的测序质量。这下我们拿到重测序输出的变异信息的时候，是不是心里更加有底了呢？