GATK项目中Pair HMM概率重比对模型详解

GATK项目中Pair HMM概率重比对模型详解

gatk Official code repository for GATK versions 4 and up 项目地址: https://gitcode.com/gh_mirrors/ga/gatk

概述

在基因组分析工具包GATK的HaplotypeCaller和Mutect工具中，Pair HMM（Pair Hidden Markov Model，配对隐马尔可夫模型）是实现高精度变异检测的核心算法之一。本文将深入解析这一概率重比对模型的技术原理和实现细节。

Pair HMM模型基础

Pair HMM模型用于计算测序reads相对于候选单倍型的比对概率，即P(R|H)，其中R代表测序read，H代表候选单倍型。这一概率是所有可能的比对路径概率之和：

P(R|H) = Σ P(R,A|H)

比对路径A被表示为一系列状态序列，每个状态包含三个要素：

• i：read中的位置
• j：单倍型中的位置
• s：状态（匹配M、插入I或缺失D）

状态转移机制

模型定义了三种基本状态及其转移规则：

1. 匹配状态(M)：read位置i+1与单倍型位置j+1匹配
2. 缺失状态(D)：单倍型位置j+1相对于read缺失
3. 插入状态(I)：read位置i+1相对于单倍型插入

状态转移概率矩阵T决定了这些状态之间的转换可能性，其中包含两个关键参数：

• δ：indel起始概率
• ε：indel延续概率

概率计算组成

比对概率由两部分组成：

1. 状态序列概率：由转移矩阵T决定
2. 碱基发射概率：由碱基质量和比对状态决定

碱基发射概率公式为：

P(b2|b1,q) = {
    ε(q)/3 (b1≠b2时)
    1-ε(q) (b1=b2时)
}

其中ε(q)=10^(-q/10)是基于Phred质量分数q计算的错误率。

动态规划实现

Pair HMM通过动态规划算法高效实现，定义三个矩阵M、I、D分别表示以匹配、插入或缺失状态结束的所有路径的总似然。

核心递推公式为：

M[i,j] = P(r_i|h_j,q_i) × (M[i-1,j-1]×T_MM + I[i-1,j-1]×T_IM + D[i-1,j-1]×T_DM)
I[i,j] = M[i-1,j]×T_MI + I[i-1,j]×T_II
D[i,j] = M[i,j-1]×T_MD + D[i,j-1]×T_DD

算法实现时采用了多项优化：

1. 初始化策略：从假想的"第0位置"开始，设置为缺失状态
2. 数值稳定性处理：使用大数初始化防止下溢
3. 计算优化：对连续相同单倍型重用部分计算结果

实际应用考虑

在实际应用中，GATK团队做出了多项工程优化：

1. 默认使用固定的indel起始质量值(Phred 45)
2. 当BQSR indel质量值不可用时使用默认值
3. 算法实现针对向量化架构进行了优化

总结

Pair HMM模型为GATK的变异检测提供了精确的比对概率计算框架，其动态规划实现既保证了计算效率，又保持了模型的准确性。理解这一模型的原理对于正确使用和优化GATK工具具有重要意义。

通过这种概率化的比对方法，GATK能够更准确地处理indel等复杂变异，提高变异检测的敏感性和特异性，特别是在低复杂度区域和重复区域。

gatk Official code repository for GATK versions 4 and up 项目地址: https://gitcode.com/gh_mirrors/ga/gatk