双序列比对：用动态规划破译基因的亲缘关系

序列比对是生物信息学的核心操作，本文用 Needleman-Wunsch 全局比对算法，带你用动态规划在 O(mn) 时间内求出两段 DNA/蛋白序列的最优对齐，并给出可运行的 Python 实现。

当你在浏览器里打开一个物种的 DNA 序列，看到一长串 A、T、C、G 时，一个自然的问题是：两段序列究竟有多相似？这种"序列比对"（sequence alignment）正是 生物信息学 最核心的操作之一——它支撑着基因注释、物种进化树构建，乃至 BLAST 这样的海量数据库搜索。

从"穷举"到"动态规划"

如果把两条长度分别为 m、n 的序列做比对，所有可能的对齐方式是指数级的，根本无法穷举。1970 年，Needleman 与 Wunsch 提出了用 动态规划（dynamic programming）在 $O(mn)$ 时间内求得 全局最优 比对的算法。其核心是一个得分矩阵 $F$ ，其中 $F(i,j)$ 表示序列 X 前 i 个字符与序列 Y 前 j 个字符的最优比对得分。

图1：得分矩阵 $F$ 的递推关系——每个格子由左、上、左上三个邻居取最大值得到

状态转移方程为：

F(i,j) = \\max\\begin{cases} F(i-1,j-1) + s(x_i, y_j) \\\\ F(i-1,j) + \\text{gap} \\\\ F(i,j-1) + \\text{gap} \\end{cases}

其中 $s(x_i, y_j)$ 是 替换矩阵（如 PAM 或 BLOSUM）给出的匹配/错配得分，gap 是空位罚分。矩阵填完后，从 $F(m,n)$ 沿取值来源反向回溯（traceback），即可还原出具体的对齐方式。

用 Python 实现

下面是一段最小可用的全局比对实现：

def needleman_wunsch(x, y, match=1, mismatch=-1, gap=-2):
    m, n = len(x), len(y)
    F = [[0] * (n + 1) for _ in range(m + 1)]
    # 初始化首行首列
    for i in range(m + 1):
        F[i][0] = i * gap
    for j in range(n + 1):
        F[0][j] = j * gap
    # 填表
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            diag = F[i-1][j-1] + (match if x[i-1] == y[j-1] else mismatch)
            up = F[i-1][j] + gap
            left = F[i][j-1] + gap
            F[i][j] = max(diag, up, left)
    return F[m][n]

print(needleman_wunsch("GATTACA", "GCATGCU"))  # 输出比对得分

图2：对 GATTACA 与 GCATGCU 的全局比对结果，连字符表示空位

从全局到局部

Needleman-Wunsch 求的是"覆盖整条序列"的全局比对；而 Smith-Waterman 算法则通过允许得分为负时归零，专找"局部最相似片段"——这正对应 BLAST 的底层思想。理解了这两个 $O(mn)$ 的表格算法，你就握住了现代 序列分析 的钥匙。