序列比对

Pattern Matching 与 Sequence Alignment

Pattern Matching（模式匹配）与 Sequence Alignment（序列比对）前者关注“存在性与定位”，后者关注“演化距离”。

核心差异

• Pattern Matching:

  • 定义：在一个长文本（Text）中查找一个较短模式（Pattern）的精确或近似出现位置。

  • 底层逻辑：字符串搜索（String Searching）。它是刚性的，通常不预设生物学模型，只关注字符的匹配状态。

  • 典型算法：KMP, Boyer-Moore algoritm, Aho-Corasick, Suffix Tree。

  • 例子：在人类基因组中寻找限制性内切酶 EcoRI 的识别位点 GAATTC。或者在 Read Mapping 阶段，将一条 150bp 的 Read 快速定位到参考基因组上（允许少量 Mismatch）。

• Sequence Alignment:

  • 定义：将两个或多个序列按线性顺序排列，通过引入空位（Gap）来最大化残基之间的相似性得分。

  • 底层逻辑：最优化问题（Optimization Problem）。它是弹性的，旨在量化两个序列从共同祖先分化后的演化距离（Evolutionary Distance）。

  • 典型算法：Needleman-Wunsch, Smith-Waterman, BLAST (HSP extension phase)。

  • 例子：比较人与小鼠的 p53 基因，确定哪些区域是保守的功能域，哪些区域发生了插入或缺失（Indel）事件。

历史坐标

序列比对的算法基石建立在两个里程碑式的工作之上：

1. Global Alignment: 1970 年，Needleman 与 Wunsch 首次将动态规划引入生物序列分析，解决了全长序列的同源性度量问题。

2. Local Alignment: 1981 年，Smith 与 Waterman 针对生物序列中普遍存在的局部保守性（如结构域），提出了局部比对算法。

这一领域的数学核心，在于如何量化序列间的差异。虽然 Hamming Distance 定义了等长字符串之间不同字符的个数，但它无法处理 Indel。为了模拟生物演化，我们需要引入 Levenshtein Distance（编辑距离），并在赋权的有向无环图中寻找最优路径。

---

构建 Alignment Graph

拓扑结构

给定两个序列 $X=x_1...x_M$ 和 $Y=y_1...y_N$，我们将比对过程映射为一个网格图（Grid Graph）：

• 节点 (Nodes)：由坐标 $(i,j)$ 表示，对应序列前缀 $X[\mathrm{1..}i]$ 和 $Y[\mathrm{1..}j]$ 的比对状态。总节点数为 $(M+1)(N+1)$。

• 边 (Edges)：连接相邻节点，代表一步演化操作。

  • 对角边 $(i-1,j-1)\to (i,j)$：代表 $x_i$ 与 $y_j$ Match（匹配）或 Mismatch（错配/突变）。权重为 $s(x_i,y_j)$。

  • 水平边 $(i,j-1)\to (i,j)$：代表在序列 $X$ 中插入空位（或 $Y$ 中的字符被保留），即 Insertion。权重为 Gap Penalty $d$。

  • 垂直边 $(i-1,j)\to (i,j)$：代表在序列 $Y$ 中插入空位（或 $X$ 中的字符被删除），即 Deletion。权重为 Gap Penalty $d$。

NOTE

虽然这种形式的alignment graph经常被绘制为表格的形式, 但是把他理解为图对于后续理解affine gap score, Local Alignment非常有好处.

全局比对 (Needleman-Wunsch)

在此 DAG 中，Global Alignment 等价于寻找从 Source $(0,0)$ 到 Sink $(M,N)$ 的最长路径（Longest Path）。这里的“长”指路径权重之和最大。

令 $S_{i,j}$ 为从 Source 到节点 $(i,j)$ 的最高得分，状态转移方程即为 DAG 中汇入节点 $(i,j)$ 的三条边的决策过程：

$S_{i,j}=\max \{\begin{array}{l}{S}_{i-1,j-1}+s(x_i,y_j)\\ S_{i-1,j}-d\\ S_{i,j-1}-d\end{array}$

复杂度的权衡

这种基础 DP 能够保证找到数学上的最优解，时间复杂度为 $O(MN)$，空间复杂度为 $O(MN)$。对于细菌基因组尚可接受，但对于真核生物染色体级别的比对，我们需要更高级的策略。

---

Gap open/extend：Affine Gap Penalty 的状态机视角

生物学矛盾

在基础 DP 中，Gap Penalty 是线性的（Linear Gap Penalty），即 $Cost(k)=k\cdot d$。这意味着一个长度为 10 的 Gap 的代价是长度为 1 的 Gap 的 10 倍。

然而，生物学事实并非如此。在分子演化中，插入/缺失事件的发生（Opening） 是小概率事件，但一旦发生，聚合酶的滑移或转座子的插入往往会导致连续的一串核苷酸改变。因此，延长一个已存在的 Gap（Extension） 的代价应当远小于 开启一个新 Gap（Opening） 的代价。

这就是 Affine Gap Penalty 模型：

$Cost(k)=\rho +k\cdot \sigma$

其中 $\rho$ 是 Gap Open Penalty，$\sigma$ 是 Gap Extension Penalty，且 $\rho \gg \sigma$。

二维表失效了？

如果我们试图在单层 DAG 矩阵上直接应用此公式，会遇到困难：当我们计算 $S_{i,j}$ 时，若选择从 $S_{i,j-1}$（水平边）移动过来，我们无法得知前一步是刚刚开启 Gap，还是已经在延伸 Gap。也就是说，简单的二维节点没有足够的“记忆”来区分 Gap 的状态。

三层图模型 (Three-State Machine)

为了解决记忆问题，我们需要将 DAG 从二维平面扩展为三层立体结构。每一层代表一种特定的比对状态：

1. Middle Layer (M)：匹配/错配层。代表当前两个字符 $x_i,y_j$ 对齐。

2. Lower Layer (X)：垂直插入层。代表当前正在扩展 $X$ 序列的 Gap（即 $Y$ 序列有残基，$X$ 为空）。

3. Upper Layer (Y)：水平插入层。代表当前正在扩展 $Y$ 序列的 Gap。

层间跃迁逻辑（边的重定义）：

• 同层移动（Extension）：

  • 在 $X$ 层内垂直移动：$X_{i-1,j}\to X_{i,j}$，代价为 $\sigma$（延伸 Gap）。

  • 在 $Y$ 层内水平移动：$Y_{i,j-1}\to Y_{i,j}$，代价为 $\sigma$（延伸 Gap）。

• 跨层跃迁（Opening/Closing）：

  • 从 $M$ 层跳入 $X$ 层：意味着开启一个 Gap，代价为 $\rho +\sigma$。

  • 从 $X$ 层跳回 $M$ 层：意味着结束 Gap，通常代价为 0（或者整合进下一步的匹配分）。

修正后的递推公式

我们需要维护三个矩阵变量 $M_{i,j},I_x(i,j),I_y(i,j)$：

$\begin{array}{rl}{I}_x(i,j)& =\max \{\begin{array}{ll}{M}_{i-1,j}-(\rho +\sigma )& \text{(Gap Open)}\\ I_x(i-1,j)-\sigma & \text{(Gap Extend)}\end{array}\\ & \\ I_y(i,j)& =\max \{\begin{array}{ll}{M}_{i,j-1}-(\rho +\sigma )& \text{(Gap Open)}\\ I_y(i,j-1)-\sigma & \text{(Gap Extend)}\end{array}\\ & \\ M_{i,j}& =\max \{\begin{array}{l}{M}_{i-1,j-1}+s(x_i,y_j)\\ I_x(i,j)\\ I_y(i,j)\end{array}\end{array}$

通过这种状态机的形式化，我们将生物学上合理的非线性罚分，成功转化为了可以在 DAG 上计算的线性递推，仅使常数项增加了 3 倍，时间复杂度依然保持在 $O(MN)$。

---

局部比对 (Local Alignment) 的图论解释

Needleman-Wunsch 强制从 $(0,0)$ 走到 $(M,N)$，这假设了两个序列在全长上都是同源的。但在基因组学中，我们经常比较两个仅在特定外显子或结构域上同源的基因，其余部分可能是完全无关的“噪声”。

Smith-Waterman 算法对 DAG 做了极其简单却深刻的修改：允许零分重置。

$S_{i,j}=\max (S_{i-1,j-1}+s,S_{i-1,j}-d,S_{i,j-1}-d,0)$

从图论的角度看，这等价于：

1. Source 的泛化：从 $(0,0)$ 到图中的每一个节点 $(i,j)$ 都连接了一条权重为 0 的虚拟边。这意味着比对路径可以从矩阵的任何位置免费开始。

2. Sink 的泛化：比对路径可以在任意位置终止。我们不再寻找 $S_{M,N}$，而是寻找矩阵中的全局最大值 ${\max }_{i,j}(S_{i,j})$ 作为终点，并从此反向回溯直到遇到 0。

这不仅解决了局部同源性问题，也自动过滤了不相关的低分区域。

---

空间/时间优化：Hirschberg 与 Banding

在处理长序列（如两整条细菌染色体）时，标准的 DP 面临两个瓶颈：空间（内存爆炸）和时间（计算太慢）。

空间优化：Hirschberg Algorithm

标准的 DP 需要存储整个 $M\times N$ 的矩阵来进行回溯（Backtracking），对于 $10^5$ 长度的序列，这需要几十 GB 的内存。

Hirschberg 算法利用了分治策略 (Divide and Conquer) 将空间复杂度从 $O(MN)$ 降至 $O(\min (M,N))$。

• 原理：计算 $S_{i,j}$ 的数值只需要上一行的信息（滚动数组），这只需要线性空间。难点在于找路径。

• Middle Edge：最优的比对路径一定会穿过矩阵的中间列（$j=N/2$）。

• 操作：

  1. 从左上角 $(0,0)$ 做 Forward DP 到中间列，得到分数向量 $Scor{e}_L$。

  2. 从右下角 $(M,N)$ 做 Backward DP 到中间列，得到分数向量 $Scor{e}_R$。

  3. 在中间列寻找索引 $k$，使得 $Scor{e}_L[k]+Scor{e}_R[k]$ 最大。连接 $(k,N/2)$ 的边即为路径上的 Middle Edge。

  4. 从 $k$ 点将矩阵切分为左上和右下两个子矩阵，递归执行上述过程。

虽然这看起来计算了两次，但根据几何级数求和 $\sum \frac{1}{2^n}=2$，总的时间复杂度仅是原来的 2 倍，换来的是指数级的空间节省。

启发式时间优化：Banding (条带化)

即便空间被优化，对于差异极小的同源序列（例如比较人和黑猩猩的同源染色体），计算整个 $N\times N$ 矩阵的大部分区域是浪费的。因为两条高度相似的序列，其最优比对路径大概率会沿着主对角线 (Main Diagonal) 附近游走。

Banding 的核心思想是：剪枝 (Pruning)。我们假设最优路径不会偏离对角线太远，因此只计算对角线周围宽度为 $k$ 的“条带”区域。

算法逻辑

设定带宽参数 $k$（Bandwidth）。对于节点 $(i,j)$，如果 $|i-j|>k$，我们直接令 $S_{i,j}=-\infty$（或不予计算）。

优缺点分析

• 优势：时间复杂度从 $O(N^2)$ 骤降为 $O(kN)$。由于 $k\ll N$（通常 $k$ 只有几百），这使得超长序列的比对成为可能。

• 风险：如果序列中存在一个长度超过 $k$ 的 Indel，真实的比对路径就会“逃逸”出这个条带。此时 Banding 算法会撞上 $-\infty$ 的边界，导致比对中断或产生错误结果。

实际应用

在实际软件（如 BLAST 或 FASTA）的某些阶段，以及长读长测序（PacBio/Nanopore）的初始比对中，Banding 是极其常用的策略。它通常与 Seeding（种子法）结合使用：先通过 K-mer 匹配找到确定的锚点（Anchor），然后在锚点之间进行 Banded DP，以确保路径不会偏离。

---

总结

序列比对算法的演进史，模型精确度与计算资源之间不断在寻找平衡点：

1. Needleman-Wunsch 将比对形式化为图的最长路径问题。

2. Affine Gap 通过扩展图的维度（三层状态机），解决了 Indel 的生物学合理性。

3. Smith-Waterman 通过引入零分重置，解决了功能域的局部搜索问题。

4. Hirschberg 利用递归分治，打破了线性增长的空间限制。

5. Banding 利用先验知识（序列相似性），通过剪枝换取了速度的飞跃。