Needleman-Wunsch 算法
Needleman-Wunsch 算法是全局比对的标准动态规划算法：它通过系统地枚举所有可能的比对路径，找到两条序列之间得分最高的全局比对。
- 掌握递推公式和回溯算法是理解的关键
- 它是所有全局比对算法的起点，后续优化都基于此框架
Needleman-Wunsch 算法是 1970 年提出的第一个用于序列比对的动态规划算法。它解决了以下问题：
给定两条序列和一个打分体系，找到使总得分最大的全局比对。
这里的"全局"意味着：两条序列必须从头到尾完全对齐，序列两端也会被解释为匹配、替换或 gap。
要解决什么生物信息学问题
Section titled “要解决什么生物信息学问题”在实际生物信息学应用中，Needleman-Wunsch 适合的场景包括：
- 比较两个长度相近的同源基因或蛋白序列
- 分析两个物种中直系同源基因的整体差异
- 作为教学示例理解动态规划在序列比对中的应用
- 为更复杂的算法提供基础框架
Needleman-Wunsch 需要定义三个参数：
- 匹配得分
s(a, a)：相同字符对齐时的奖励 - 错配得分
s(a, b)：不同字符对齐时的惩罚（通常为负值） - Gap 罚分
g：插入或删除一个字符的代价（通常为负值）
动态规划矩阵定义
Section titled “动态规划矩阵定义”设序列 X = x₁x₂...x_m，序列 Y = y₁y₂...y_n。
定义 F[i][j] 为：序列 X 的前 i 个字符与序列 Y 的前 j 个字符的最优比对得分。
边界条件的含义：
- 两个空串对齐得分为 0
- 前
i个字符与空串对齐，需要插入i个 gap - 空串与前
j个字符对齐，需要删除j个字符
核心递推公式
Section titled “核心递推公式”当 i > 0 且 j > 0 时：
递推直觉图解
Section titled “递推直觉图解”当前单元格 F[i][j] 的值来自三个方向：对角线（F[i-1][j-1] + s(xi, yj)，匹配/替换）、上方（F[i-1][j] + g，X 序列 Gap）和左方（F[i][j-1] + g，Y 序列 Gap）。
这个递推体现了动态规划的核心思想：当前位置的最优解来自三个更小子问题中的最优选择。
步骤 1：初始化矩阵
Section titled “步骤 1：初始化矩阵”- 创建一个
(m+1) × (n+1)的矩阵F - 填充边界条件
步骤 2：填充矩阵
Section titled “步骤 2：填充矩阵”按照行优先或列优先的顺序，从 F[1][1] 开始，使用递推公式填充整个矩阵。
步骤 3：回溯得到比对
Section titled “步骤 3：回溯得到比对”从 F[m][n] 开始，回溯到 F[0][0]。在每一步，检查当前值来自哪个方向（对角线 → 匹配/替换，上方 → X 序列 Gap，左方 → Y 序列 Gap），然后沿该方向移动。
- 如果
F[i][j]来自F[i-1][j-1] + s(x_i, y_j)：- 将
x_i和y_j对齐 - 移动到
F[i-1][j-1]
- 将
- 如果
F[i][j]来自F[i-1][j] + g：- 将
x_i与 gap 对齐 - 移动到
F[i-1][j]
- 将
- 如果
F[i][j]来自F[i][j-1] + g：- 将 gap 与
y_j对齐 - 移动到
F[i][j-1]
- 将 gap 与
回溯完成后，将路径反转得到最终比对。
考虑两条 DNA 序列：
X = GATTACA（长度 m = 7）Y = GCATGCU（长度 n = 7）
使用以下打分体系：
- 匹配得分：
+1 - 错配得分：
-1 - Gap 罚分：
-2
步骤 1：初始化
Section titled “步骤 1：初始化”| | ε | G | C | A | T | G | C | U | |---|---|---|---|---|---|---|---|---| | ε | 0 | -2 | -4 | -6 | -8 | -10 | -12 | -14 | | G | -2 | | | | | | | | | A | -4 | | | | | | | | | T | -6 | | | | | | | | | T | -8 | | | | | | | | | A | -10 | | | | | | | | | C | -12 | | | | | | | | | A | -14 | | | | | | | |
步骤 2：填充矩阵
Section titled “步骤 2：填充矩阵”以 F[1][1] 为例（G 对 G）：
- 匹配：
F[0][0] + s(G,G) = 0 + 1 = 1 - X 插入 gap：
F[0][1] + g = -2 + (-2) = -4 - Y 插入 gap：
F[1][0] + g = -2 + (-2) = -4 - 取最大值：
F[1][1] = 1
再以 F[3][4] 为例（T 对 T）：
- 匹配：
F[2][3] + s(T,T) = 0 + 1 = 1 - X 插入 gap：
F[2][4] + g = -2 + (-2) = -4 - Y 插入 gap：
F[3][3] + g = -1 + (-2) = -3 - 取最大值：
F[3][4] = 1
完整矩阵如下：
| | ε | G | C | A | T | G | C | U | |---|---|---|---|---|---|---|---|---| | ε | 0 | -2 | -4 | -6 | -8 | -10 | -12 | -14 | | G | -2 | 1 | -1 | -3 | -5 | -7 | -9 | -11 | | A | -4 | -1 | 0 | 0 | -2 | -4 | -6 | -8 | | T | -6 | -3 | -2 | -1 | 1 | -1 | -3 | -5 | | T | -8 | -5 | -4 | -3 | 0 | 0 | -2 | -4 | | A | -10 | -7 | -6 | -3 | -2 | -1 | -1 | -3 | | C | -12 | -9 | -6 | -5 | -4 | -3 | 0 | -2 | | A | -14 | -11 | -8 | -5 | -6 | -5 | -2 | -1 |
最终得分：F[7][7] = -1
步骤 3：回溯
Section titled “步骤 3：回溯”从 F[7][7] = -1 开始，逐步检查每个单元格的值来自哪个方向：
| 当前 | 来源 | 计算 | 操作 |
|------|------|------|------|
| F[7][7] = -1 | F[6][6] + s(A,U) | 0 + (-1) = -1 | A 对 U（错配） |
| F[6][6] = 0 | F[5][5] + s(C,C) | -1 + 1 = 0 | C 对 C（匹配） |
| F[5][5] = -1 | F[4][4] + s(A,G) | 0 + (-1) = -1 | A 对 G（错配） |
| F[4][4] = 0 | F[3][3] + s(T,T) | -1 + 1 = 0 | T 对 T（匹配） |
| F[3][3] = -1 | F[2][2] + s(T,A) | 0 + (-1) = -1 | T 对 A（错配） |
| F[2][2] = 0 | F[1][1] + s(A,C) | 1 + (-1) = 0 | A 对 C（错配） |
| F[1][1] = 1 | F[0][0] + s(G,G) | 0 + 1 = 1 | G 对 G（匹配） |
回溯路径全部沿对角线，说明本例的最优全局比对不需要引入 gap：
G A T T A C A| × × | × | ×G C A T G C U验证得分：+1 -1 -1 +1 -1 +1 -1 = -1 ✓
这个结果有一个教学要点：当两条序列长度相同且整体相似度尚可时，最优全局比对可能恰好是逐位对齐，不引入任何 gap。gap 的引入需要"错开"后续所有位置，在 gap 罚分较重（本例 g = -2）时代价很高。
- 填充矩阵：
O(mn)，每个单元格需要常数时间计算 - 回溯：
O(m + n)，最多走m + n步 - 总时间复杂度：
O(mn)
- 存储完整矩阵：
O(mn) - 如果只需要得分不需要回溯：可以优化到
O(min(m, n))（只保留两行或两列）
与真实工具或流程的连接
Section titled “与真实工具或流程的连接”Needleman-Wunsch 算法本身很少直接用于大规模数据分析，主要因为：
- 时间和空间复杂度对长序列过高
- 现代工具通常使用更复杂的打分模型（如 affine gap penalty）
但它在以下方面仍然重要：
- 作为教学示例理解动态规划在序列比对中的应用
- 作为更复杂算法（如带仿射 gap 罚分、多序列比对）的基础框架
- 用于小规模序列的精确比对
现代工具如 BLAST、BWA 等虽不直接使用 Needleman-Wunsch，但其扩展阶段的精确比对仍然基于类似的动态规划原理。
算法优化与变体
Section titled “算法优化与变体”Hirschberg 算法
Section titled “Hirschberg 算法”将空间复杂度从 O(mn) 优化到 O(min(m, n))，同时仍能得到最优比对。详见 Hirschberg 算法。
Banded Dynamic Programming
Section titled “Banded Dynamic Programming”如果已知两条序列的相似度较高，可以限制搜索带宽，将时间复杂度从 O(mn) 降低到 O(b·min(m,n))，其中 b 是带宽。详见 带状动态规划。
- Needleman, S. B., & Wunsch, C. D. (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of molecular biology, 48(3), 443-453.
- Durbin, R., Eddy, S., Krogh, A., & Mitchison, G. (1998). Biological sequence analysis: probabilistic models of proteins and nucleic acids. Cambridge university press.