Profile HMM

快速概览

Profile HMM 是 PWM/PSSM 的强大扩展：它通过 Match、Insert、Delete 三类状态，能够显式建模序列家族中的插入和缺失，适用于蛋白家族和结构域分析。

核心是 M/I/D 三类状态组成的线性结构
用 Viterbi 算法找最优比对路径，用 Forward 算法评估匹配程度
它是蛋白家族数据库（如 Pfam）的核心建模方法

所属板块 核心方法

索引、比对、组装与概率模型构成的核心算法主轴。

阅读目标 帮助建立阅读上下文

先判断这页与你当前问题的关系，再决定是否深入展开。

建议前置 先建立相关基础对象与方法直觉

建议先建立相关基础对象与方法直觉，再进入本页。

是什么

PWM/PSSM 很适合描述一个固定长度的 motif：

每个位置都有字符偏好；
候选片段通常与 motif 长度相同；
不擅长显式表示插入和缺失。

而 Profile HMM 则进一步解决：

如果一个序列家族在不同成员之间存在插入、缺失和长度变化，如何仍然用概率模型刻画它的”家族结构”？

数学模型

状态定义

Profile HMM 由三类状态组成：

Match 状态（M_i）: 表示家族中的保守位置 i，发射一个字符。
Insert 状态（I_i）: 表示在位置 i 之后插入额外字符，可以循环发射多个字符。
Delete 状态（D_i）: 表示跳过保守位置 i，不发射字符。

状态转移结构

对于长度为 L 的 motif，Profile HMM 的状态序列为：

M_0 → M_1 → M_2 → ... → M_L → M_\{L+1\}
  ↓     ↓     ↓           ↓     ↓
  I_0   I_1   I_2         I_L   I_\{L+1\}
        ↓     ↓           ↓
        D_2   D_3         D_\{L+1\}

其中：

M_0 和 M_{L+1} 是开始和结束状态
每个 M_i 可以转移到 M_{i+1}、I_i 或 D_{i+1}
每个 I_i 可以转移到 M_{i+1} 或 I_i（自身循环）
每个 D_i 可以转移到 M_{i+1} 或 D_{i+1}

概率参数

转移概率 $a_{uv}$ ：从状态 u 转移到状态 v 的概率
发射概率 $e_i(c)$ ：Match 或 Insert 状态 i 发射字符 c 的概率

Viterbi 算法与特殊状态处理

在 Profile HMM 中，Viterbi 算法需要处理三种状态的转移。特别需要注意的是 Delete 状态（D），因为它不发射任何字符（称为”静默状态”）。

定义 $v^M_j(i)$ 为观测序列前缀 $x_1 \dots x_i$ 与模型前缀（直到状态 $M_j$ ）对齐的最优对数似然得分。

Match 状态（ $M_j$ ）： $v^M_j(i) = \log \frac{e_{M_j}(x_i)}{p(x_i)} + \max \begin{cases} v^M_{j-1}(i-1) + \log a_{M_{j-1}M_j} \\ v^I_{j-1}(i-1) + \log a_{I_{j-1}M_j} \\ v^D_{j-1}(i-1) + \log a_{D_{j-1}M_j} \end{cases}$
Insert 状态（ $I_j$ ）： $v^I_j(i) = \log \frac{e_{I_j}(x_i)}{p(x_i)} + \max \begin{cases} v^M_j(i-1) + \log a_{M_jI_j} \\ v^I_j(i-1) + \log a_{I_jI_j} \\ v^D_j(i-1) + \log a_{D_jI_j} \end{cases}$
Delete 状态（ $D_j$ ）：由于不消耗字符，得分直接从上一个位置转移： $v^D_j(i) = \max \begin{cases} v^M_{j-1}(i) + \log a_{M_{j-1}D_j} \\ v^I_{j-1}(i) + \log a_{I_{j-1}D_j} \\ v^D_{j-1}(i) + \log a_{D_{j-1}D_j} \end{cases}$

注意：Delete 状态的计算不需要 $i-1$ ，这在动态规划填表时需要特别注意计算顺序，确保在计算 $v^D_j(i)$ 时， $v^M_{j-1}(i)$ 等同一列的状态已经更新。

为什么 PWM 不够用

对于简单 DNA motif，PWM 已经足够强大。

但对于蛋白家族或更复杂的序列模式，常见问题是：

某些位置可以缺失；
某些位置允许插入额外残基；
不同成员的长度并不完全一致；
只靠固定长度窗口会把这些变化都当成”糟糕匹配”。

这时，需要一个能显式表示”匹配 / 插入 / 缺失”的模型。

算法步骤

构建 Profile HMM

算法1：从多序列比对构建 Profile HMM

输入：多序列比对 MSA，motif 长度 L
输出：Profile HMM 参数

1. 初始化转移和发射计数矩阵
2. for each 序列 in MSA:
     确定状态路径（M/I/D）
     更新转移计数
     更新发射计数（M 和 I 状态）
3. 添加伪计数避免零概率
4. 归一化得到概率
5. return 转移概率 a_ij 和发射概率 e_i(c)

时间复杂度：O(N · L)，N 是序列数量

Viterbi 对齐

算法2：使用 Viterbi 算法寻找最优状态路径

输入：序列 s，Profile HMM 参数
输出：最优状态路径和得分

1. 初始化：V_0(M_0) = 1, 其他状态 = -∞
2. for t = 1 to |s|:
     for each 状态 j:
        if j 是 M 或 I 状态:
           V_t(j) = max_i [V_\{t-1\}(i) · a_ij · e_j(s_t)]
           记录回溯指针
        else if j 是 D 状态:
           V_t(j) = max_i [V_t(i) · a_ij]  // 不发射字符
           记录回溯指针
3. 回溯得到最优状态路径
4. return 路径和得分

时间复杂度：O(|s| · L · K)，K 是状态数（通常为 3L）

示例

示例：简单蛋白家族

假设我们有一个简单的蛋白家族，保守骨架为 3 个位置：

序列 1: A G T
序列 2: A Q G R T  (在 G 和 T 之间插入了 R)
序列 3: A T         (跳过了 G)

步骤 1：构建状态结构

M_0 → M_1 → M_2 → M_3 → M_4
  ↓     ↓     ↓     ↓     ↓
  I_0   I_1   I_2   I_3   I_4
        ↓     ↓     ↓
        D_2   D_3   D_4

步骤 2：估计参数

从对齐中统计：

发射概率（简化示例）：
- M_1: A=1.0 (非常保守)
- M_2: G=0.67, T=0.33
- M_3: T=1.0
- I_2: R=1.0 (插入位置偏好 R)
转移概率（简化）：
- M_1 → M_2: 0.67, M_1 → D_2: 0.33
- M_2 → M_3: 0.67, M_2 → I_2: 0.33
- I_2 → M_3: 1.0

步骤 3：对新序列评分

序列：A G R T

可能路径 1：M_1 → M_2 → I_2 → M_3

得分 = P(M_1→M_2) · P(A|M_1) · P(M_2→I_2) · P(G|M_2) · P(I_2→M_3) · P(R|I_2) · P(T|M_3)
= 0.67 · 1.0 · 0.33 · 0.67 · 1.0 · 1.0 · 1.0 = 0.148

可能路径 2：M_1 → M_2 → M_3

得分 = 0.67 · 1.0 · 0.67 · 0.67 · 1.0 = 0.301

路径 2 得分更高，说明 A G R T 更可能通过匹配 G 和 T 来对齐，而 R 被视为错配。

复杂度分析

时间复杂度

构建 Profile HMM：O(N · L)，N 是序列数量，L 是 motif 长度
Viterbi 对齐：O(|s| · L · K)，其中 K ≈ 3 是每个位置的状态数
Forward/Backward 评估：O(|s| · L · K)

空间复杂度

存储参数：O(L · K · |Σ|)，|Σ| 是字母表大小（如 20 个氨基酸）
Viterbi 矩阵：O(|s| · L · K)
可优化到：O(L · K)（只保留当前行）

Profile HMM 的核心结构

Profile HMM 通常由三类状态构成：

Match states (M)：表示家族中的保守列；
Insert states (I)：表示相对于保守骨架插入的字符；
Delete states (D)：表示某个保守位置被跳过。

因此，一条序列与 profile 的关系不再只是：

“每个位点得多少分”；

而是：

“它沿着这条 match/insert/delete 状态路径走了一遍”。

一个直观例子

假设一个蛋白家族的保守骨架大致是：

A - G - T

其中第二和第四位置之间可能允许插入一些残基，那么 Profile HMM 可以表示为：

M1：偏好 A
M2：偏好 G
M3：偏好 T
在相邻 match 状态之间插入 I 状态
对某些 match 位置允许通过 D 状态跳过

这样，序列：

AGT
AQGRT
AT

都可以被统一地解释为”属于这个家族”，只是经过的状态路径不同。

与 Viterbi / Forward 的关系

Profile HMM 的推断方式和一般 HMM 一样，也依赖：

Viterbi：找最可能的 match/insert/delete 路径；
Forward/Backward：评估序列和 profile 的总体兼容程度；

因此，Profile HMM 可以被看作：

PWM/PSSM 的更强表达形式；
又是一般 HMM 思想在序列家族问题中的专门化版本。

与 gene prediction 的区别

Profile HMM：
- 强调一个序列是否符合某个家族/域的 profile；
- 常用于蛋白家族比对、结构域搜索；
Gene prediction HMM：
- 强调把一段长基因组序列切分成 exon/intron/intergenic 等状态；
- 更关注整条基因结构的解释。

两者都属于 HMM，但应用目标不同。

与真实工具或流程的连接

Profile HMM 在蛋白家族分析中非常重要，典型用途包括：

family/domain search；
远缘同源检测；
蛋白功能注释；
结构域边界判断。

它是从”简单 motif 打分”走向”序列家族概率建模”的关键一步。