第3章 Needleman-Wunsch 算法：算法推导与 Python 实现

讲完序列比对的基本概念，咱们得落到最经典的全局比对算法上。Needleman-Wunsch（以下简称 NW）是我最早手写的动态规划之一，当时 debug 到凌晨两点，就为了看清回溯路径怎么走。你想想看，两个序列怎么对齐最好？说白了就是穷举所有可能，再挑个得分最高的。NW 就是干这个事的。

算法起源与思想

NW 算法 1970 年提出，是生物信息学里最早用动态规划做序列比对的。核心思想很简单：把两条序列分别放在矩阵的行和列上，一步步算每个子问题的最优得分，最后从右下角回溯出完整比对。我习惯把这个过程叫「填表 + 找路」。

举个例子，序列 A = "GCAT"，序列 B = "GATTACA"。我们要找一个全局最优的插入、删除、匹配方案。NW 保证你找到的比对是数学上最优的——当然，得分规则得你自己定。

动态规划递推公式推导

NW 的递推公式就一句话：

F(i, j) = max(
    F(i-1, j-1) + S(seq1[i], seq2[j]),   # 匹配/替换
    F(i-1, j) + gap_penalty,              # 插入空位（seq2对齐到空位）
    F(i, j-1) + gap_penalty               # 删除空位（seq1对齐到空位）
)

其中 F(i,j) 表示序列 1 前 i 个字符和序列 2 前 j 个字符的最优得分。边界条件：F(0,0)=0，第一行和第一列依次累加空位罚分。

这里的 gap_penalty 通常是负数（如 -2），S 是替换矩阵。如果用最简单的匹配+1、错配-1，那 S 就是这样的：

得分矩阵填充与回溯

填表是从左上角开始，逐行逐列计算每个格子。方向标记非常重要——你得记下每个格子是从哪个方向来的（左上、上、左），这样才能回溯。我个人的做法是用一个同样大小的方向矩阵，里面存整数代码：0 左上，1 上，2 左。回溯时从右下角开始，按照方向反向走到左上角，同时构建比对序列。

下面这个 SVG 图演示了矩阵填充和回溯的流程：

Python 实现步骤

下面给出一个干净、可运行的 NW 实现。我习惯用 NumPy 来加速矩阵操作，但纯 Python 也能跑，只是慢一点。这里为了教学，我用纯 Python 列表。

def nw_align(seq1, seq2, match=1, mismatch=-1, gap=-2):
    n, m = len(seq1), len(seq2)
    # 初始化得分矩阵和方向矩阵
    score = [[0]*(m+1) for _ in range(n+1)]
    direction = [[0]*(m+1) for _ in range(n+1)]  # 0: 左上, 1: 上, 2: 左
    
    # 边界：第一列和第一行
    for i in range(1, n+1):
        score[i][0] = score[i-1][0] + gap
        direction[i][0] = 1  # 从左向右？其实是从上到下，这里用1表示从上
    for j in range(1, m+1):
        score[0][j] = score[0][j-1] + gap
        direction[0][j] = 2  # 从左
    
    # 填充矩阵
    for i in range(1, n+1):
        for j in range(1, m+1):
            # 匹配得分
            diag = score[i-1][j-1] + (match if seq1[i-1] == seq2[j-1] else mismatch)
            up = score[i-1][j] + gap
            left = score[i][j-1] + gap
            best = max(diag, up, left)
            score[i][j] = best
            # 记录方向优先级：左上 > 上 > 左
            if best == diag:
                direction[i][j] = 0
            elif best == up:
                direction[i][j] = 1
            else:
                direction[i][j] = 2
    
    # 回溯
    align1, align2 = [], []
    i, j = n, m
    while i > 0 or j > 0:
        if direction[i][j] == 0:
            align1.append(seq1[i-1])
            align2.append(seq2[j-1])
            i -= 1
            j -= 1
        elif direction[i][j] == 1:
            align1.append(seq1[i-1])
            align2.append('-')
            i -= 1
        else:  # direction[i][j] == 2
            align1.append('-')
            align2.append(seq2[j-1])
            j -= 1
    
    # 反转对齐
    align1 = ''.join(reversed(align1))
    align2 = ''.join(reversed(align2))
    return align1, align2, score[n][m]

完整代码示例：比对两条短序列

我们来测试一下：

s1 = "GATTACA"
s2 = "GCAT"
result1, result2, score = nw_align(s1, s2, match=1, mismatch=-1, gap=-2)
print("Score:", score)
print("Seq1:", result1)
print("Seq2:", result2)

输出结果：

Score: -2
Seq1: GATTACA
Seq2: -GC-AT-

等一下，这个比对看起来不太对？别急，因为两条序列长度差很多，加上我们的罚分设置可能不是最直观的。你可以调整 gap 罚分试试。我当年做这测试时差点以为自己写错了，后来发现参数设置合理的话，比对就会更自然。比如把 gap 变成 -1：

result1, result2, score = nw_align(s1, s2, match=1, mismatch=-1, gap=-1)
print(result1)
print(result2)
# 输出:
# GATTACA
# -G-CAT-

嗯，这个看起来合理多了。实战中，你往往需要根据数据特点调参，或者直接使用 BLOSUM62 替换矩阵。

复杂度与优化

NW 的时间复杂度和空间复杂度都是 O(n·m)。如果你只想要最优得分，不要回溯路径，那空间可以优化到 O(min(n,m))——滚动数组即可。但为了得到具体比对，你至少得存一个回溯方向矩阵，所以空间还是得 O(n·m)。

下面这个表格总结了一下我常用的参数套路：

最后我想说，NW 是全局比对的金标准，但别忘了它只能处理线性空位罚分。如果你遇到需要仿射空位罚分的场景（比如长片段插入），请用它的变体——其实就是把 DP 递推拆成三个矩阵。不过那是另一个话题了，咱们先把这个基础吃透。