序列比对原理

序列比对是非常基础的工作，生物信息学分析大多涉及到这一步骤，虽然大部分人不做相关算法或软件的开发，但是了解总体的算法原理还是需要的。
本文主要分为两部分。第一部分介绍简单的 2 条序列比对(多重序列比对不考虑)，了解动态规划策略如何寻找全局比对的解；第二部分介绍二代测序的 BWT + FM-index 策略如何实现降低资源消耗下的大规模数据比对。

Needleman–Wunsch

Needleman–Wunsch 是用动态规划(dynamic programming)寻找最优序列全局比对的算法。他将大问题(如：全局序列比对)分解为小问题，通过最优化解决小问题实现大问题的最优解。对于给定的 2 个比对序列，通过计算得分和回溯矩阵，通过回溯路径得到最佳的全局比对。

假设需要比对的序列是 ATTGTC 和 ATCTC. 先定义匹配得分规则，假设

正确比对：得分 +1
错误比对：得分 -1
Indel: 得分 -1

定义得分规则是关键， 定义不同的得分规则将产生不同的比对结果。 实际操作得分规则的设定是很复杂的，不仅仅考虑匹配不匹配得分，还会不同的匹配赋予不同得分，indel 长度不同也赋予不同得分。下图是蛋白比对的 Blosum50 得分系统。

Blosum50

如下图所示，将 2 序列分别作为矩阵行和列表头。矩阵每个格子分数可能来源有 3 个，分别是对角线、左边和上面，分别计算这 3 种来源的分数，然后取最大值作为得分，且记下这个来源后续用于回溯。假设用 q 表示不同来源计算的本格得分，用 C(i, j) 表示第 i 行第 j 列格得分，用 S(i, j) 表示本格匹配得分，用 g 表示 gap 的惩罚分数。那么 3 种方法本格得分计算公式分别是：

$q^{diag} = C(i-1, j-1) + S(i, j)$
$q^{left} = c(i, j - 1) + g$
$q^{up} = c(i-1, j) + g$

第一行和第一列比较特殊，格子得分只能来源于一个方向，所以从第一行第一列格子分数为 0 可以先把第一行和第一列其余得分先计算。

FirstRowCol

然后考虑下图问号格子得分计算，这里开始格子得分就有 3 种来源，按照之前公式分别计算每一种得分然后取最高分。

3Path

按照公式计算：

$q^{diag} = 0 + 1 = 1$
$q^{left} = -1 -1 = -2$
$q^{up} = -1 -1 = -2$

所以问号格子分数取对角线来源得分，那么记下这个分数并且记录来源。

Diag

以此类推能计算出完整的得分和回溯矩阵。如下图所示。

ScoreMatrix

最后从右下角开始回溯，如果格子得分来源是对角线那么该位置匹配，如果来源于左侧那么左边序列有 gap, 如果来源于上方那么上方序列有 gap.
从最右下角的格子开始，该格子来源于对角线，属于匹配。

C
|
C

回溯到他的对角线，也是来源于对角线，那么也是匹配。

TC
||
TC

再往上回溯这时候有 2 种分数来源相等，其实等于说此时局部没有完全的最优解，那么此时可能有 2 种情况。

# 情况一，来源左侧，那么左侧序列有 gap
GTC
 ||
-TC

# 情况二，来源于对角线，那么属于匹配
GTC
 ||
CTC

以此类推如图中绿色格子所示，这里也刚好出现路径不唯一的情况，也即是这种方法没有找到唯一最优解。按照图中绿色回溯路径，总共有以下几种比对结果。

# 1
ATTGTC
||  ||
ATC-TC

# 2
ATTGTC
| | ||
A-TCTC

# 3
ATTGTC
||  ||
AT-CTC

在这里刚好出现多个最优解，主要是举例子的得分体系设置太简单了，设置好合理的得分体系加上现实一般比对序列比较长，不会这么容易出现多个最优解情况。

Smith-Waterman

Smith-Waterman 算法是从 Needleman-Wunsch 算法发展而来的适用于局部比对(local alignment)的算法。其计算得分依旧是取 3 个方向来源的最大值，但如果出现负分则修改为 0. 所以其第一行和第一列全部值为 0. 后面回溯也不从右下角回溯到左上角，而是从得分最高的格子回溯到第一个得分为 0 的格子。从而在长序列中寻找最佳匹配的子序列比对。这在生物学上很有用，寻找到高度相似的子序列，可能就是同源序列部分。

Word method

也被称为 k-tuple 方法，是一种启发式算法(heuristic)，不追求最优比对但是比动态规划效率高很多，尤其适用于在大规模数据库搜索匹配序列，大名鼎鼎的 BLAST 就是用此方法。
这种方法将要检索的序列分割成多个短序列(word)，然后用所有的短序列在数据库进行匹配，如果有短序列能匹配，那么将该序列往两边继续延申继续匹配，这样能快速取得局部比对的结果。

二代测序的比对

二代测序的比对面临特殊的的困难，一是 reads 数目多且方向未知，二是参考基因组长，三是有测序错误、基因突变、结构变异等。所以比对要求能检测出各种变异，还要求能尽可能低的资源消耗(时间，计算资源)下工程实现。
二代测序也有许多不同的策略方案，甚至有 alignment-free 方法，本文介绍 bowtie 和 BWA 用的 BWA + FM-index 方案总体原理。

BWT(Burrows-Wheeler Transform) 算法将字符串相同字符排列在相邻或相近位置，而且是可逆的，常常被用以压缩数据，比如 bzip2 软件。参考基因组往往很大(如人种约 30 亿碱基)，直接拿去比对消耗太大量计算机资源，BWT转换压缩能节约大量计算机资源。BWT 算法转换过程如下图所示。

BWT

给定一个字符串先进行循环右移，这里要加入一个不在原字符串的标识符 $(也可以用其他字符做标识符)，然后按照将所有的循环右移序列排序，取最后一列字符为输出结果(上图红色)。
“循环右移”听着很不直观，用代码帮助理解：

# 这里字符串 t 已经添加了标识符 $
# 把字符串复制连接后以字符串长度为窗口，滑窗取子字符串
def rotations(t):
    tt = t * 2
    return [ tt[i:i+len(t)] for i in xrange(0, len(t)) ]
    
# 进行排序
sort(rotations(t))

BWT 转换后对于除第一行的每一行，其开头字符是结尾字符的下一个字符，在上图分别用 F,L 2列表示开头与结尾字符，对于同一行 i 来说在原字符串的字符顺序是 L[i]F[i].
BWT 转换 F 列和 L 列是相关的。把他们的映射关系表示为 LF，C(c) 表示在 T 中小于 c 字符的字符数，Occ(c, i) 表示在 $T^{bwt}[1, i]$ 中字符 c 出现数目， $T^{bwt}[i]$ 字符在 F 列出现位置为：
$LF(i) = C(T^{bwt}[i]) + Occ(T^{bwt}[i], i)$
以下图序列为例，假设 i 为 3 那么 $T^{bwt}[3]$ 为字符 G，此时 C(G) 为 5 Occ(G, 3) 为 $T^{bwt}[1, 3]$ 中 G 字符出现数目——等于 1. 所以此时 LF(3) = 5 + 1 =6 ——最后一列的第 3 位 G 映射到第一列的第 6 位。

FM-index(full-text minute index) 允许压缩字符串同时可以快速检索子字符串。他可以高效查找到压缩文本的模式匹配字符串(pattern)出现次数以及相应的位置。
假设原字符串长度 n 需要检索的子字符串 P 长度为 m, 通过向后搜索的循环能找出排序矩阵(M)中 P[i, m] 为前缀的行数范围。循环算法 BackwardSearch 如下。

sp <- 1, ep <- n
for i <- m to 1 do
    sp <- C(P[i]) + Occ(P[i], sp - 1) + 1
    ep <- C(P[i]) + Occ(P[i], ep)
    if sp > ep then return Φ
end
return [sp, ep]

用上图序列为例搜索 GAC 子字符串其过程在上图用红色箭头展示，其循环迭代如下图。

BackwardSearch

因为矩阵(M)是排序的所以相同前缀的字符串是相邻的，这里 sp, ep 都为 6，所以 GAC 子字符串出现次数为 6 - 6 + 1 = 1 次。优化后这个向后搜索的时间复杂度与子字符串长度 m 线性相关。
确定子字符串出现次数后需要进行定位(在二代测序就是确定 reads 匹配序列在基因组的位置)。这是通过存储部分(标记的)矩阵M的行的文本偏移量，然后将需要定位的行(子字符串行)按照 LF 映射进行跳转，直到遇到标记的行，根据标记行位置和跳转次数确定子字符串的位置。
上面是完全匹配的子序列，实际测序的比对是允许部分错配的，假设允许错配的碱基数为 k 那么这个匹配检索的算法 kMismatchSearch 如下。

kMismatchSearch(P, k, j, sp, ep):

if sp > ep then return Φ
if j = 0 then return [sp, ep]
foreach c ∈ {A, C, G, T} do
    sp' <- C(c) + Occ(c, sp - 1) + 1
    ep' <- C(c) + Occ(c, ep)
    if P[j] != c then k' <- k - 1 else k' <- k
    if k' >= 0 then
        kMismatchSearch(P, k', j - 1, sp', ep')
    end
end

与之前的 BackwardSearch 主要区别是每一个位置的所有可能碱基都考虑，通过允许的错配次数 k 控制迭代。
实际的工程实现还有非常多的技术细节需要考虑和优化，本文不多介绍，比如 kMismatchSearch 在出现 gap 时表现很差，bowtie2 因此先将 reads 生成小的 seeds, 因为够短从而比较不会被 indel 影响，用这些 seeds 先去找可能的比对位置，再进行延伸到 reads 全长比对，实现 gap 友好同时速度较快的比对。

参考资料

Canzar, Stefan, and Steven L. Salzberg. "Short read mapping: An algorithmic tour." Proceedings of the IEEE 105.3 (2015): 436-458.
Langmead, Ben, and Steven L. Salzberg. "Fast gapped-read alignment with Bowtie 2." Nature methods 9.4 (2012): 357.
Sequence alignment - Wikipedia
Needleman–Wunsch Algorithm
Needleman–Wunsch algorithm - Wikipedia
用Python从头实现Needleman-Wunsch序列比对算法 - 知乎
 从零开始生物信息学(3):序列比对-Smith–Waterman算法 - 知乎
 BLAST (biotechnology) - Wikipedia
BWT (Burrows–Wheeler_transform)数据转换算法 - 旭东的博客 - 博客园
 FM-index - Wikipedia
blosum50
Teaching Materials |