【原】原理篇：HaplotypeCaller变异检测

BIOINFO_J 2020-07-10

展开全文

GATK的HaplotypeCaller模块(以下简称HC)是用来进行snp及小indel变异检测的，以下是基本的原理。

1. 总览

HC使用的是预组装的方法，能提高变异检测的准确度，但是在某种程度上增加了资源的消耗和分析时长。在分析时，该模块并不会在基因组范围进行全局的变异检测，而是划定高变区间检测，以下是具体的步骤及原理图。

1. 定义变异区间：首先寻找确定进行变异检测的区间
2. 数据组装：对区间内的数据及基因组进行组装(De Bruijnlike graph)并预估单倍型
3. 根据比对结果，确定单倍型的最大似然值：将reads回帖单倍型数据(PairHMM algorithm)，计算各个单倍型出现的最大似然值
4. 将单倍型比对结果转化为在特定位置位置处的分型

左右滑动查看全部内容。

2. 模块分步分析说明

2.1 定义区间

为了加速程序的处理速度，gatk并不会对全基因组范围内的所有位点进行变异检测，只会对高变区进行检测，也就是所谓的Active Regions。此处定义的区间选择，来源于比对结果，包括其中的错配、插入、缺失以及soft clipping。

计算时，会计算每个位点出现突变的概率，这一步会同时考虑分型最大似然值以及杂合率，之后通过高斯核做卷积来扩大概率值。当概率值超过0.002时，该区间就会被定义为高变区。在计算过程中，程序会对毗邻的高变区(100bp)进行合并处理，设定的区域大小最小50bp，最长300bp，如果合并区超过300bp，会被截断形成2个或多个高变区，并且存在overlap的reads在相邻的两个高变区中都参与运算。

2.2 数据组装

高变区的reads以及对应的基因组区域会被切割成存在overlap的小片段，若参考基因组对应的片段集合中存在重复，会将短片段的长度递增直到没有重复或者达到最大长度限制(65nt)，若存在65nt的重复短序列，组装停止。序列切割完成后使用de-Bruijn-like graphs算法进行组装。组装获得的边缘根据比对到的reads数目分配权重，比对数量较少的边会从图形中删除，以此来对图形进行精简并形成单倍型数据。

组装完成后会给每个单倍型生成CIGRA信息，以此为依据进行位点的分型。

2.3 隐马尔科夫模型预估最大似然值

组装完成后，将原有的reads回帖到获得的单倍型序列，之后根据比对结果预估在不同单倍型中的最大似然值(likelihood)。具体原理是使用隐马尔科夫模型(pair-HMM，paired Hidden Markov Model)对reads和单倍型之间进行最大似然值的预估，pair-hmm模型原理如下图所示。

其中δ表示gap发生罚分，ε表示gap延伸罚分，涉及M(match，得分会考虑比对质量)、I(insertion)以及D(deletion)三种类型，具体计算如下。

M_ij =P(r_i|h_j, q_i) (M_i−1,j−1T_MM + I_i−1,j−1T_IM + D_i−1,j−1T_DM)
I_ij =M_i−1,jT_MI + I_i−1,jT_II
D_ij =M_i,j−1T_MD + D_i,j−1T_DD

其中i表示reads上碱基的位置，j表示单倍型中碱基的位置，T表示相应的罚分，具体对应关系见上图。

2.4 分型结果判定

根据比对之后的cigra信息进行分型结果的判定。判定完成后，会与各单倍型进行比对分析，以上一步最大似然值为依据，计算分型的最大似然值，公式如下。

以上是获取某一分型出现的最大似然值。根据得到的最大似然值，可以计算分型出现的后验概率，公式如下。

之后，通过比较不同分型出现的概率来确定对应样本在某一位点处的分型结果。

3. 猜你喜欢

[1] 动植物重测序分析
[2] gatk进行变异检测
[3] gvcf文件信息及处理

4. 参考信息：

[1] https://gatk./hc/en-us/articles/360035531412
[2] https://gatk./hc/en-us/articles/360035531812-GVCF-Genomic-Variant-Call-Format
[3] Scaling accurate genetic variant discovery to tens of thousands of samples. doi: https:///10.1101/201178