企业官网建设流程全解析

序列联配（Sequence Alignment）是生物信息学里几乎绕不开的入门问题，动态规划则是解决这类问题最经典的框架。这篇先聊最基本的一类：全局最优对齐——给定两条序列，通过引入空位让它们尽可能对齐，同时让所有配对位置的得分总和最高。很多人第一次看算法书，会被一堆 dp[i][j] 公式劝退，其实整个思路可以用一句话概括：把一个大的对齐问题拆成无数个小前缀的对齐问题，每个小问题只关心三种选择。我当初是手动把一个 5×4 的矩阵算完整，才真正看懂状态转移和回溯之间的关系，所以这篇文章也会从最直观的例子推演到底，再给一段能直接改着用的 Python 代码。适合刚入门算法或生信分析的同学，也适合写过代码但一直对回溯逻辑模棱两可的人。

1. 从一个最朴素的问题出发：两条序列到底怎么比

1.1 序列联配在比什么

先别急着上公式，我们看两条短序列 ACGT 和 AGT。肉眼扫一遍，谁都能看出它们大概“很像”，但真要问一句：它们哪里相同、哪里不同？你是说不清具体坐标的。序列联配要做的，就是把两条序列摆成上下两排，相同的区域对齐到同一列，用空位填掉那些“多出来”的字符。比如：

A C G T A - G T

这是一个可行的对齐。另一个可行方案是：

A C G T A G T -

两个方案都能接受，但哪一个更好？这里就需要引入一个量化标准——打分规则。

每个配对位置只有三种状态：两个真实字符且相同，叫匹配；两个真实字符但不同，叫错配；一边是空位，叫插入/缺失（统称 indel）。我们会给匹配一个正分，给错配和空位一个负分，一条对齐方式的总分就是所有列的和。“最优对齐”本质上是在所有合法排列方式里，找到总分最大的一种。

1.2 空位为什么要存在

为什么不能老老实实逐字符比对，非要引入空位这种操作？因为真实数据里，一条序列相对于另一条，经常会有插入或缺失。DNA 在进化过程中会发生碱基插入/缺失，测序过程也可能多测或少测一个碱基；哪怕是文本比较，版本之间也会增删字符。如果强行不允许空位，ACGT 和 AGT 从第一个不同字符开始就全是错配，两条明明高度相似的序列得分会非常难看。

空位本质上是一种“对齐代价的缓冲”，让你能区分“这一处真的不一样”和“这一处只是多了一段”。理解这一点，你就知道空位罚分为什么必须存在，也知道它不能被设成 0——设成 0 的话，随便乱插空位就能得到高分，比对就失去意义了。

1.3 暴力枚举为什么走不通

既然定义清楚了，最笨的办法就是把所有可能的空位排列都试一遍，挑得分最高的。问题是数量爆炸。两条长度分别为 n 和 m 的序列，对齐后的总长度可以是 0 到 n+m 之间的任意值，把空位分配到两条序列里，方式数大约是组合级别的规模，整体复杂度可以写成指数级。

哪怕只看 30 对 30 的短序列，2 的 30 次方就已经超过十亿，真实场景下 DNA 片段动辄上千个碱基，暴力枚举在计算机上根本等不到结果。这也是为什么需要用动态规划：它不枚举全部方案，而是通过保存已算好的小前缀结果，避免重复计算，把指数级问题压缩成多项式级问题。

2. 动态规划如何把“全局最优”拆成小决策

2.1 最优子结构：末尾字符只有三种可能

动态规划能成立，依赖一个叫最优子结构的性质。用大白话说就是：一个全局最优对齐，切掉最后一列之后，剩下的部分也必须是那两条前缀的最优对齐；否则，如果前缀部分存在得分更高的对齐，把它替换回来，整体得分还能更高，原方案就不是最优了。

于是问题变成：当我们考虑 A 的前 i 个字符和 B 的前 j 个字符时，最后一列到底怎么配？答案永远只有三种可能：A 的第 i 个字符对着 B 的第 j 个字符；A 的第 i 个字符对着空位；B 的第 j 个字符对着空位。没有第四种情况，因为任何对齐的最后一列必须有东西，而且只有这一个位置。

这个“只有三种可能”看起来简单到近乎废话，但整个递推就建立在这三选一上。

2.2 状态定义与状态转移方程

定义 dp[i][j] 为 A 的前 i 个字符与 B 的前 j 个字符能得到的最优比对得分。注意 i 和 j 从 0 开始，dp[0][0] 表示两条空前缀，得分为 0。那么 dp[i][j] 就是下面三个候选的最大值：

dp[i][j] = max( dp[i-1][j-1] + score(A[i-1], B[j-1]), # 对角：字符对齐字符 dp[i-1][j] + gap, # 上方：A[i-1] 对齐空位 dp[i][j-1] + gap # 左方：B[j-1] 对齐空位 )

这里 gap 是空位罚分，通常是一个负数。score 是匹配/错配的得分，比如相同 +2、不同 -1。之所以取 max 而不是别的，是因为我们寻找的是分数最高的一条路径，而每个子问题只需要记住最优值，后续所有决策只用这个最优值就够了。

这三个方向的名字需要记牢，因为它们和后面的回溯直接挂钩。从 dp[i-1][j-1] 过来，表示当前列是字符对字符；从 dp[i-1][j] 过来，表示 A 的第 i 个字符多出来了，对着一个空位；从 dp[i][j-1] 过来，表示 B 的第 j 个字符多出来了，对着一个空位。

2.3 第一行第一列：地基决定上层建筑

递推公式里的 i-1、j-1 都要求有边界值，所以 dp[i][0] 和 dp[0][j] 必须手动初始化，不能靠公式算。dp[i][0] 表示 B 是空序列，A 的 i 个字符只能全部对齐空位，于是得分是 i × gap；同理 dp[0][j] = j × gap。

有人会想：为什么不从 0 开始让边界一直等于 0？因为空位罚分是负的，初始化成 0 等于默认“缺失的前缀不算分”，会让全局高分方案大量；少初始化一行，边界状态就会乱掉，整个矩阵算出来都是错的。这块地基没有打好，后面每一格都可能计算错误，而且回溯时更容易出问题。

3. 手工推演：ACGT 与 AGT 的完整算例

3.1 打分规则

为了能手动复现，我们先把打分规则定成最经典的一组：匹配 +2，错配 -1，空位罚分 -2。这个组合常用在两个长度相近的 DNA 短序列上，后面第 5 节会讲参数对结果的影响。

3.2 填充得分矩阵

令 A = ACGT，B = AGT。用行表示 A 的前缀，用列表示 B 的前缀，那么矩阵是 5 行 4 列。先用初始化公式把第 0 行和第 0 列填好：

空 A G T 空 0 -2 -4 -6 A -2 ? ? ? C -4 ? ? ? G -6 ? ? ? T -8 ? ? ?

接下来逐个填。看 dp[1][1]：A[0]=A，B[0]=A，匹配。对角候选是 dp[0][0]+2 = 2；上方候选是 dp[0][1]+(-2) = -4；左方候选是 dp[1][0]+(-2) = -4。取最大值 2。

再看 dp[1][2]：A[0]=A，B[1]=G，错配。对角候选是 dp[0][1]-1 = -3；上方候选是 dp[0][2]-2 = -6；左方候选是 dp[1][1]-2 = 0。这里三个候选来源有竞争，真正胜出的是左方，说明从“A 与 A 已经配好”这个状态出发，把 B 的 G 对齐一个空位是当前局部最优。

继续填到中间，看 dp[2][2]：A[1]=C，B[1]=G，错配。对角候选是 dp[1][1]-1 = 1；上方候选是 dp[1][2]-2 = -2；左方候选是 dp[2][1]-2 = -2。最终取对角 1，说明此时“C 对 G 错配”虽然扣分，但比在 G 前再开一个空位更划算。

最后一个格子 dp[4][3]：A[3]=T，B[2]=T，匹配。对角候选是 dp[3][2]+2 = 4；上方候选是 dp[3][3]-2 = -2；左方候选是 dp[4][2]-2 = -2。取 4。整个矩阵填完后是这样：

空 A G T 空 0 -2 -4 -6 A -2 2 0 -2 C -4 0 1 -1 G -6 -2 2 0 T -8 -4 0 4

右下角 4 就是全局最优得分。

3.3 回溯：从右下角反推回去

算分容易，真正让很多人困惑的是回溯。回溯就是从右下角开始，每一步看当前格子的值到底来自三个候选中的哪一个，然后沿着来源方向往左上走，同时把对应的字符/空位记录下来。

从 dp[4][3]=4 出发，4 来自对角 dp[3][2]+2，所以这一列是 T/T，走到 (3,2)。dp[3][2]=2，来自对角 dp[2][1]+2，所以是 G/G，走到 (2,1)。dp[2][1]=0，候选里对角是 -3，左方是 -6，上方是 dp[1][1]+(-2)=0，所以它来自上方。上方表示 A[1]=C 对齐空位，于是记录 C/-，走到 (1,1)。dp[1][1]=2，来自对角 dp[0][0]+2，记录 A/A，走到 (0,0)。

把所有记录倒过来，得到对齐结果：

A C G T A - G T

用打分规则验算：A/A +2，C/- -2，G/G +2，T/T +2，总分 4，和矩阵右下角完全一致。到这里，最优对齐的“算”和“读”已经打通了。

3.4 结果说明了什么

这个例子里，C 被解释成 B 序列中的一段插入缺失，而不是三个错配。原因很直观：把 C 当空位只罚 -2，但若让 C 分别去对 A、G、T，会累积三个错配共 -3，更何况还会打断后续 G、T 的匹配。这就是动态规划自动权衡出来的结果。

如果换一组打分，例如把错配罚分改成 -1 而空位改成 -4，也许最优解就变成 C/G 错配而不是插入空位了。参数不同，最优路径就不同，这是序列联配里最需要记住的一点。

4. 落地成代码：Python 完整实现与细节坑

4.1 基础版本：一张得分表加一张方向表

手工算完后，代码不过是把上面的过程翻译成循环。实践中为了回溯方便，除了 dp 矩阵，通常还会同步维护一个 trace 矩阵，方向用 0 表示对角、1 表示上方、2 表示左方。

def global_align(seq_a, seq_b, match=2, mismatch=-1, gap=-2): """ 全局最优对齐。 注意：这里的 gap 直接传负值，例如 -2。 返回得分矩阵 dp 和方向矩阵 trace。 """ n, m = len(seq_a), len(seq_b) dp = [[0] * (m + 1) for _ in range(n + 1)] trace = [[0] * (m + 1) for _ in range(n + 1)] # 初始化第 0 列：B 为空序列，A 的每个字符都对齐空位 for i in range(1, n + 1): dp[i][0] = dp[i - 1][0] + gap trace[i][0] = 1 # 1 表示从上方来 # 初始化第 0 行：A 为空序列，B 的每个字符都对齐空位 for j in range(1, m + 1): dp[0][j] = dp[0][j - 1] + gap trace[0][j] = 2 # 2 表示从左方来 # 填充剩余格子 for i in range(1, n + 1): for j in range(1, m + 1): score = match if seq_a[i - 1] == seq_b[j - 1] else mismatch diag = dp[i - 1][j - 1] + score up = dp[i - 1][j] + gap left = dp[i][j - 1] + gap best = max(diag, up, left) dp[i][j] = best if best == diag: trace[i][j] = 0 # 字符对齐字符 elif best == up: trace[i][j] = 1 # seq_a[i-1] 对齐空位 else: trace[i][j] = 2 # seq_b[j-1] 对齐空位 return dp, trace

这里有一个容易被忽略的设计：当 diag、up、left 三个值相等时，代码优先选择对角。这个优先级会影响并列最优解的最终结果，后面会专门展开。

4.2 回溯代码

方向矩阵有了，回溯就是沿着方向走回去：

def backtrack(seq_a, seq_b, trace): i, j = len(seq_a), len(seq_b) align_a, align_b = [], [] while i > 0 and j > 0: if trace[i][j] == 0: align_a.append(seq_a[i - 1]) align_b.append(seq_b[j - 1]) i -= 1 j -= 1 elif trace[i][j] == 1: align_a.append(seq_a[i - 1]) align_b.append('-') i -= 1 else: align_a.append('-') align_b.append(seq_b[j - 1]) j -= 1 # 如果 i 或 j 还没走完，把剩余的字符补到空位后面 while i > 0: align_a.append(seq_a[i - 1]) align_b.append('-') i -= 1 while j > 0: align_a.append('-') align_b.append(seq_b[j - 1]) j -= 1 return ''.join(reversed(align_a)), ''.join(reversed(align_b))

调用一下：

dp, trace = global_align("ACGT", "AGT") a, b = backtrack("ACGT", "AGT", trace) print(a) print(b) print(dp[-1][-1])

输出：

ACGT A-GT 4

这里 A-GT 就是第二行对齐结果，第一行是 ACGT，逐列看就是 A/A、C/-、G/G、T/T，和手工推演的结论一样。

4.3 新手最容易踩的三个坑

第一个坑是边界处理。很多人写 while i > 0 and j > 0 循环结束后，忘了补上剩余的 i 或 j，结果回溯出来的序列开头少了一大截。比如 B 明显比 A 短的时候，A 开头那几个多出来的字符就会在倒数第二步被丢掉。一定要记得把两个 while 补字符循环加上。

第二个坑是方向矩阵初始化。trace 的第 0 行和第 0 列如果不赋值，默认值是 0，回溯时一到边界就会被当成“对角路径”，继续让 i 和 j 同时减一，轻则路径错误，重则越界或死循环。初始化的顺序也很重要：先填 dp 边界，再填中间格子，方向矩阵也要同步填。

第三个坑是并列最优解。比如某个格子 diag、up、left 有两个甚至三个值相同，代码只能记一个方向。这会导致回溯结果不是唯一，但得分仍然是最优的。如果你需要的是“所有”最优对齐，而不是“某一个”，那就要额外记录多个来源，或者用递归树去枚举。实际分析中，这一点往往比想象中更重要，尤其是碰到重复序列区域的时候。

5. 从矩阵走向真实数据：参数、变体与内存算账

5.1 罚分参数不是随便拍的

很多人刚学会跑通代码，就会问：匹配 +2、错配 -1、空位 -2 是哪来的？答案是，这只是一个经典默认组合，不一定是你的数据最适合的组合。真实的替换矩阵里，DNA 和蛋白质的打分策略完全不同；蛋白序列通常会用更细的替换矩阵，不同氨基酸之间的替换代价不一样，而不是简单分“相同/不同”两类。

空位罚分更要小心。罚分太重，真正的插入缺失被强行拆成若干个错配，序列之间的真实同源关系会被掩盖；罚分太轻，对齐里到处都是空位，把本来连续的匹配区域打碎，结果同样没法看。我自己的习惯是：先拿一小撮人工标注过的高置信序列跑几组参数，比较输出对齐的连续性和得分分布，再决定正式分析用哪一套。不要迷信任何“通用参数”。

5.2 全局对齐和局部对齐的岔路口

这篇讲的是全局最优对齐，要求两条序列从头到尾都比完。但真实场景里，很多时候你并不关心两端对不齐的长尾巴，只想知道中间有没有一段高度相似的区域。比如从一个长序列里找一个短 motif，全局对齐会把两端的无关区域也拉进来扣分，反而淹没了真正的信号。

这时候要用局部对齐。核心改动只有一个：状态转移的四个候选中多了一个 0。也就是说，一旦前缀得分变成负数，索性从当前这个位置重新开始算，不再背着前面欠下的分。这个“归零”的改动，解决的问题完全不同。下一篇如果要写第二部分，我最想展开的就是这个，以及它对应的空位罚分变体。

5.3 空间复杂度：拿到分数容易，回溯难

二维 dp 矩阵的时间复杂度是 O(nm)，空间复杂度也是 O(nm)。对于几千碱基的序列，这完全没问题；但人类基因组级别的序列或者长读长数据，nm 动辄上亿甚至更多，一张完整矩阵根本存不下。

如果只想要最优得分，非常简单：滚动数组，只保留上一行和当前行，空间降到 O(m)。但代价是方向信息没了，回溯做不了。要同时保住低内存和回溯，业界常用的思路是分治：在中间某一列算出一个交叉点，把问题切成左右两半递归去做。空间降到线性，时间仍然维持在 O(nm) 量级。当你开始处理长序列时，这个方向值得专门去啃。

5.4 多个最优解：你的回溯优先级早就替你做了决定

前面代码里，当三个候选值相同时我默认选对角。换一个方向优先级，比如优先上方、优先左方，结果往往不同。以一段低复杂度序列为例，比如 AAAA 和 AA，无论空位插在哪里，得分都一样，但不同工具会输出不同的空位位置。

这会造成什么实际影响？如果你接下来基于对齐坐标做下游分析，比如统计突变位点、抽取保守区域，空位位置不同，坐标就不同，结论可能有细微偏差。很多工具不会主动告诉你它的方向优先级，但你在阅读输出时，最好知道这一点。遇到并列最优的区域，至少要用不同的回溯策略跑一遍，看看结论是否稳定。

我自己在实际项目里已经很少手写这个二维矩阵了，但每次调试比对工具或者读别人的源码，都还是要回到这张表上。把状态转移和回溯这两条线彻底理顺之后，后面看半全局对齐、局部对齐、仿射空位罚分，都会顺畅很多。下一篇我打算接着写局部比对的“归零”逻辑和它在长序列里的表现，到时候见。