Skip to main content

2.2.1 PairHMM model

这篇论文中的 PairHMM（Pairwise Hidden Markov Model）数学模型，是用来计算一条测序读段（Read）由某个候选单倍型（Haplotype）产生的似然概率。它本质上是一个概率生成模型，通过隐马尔可夫链的动态规划来计算。

Haplotype 不是指单个碱基，而是指一段连续的 DNA 序列（一串碱基）。

Haplotype（单倍型）：一段连续的 DNA 序列，代表同一条染色体上特定区域的一组变异组合。例如：A C G T A 这五个碱基构成一个 Haplotype。

一、核心变量定义​

在理解公式前，需要明确各个符号的含义：

  • ii：Read 序列上的位置（1≤i≤Nr1 \le i \le N_r）。
  • jj：Haplotype 序列上的位置（1≤j≤Nn1 \le j \le N_n）。
  • fM(i,j)f^M(i,j)：在位置 (i,j)(i,j) 处于 Match（匹配）状态的概率（或得分）。
  • fI(i,j)f^I(i,j)：在位置 (i,j)(i,j) 处于 Insertion（插入）状态的概率。
  • fD(i,j)f^D(i,j)：在位置 (i,j)(i,j) 处于 Deletion（缺失）状态的概率。
  • ρ\rho：发射概率（Emission Probability）。通常表示为 ρ(X(i),Y(j))\rho(X(i), Y(j))，即给定 Haplotype 碱基 Y(j)Y(j)，测序仪观测到 Read 碱基 X(i)X(i) 的概率（通常基于测序质量值计算）。
  • αXY\alpha_{XY}：状态转移概率（Transition Probability）。下标 X→YX \to Y 表示从状态 X 转移到状态 Y 的概率。例如 αmm\alpha_{mm} 是从 M 到 M 的转移概率。

二、公式逐行解释​

1. Match 状态：fM(i,j)f^M(i,j)​

fM(i,j)=ρ⋅(αmmfM(i−1,j−1)+αimfI(i−1,j−1)+αdmfD(i−1,j−1))f^M(i,j) = \rho \cdot (\alpha_{mm}f^M(i-1,j-1) + \alpha_{im}f^I(i-1,j-1) + \alpha_{dm}f^D(i-1,j-1))
  • 物理意义：Read 的碱基 X(i)X(i) 和 Haplotype 的碱基 Y(j)Y(j) 在此处“对齐”（匹配或错配）。
  • 位置依赖：因为消耗了 Read 和 Haplotype 各一个碱基，所以前一个状态必须是 (i−1,j−1)(i-1, j-1)。
  • 来源状态：前一个位置可能是 Match（αmm\alpha_{mm}）、Insertion（αim\alpha_{im}）或 Deletion（αdm\alpha_{dm}）。这三种情况都会转移到当前的 Match 状态。
  • 发射概率：进入 Match 状态后，会“发射”一对碱基 X(i)X(i) 和 Y(j)Y(j)，因此乘以 ρ\rho。
  • 直观理解：当前匹配的概率 = （前一步所有可能状态的概率 × 对应的转移概率）之和 × 当前碱基的观测概率。

2. Insertion 状态：fI(i,j)f^I(i,j)​

fI(i,j)=αmifM(i−1,j)+αiifI(i−1,j)f^I(i,j) = \alpha_{mi}f^M(i-1,j) + \alpha_{ii}f^I(i-1,j)
  • 物理意义：Read 上多出了一个碱基 X(i)X(i)，而 Haplotype 在此处没有对应的碱基（即 Read 相对于 Haplotype 发生了插入）。
  • 位置依赖：因为只消耗了 Read 的一个碱基，Haplotype 位置不变，所以前一个状态是 (i−1,j)(i-1, j)。
  • 来源状态：只能从 Match（αmi\alpha_{mi}）或 Insertion（αii\alpha_{ii}）转移过来。通常不允许直接从 Deletion 转移到 Insertion（即 αdi=0\alpha_{di} = 0），因为那在生物学上不自然。
  • 无 ρ\rho 项：在这个简化模型中，Insertion 状态没有显式写出 ρ\rho。实际实现中（如 GATK），插入状态可能使用一个固定的发射概率，或者将其合并到了转移概率中。
  • 直观理解：当前插入的概率 = （前一步 M 或 I 状态的概率 × 转移概率）之和。

3. Deletion 状态：fD(i,j)f^D(i,j)​

fD(i,j)=αmdfM(i,j−1)+αddfD(i,j−1)f^D(i,j) = \alpha_{md}f^M(i,j-1) + \alpha_{dd}f^D(i,j-1)
  • 物理意义：Haplotype 上多出了一个碱基 Y(j)Y(j)，而 Read 在此处没有对应的碱基（即 Read 相对于 Haplotype 发生了缺失）。
  • 位置依赖：只消耗了 Haplotype 的一个碱基，Read 位置不变，所以前一个状态是 (i,j−1)(i, j-1)。
  • 来源状态：只能从 Match（αmd\alpha_{md}）或 Deletion（αdd\alpha_{dd}）转移过来。
  • 直观理解：当前缺失的概率 = （前一步 M 或 D 状态的概率 × 转移概率）之和。

4. 最终结果​

Result=fM(Nr,Nn)+fD(Nr,Nn)Result = f^M(N_r, N_n) + f^D(N_r, N_n)
  • 物理意义：整条 Read 和整条 Haplotype 比对完成后的总似然概率。
  • 为什么没有 fIf^I：因为 Read 和 Haplotype 的末端必须对齐。如果最后处于 Insertion 状态，意味着 Read 末尾有碱基没有比对到 Haplotype 上，这通常是不允许的（或者需要额外处理）。因此，最终状态只能是 Match 或 Deletion（表示 Haplotype 末尾可能有未比对的碱基）。

三、与 BSW 数学模型的对比​

理解 PairHMM 时，可以将其与之前讨论的 BSW 进行对比：

特性BSWPairHMM
本质相似性得分最大化概率似然计算
核心运算加法与取最大值（max⁡\max）乘法与加法（概率相乘相加）
状态H, E, FM, I, D
对角线H(i−1,j−1)+S(X(i),Y(j))H(i-1,j-1) + S(X(i),Y(j))αmmfM(i−1,j−1)\alpha_{mm} f^M(i-1,j-1) 等
上方E(i,j)E(i,j)（gap 罚分）fI(i,j)f^I(i,j)（插入概率）
左方F(i,j)F(i,j)（gap 罚分）fD(i,j)f^D(i,j)（缺失概率）
初始值01（或特定的边界概率）
精度8/16-bit 整数浮点数（或对数域定点数）

在 PairHMM 中，由于概率值通常非常小（大量小于 1 的数相乘），直接计算会导致浮点数下溢。因此，实际实现（包括 GenDP 论文中提到的优化）通常会将公式转换到对数域，将乘法变成加法，即：

log⁡fM(i,j)=log⁡ρ+log⁡(αmmelog⁡fM(i−1,j−1)+… )\log f^M(i,j) = \log \rho + \log(\alpha_{mm} e^{\log f^M(i-1,j-1)} + \dots)

这涉及到 log_sum 运算（即论文 Table 4 中的 Log_sum LUT），GenDP 为此设计了专门的查找表指令来加速。

四、GenDP 论文中的上下文​

在论文 2.3 节中，作者提到：

The likelihood score is computed by the formula shown in Figure 2b, where fM,fIf^M, f^I and fDf^D represent match, insertion and deletion probabilities... The computation in PairHMM uses floating-point arithmetic.

这段描述正是对应你图片中的公式。论文指出 PairHMM 在变异检测阶段耗时占比高达 70%，是极其耗时的核心算法。

为了加速，GenDP 采用了以下策略：

  1. 剪枝（Pruning）：在 DP 表中提前终止概率极低的路径，只计算有希望的单元格（论文评估时，scan 阶段占 97.7%）。
  2. 对数域与定点数近似：用对数运算代替浮点乘法，用定点数表示对数概率，降低硬件复杂度。
  3. 波前并行：与 BSW 类似，PairHMM 也是 2D DP 表，依赖前两个波前，GenDP 的 systolic array 可以高效流水化计算。

一句话总结：PairHMM 的数学模型是一个基于 M/I/D 三状态隐马尔可夫链的二维动态规划，通过概率转移和发射来评估 Read 与 Haplotype 的匹配程度，是变异检测中最核心也最耗时的计算环节。