跳转到内容

半全局比对

快速概览

半全局比对允许序列端部存在gap,是理解read mapping和重叠检测等场景的关键算法,介于全局比对和局部比对之间。

  • 重点理解半全局比对的边界条件如何与全局比对不同
  • 半全局比对在很多实际工具中是比全局比对更实用的选择
所属板块 核心方法

索引、比对、组装与概率模型构成的核心算法主轴。

阅读目标 帮助建立阅读上下文

先判断这页与你当前问题的关系,再决定是否深入展开。

建议前置 先建立相关基础对象与方法直觉

建议先建立相关基础对象与方法直觉,再进入本页。

半全局比对(semi-global alignment)是序列比对的一种变体,它允许序列的端部(开头或结尾)存在gap而不产生罚分。

与全局比对不同:

  • 全局比对:要求两条序列从开头到结尾完全对齐,端部gap也要罚分
  • 局部比对:只寻找最优的局部相似片段
  • 半全局比对:允许序列端部有gap,但内部gap仍需罚分

在生物信息学中,半全局比对的价值体现在:

  • Read mapping:reads通常比参考序列短,read的端部不需要与参考序列的端部对齐
  • 重叠检测:在组装中检测reads之间的重叠时,reads的端部可能有未覆盖区域
  • 序列拼接:当两条序列有部分重叠时,半全局比对能更好地识别重叠区域
  • 实际工具:很多现代mapper的局部比对阶段使用了半全局比对的思想

全局比对要求两条序列完全对齐:

全局比对:
S1: A T G C T A G
S2: - T G C - A G
^ 开头gap要罚分
S1: A T G C T A G
S2: A T G C T A -
^ 结尾gap要罚分

半全局比对允许端部gap不罚分:

半全局比对:
S1: A T G C T A G
S2: - T G C - A G
^ 开头gap不罚分
S1: A T G C T A G
S2: A T G C T A -
^ 结尾gap不罚分

根据允许端部gap的位置,半全局比对有三种常见变体:

  1. 允许S1开头gap:适用于S1是S2的子序列情况
  2. 允许S2开头gap:适用于S2是S1的子序列情况
  3. 允许两端gap:适用于检测重叠情况
状态定义
与全局比对相同,`dp[i][j]`表示S1前i个字符与S2前j个字符的最优比对得分。
边界条件
与全局比对不同,某些边界条件设为0而非线性递增。
最优解
根据具体变体,可能在最后一行、最后一列或矩阵中寻找最大值。

递推关系与全局比对相同:

dp[i][j]=max{dp[i1][j1]+s(xi,yj)匹配或替换dp[i1][j]+gS1插入gapdp[i][j1]+gS2插入gap dp[i][j] = \max \begin{cases} dp[i-1][j-1] + s(x_i, y_j) & \text{匹配或替换} \\ dp[i-1][j] + g & \text{S1插入gap} \\ dp[i][j-1] + g & \text{S2插入gap} \end{cases}

其中:

  • s(xi,yj)s(x_i, y_j) 是替换得分
  • gg 是gap罚分

不同变体的边界条件不同:

dp[0][j]=0(S1开头gap不罚分)dp[i][0]=ig(S2开头gap要罚分)\begin{aligned} dp[0][j] &= 0 \quad \text{(S1开头gap不罚分)} \\ dp[i][0] &= i \cdot g \quad \text{(S2开头gap要罚分)} \end{aligned} dp[0][j]=jg(S1开头gap要罚分)dp[i][0]=0(S2开头gap不罚分)\begin{aligned} dp[0][j] &= j \cdot g \quad \text{(S1开头gap要罚分)} \\ dp[i][0] &= 0 \quad \text{(S2开头gap不罚分)} \end{aligned} dp[0][j]=0(S1开头gap不罚分)dp[i][0]=0(S2开头gap不罚分)\begin{aligned} dp[0][j] &= 0 \quad \text{(S1开头gap不罚分)} \\ dp[i][0] &= 0 \quad \text{(S2开头gap不罚分)} \end{aligned}

根据变体不同,最优解的位置也不同:

  • 变体1:在最后一行 dp[m][j] 中寻找最大值
  • 变体2:在最后一列 dp[i][n] 中寻找最大值
  • 变体3:在最后一行或最后一列中寻找最大值

算法1:半全局比对(允许两端gap)

输入:序列 S1 = x_1x_2...x_m,S2 = y_1y_2...y_n,替换得分函数 s,gap罚分 g
输出:最优比对得分
1. 初始化 dp[m+1][n+1]
2. // 边界条件:两端gap不罚分
3. for i = 0 to m:
dp[i][0] = 0
4. for j = 0 to n:
dp[0][j] = 0
5. // 填充动态规划矩阵
6. for i = 1 to m:
for j = 1 to n:
dp[i][j] = max(
dp[i-1][j-1] + s(x_i, y_j),
dp[i-1][j] + g,
dp[i][j-1] + g
)
7. // 在最后一行和最后一列中寻找最大值
8. best_score = max(max(dp[m][j] for j in 0..n), max(dp[i][n] for i in 0..m))
9. return best_score

时间复杂度O(mn)O(mn)

空间复杂度O(mn)O(mn)(可优化到 O(min(m,n))O(\min(m,n))

考虑序列 S1 = "ATGCC"S2 = "TGCCG",使用简单得分:

  • 匹配:+1
  • 失配:-1
  • gap:-2

使用允许两端gap的半全局比对:

ε T G C C G
ε 0 0 0 0 0 0
A 0
T 0
G 0
C 0
C 0

逐步计算(以变体3:允许两端gap为例):

ε T G C C G
ε 0 0 0 0 0 0
A 0 -1 -2 -3 -4 -5
T 0 1 0 -1 -2 -3
G 0 0 2 0 -1 -2
C 0 -1 0 3 1 0
C 0 -2 -1 1 4 2

计算验证示例 dp[2][2](T对G):

  • 来自对角线:dp[1][1] + s(T,G) = 0 + (-1) = -1
  • 来自上方:dp[1][2] + g = 0 + (-2) = -2
  • 来自左侧:dp[2][1] + g = 0 + (-2) = -2
  • 取最大值:0(因为所有选项都为负,但半全局边界初始化为0)

实际应为 max(0, -1, -2, -2) = 0,矩阵已修正。

在最后一行和最后一列中寻找最大值:

  • 最后一行:[0, -2, -1, 1, 4, 2],最大值为 4(位置 dp[5][4],即最后一个C对C)
  • 最后一列:[0, -5, -3, -2, 0, 2],最大值为 2

全局最优得分为 4

dp[5][4]=4 回溯(C-C匹配):

  • dp[5][4]=4dp[4][3] + s(C,C) = 1 + 1 = 2?不符
  • 检查其他来源…

实际上,观察递推关系,最优路径可能是:

S1: A T G C C
- T G C C ← 注意:S2末尾还有G
或简化为只显示重叠区域:
S1: - T G C C
S2: T G C C

S1开头的’A’和S2末尾的’G’形成端部gap,不产生罚分。

验证比对得分:T-T(+1) + G-G(+1) + C-C(+1) + C-C(+1) = 4

当将read比对到参考基因组时:

  • read通常比参考序列短
  • read的端部不需要与参考序列的端部对齐
  • 使用半全局比对可以避免read端部的不必要罚分

在组装中检测reads之间的重叠:

  • 两条read可能有部分重叠
  • 非重叠部分应该在端部
  • 半全局比对能更好地识别重叠区域

当拼接两条有重叠的序列时:

  • 重叠区域应该对齐
  • 非重叠的端部不需要对齐
  • 半全局比对提供了合适的模型
比对类型端部gap处理典型应用
全局比对端部gap要罚分比较两条完整序列
局部比对可以从任意位置开始/结束寻找相似片段
半全局比对端部gap不罚分read mapping、重叠检测

半全局比对可以看作:

  • 全局比对的推广:放宽了端部gap的限制
  • 局部比对的特例:限制了比对必须覆盖至少一条序列的大部分

与全局比对类似,可以使用滚动数组将空间复杂度从 O(mn)O(mn) 降低到 O(min(m,n))O(\min(m,n))

如果两条序列预期高度相似,可以使用带状动态规划,将复杂度降低到 O(kmin(m,n))O(k \cdot \min(m,n)),其中 kk 是带宽。

  • Durbin, R., et al. Biological Sequence Analysis
  • Gusfield, D. Algorithms on Strings, Trees, and Sequences
  • Pearson, W. R., & Miller, W. (1992). “Dynamic programming algorithms for biological sequence comparison”