原文：序列比对那点事儿-生物知识学习 (biotechknowledgestudy.com)

序列比对

关于序列比对，首先我们来完成一个任务，手工比对一下ATGGCGT和ATGAGT：

我们是如何得出这样的结果的呢？可能的思路应该是从第一个字母开始，查找两个序列一致的部分，然后分割成小段，把一致的部分都对齐，然后把不一致的部分补齐。

但如果需要使用计算机来完成这一过程，似乎有些不太合理，因为计算量的问题。好吧，简单的，计算机是如何来解决这一问题的呢？

首先我们设定几个比对结果的准则，

两序列比对的结果必须遵守它们原有的序列顺序。

可以引入空格。空格的意思就是在其中一个序列当中插入一个空的字符，表示该位置相对的另一序列的碱基对应一个空碱基。

空格对空格是没有意义的，也是不被允许的。

有了上面的几条，接下来的思路还是评分体系。如果两序列对应的位置一致，那就得1分，如果是空格对碱基，那就不得分，如果完全不匹配，那就得-1分。总得分就是把所有的得分都加起来。最高分的就是最好的比对结果。

比如说：

以上两个结果来说，当然是4分的那个好些了。

这用计算机实现起来就是，1，因为两序列ATGAGT比ATGGCGT要少一个碱基，所以给少的那个加入一个空格。2，移动空格，计算比对得分。3，取得其中的最高分，做为比对结果。

但是，这样的比对总是需要去比较两者相应位置碱基是否一致。在速度上慢了。如果有个索引值的话就会快一些。于是对于DNA来说人们就建立了一个4×4的评分表（如果是蛋白质的话就是20×20）。

得分表：

评分表的使用，类似map[s1[i]][s2[i]]这样就可以拿到评分，速度自然会快一些。

这样的得分表效果还是太差，因为碱基和碱基还是不一样，其中A,G是嘌呤，T,C是嘧啶，所以不应该说C对T和C对A是一个评分。于是人们又对得分表进行的改进：

当然还有更复杂的，就是依据统计结果用马尔可夫链来生成评分表，这个马尔可夫链可以是一维的，就是说我们只考虑CTAG每个字符出现的概率。当然也可以是二维的，我们需要区分在A后面出现A可能性，出现G的可能性，出现T的可能性，出现C的可能性，等等。

对于蛋白质的评分表，先后出过很多，比如PAM(point accepted mutation)表(包括PAM250, PAM120等等)，BLOSUM(Blocks substitution matrix)表(包括BLOSUM62, BLOSUM50等等)。BLOSUM相对较新，也是相对比较好的一个表。对于蛋白质评分表来说，需要考虑的事情就比较多了，比如从氨基酸的物理化学特性，大小，极性，电核，亲水性等等。

这个评分表似乎没有考虑空格啊。我们把空格就想象成一个碱基的插入或者敲除。其生物学意义是移动了DNA开放阅读框(open reading frame)。所以，空格当然是非常不欢迎的。所以有一个原则，那就是如果需要插入空格的话，空格最好是比较集中地插入，而不是分散地插入。

有了评分表，我们分两步走来比对序列，第一步，生成所有可能的比对。第二步，计算得分并选出最高分。

那么，生成所有可能的比对应该是多少对呢？对于两条长度为Ｎ的序列来说，可能的比对有2^2N/sqrt(PI*N)那么多。比如长度为250的序列，可能的比对就是约10^149之多。好象是太多了点啊。

为了降低这个可能的比对数，还能不错过最佳的比对结果，于是就有了两种非常聪明的算法：Needleman-Wunsch算法(1970年，准确的讲是1969年提出的,70年发的paper)及在此算法基础上演变而来的Smith-Waterman算法(1981年)。这两个算法都涉汲到动态编程。其思路就是我前面讲的人的正常思维，就是先把大问题划成小问题，小问题解决了之后来拼起来解决大问题。具体到序列 上就是先把大序列从头看成小片段，在优化了小片段的结果之后，再逐步得出整个比对结果。好象是很难懂的样子，其实前面的字符串比对就已经用到了这个思路了，只是这里给出了一个名字而已。所以很多生物方面的进步，真的只是将其它领域的成功移植到生物领域而已。

首先讲一下Needleman-Wunsch算法,分两步走，第一步，构建得分表，第二步，从得分表逆向找回比对结果(traceback)。第一步，和之前写的levenshtein distance大同小异，而第二步就是前文没有提到过的。

在构建得分表的过程中，基本步骤是，

获取字符串s1的长度为n，获取字符串s2的长度为m，m!=0,n!=0。

假设空格的罚分为d。初始化一个二维数组D[m+1,n+1]，将第一行填充为0*d,..,j*d,..,n*d，将第一列填充为0*d,..,i*d,..,m*d。注意这里的填充不是0。

循环比对s1的每一个字符与s2的每一个字符。对于字符s[j](j取值1至n)与s2[i](i取值1至m)。这时比较二维数组D中的三个值：D[i-1,j]+d（左侧值）, D[i,j-1]+d（顶部值）, d[i-1,j-1]+s[s2[i]][s1[j]]（左上角）；这三个值分别代表删除，插入或者替换。取三个值中的最大值填充D[i,j]。注意是最大值。上面的s[s2[i]][s1[j]]就是指的ATGC的评分表中对应的值。

接下来的步骤是traceback，找回最佳的比对结果。从D矩阵最右下角出发，一直找到最左上角。移动的过程有三个可能的方向，左，上，以及左上。在这三个可能的方向上找出最大值。然后前进至下一步。

关于traceback，出示几个图可能更清楚一点，同时也可以进一步清晰为什么得分表就可以得出是删除还是插入了。我们在计算得分表时，总是取左侧+d，顶部+d，以及左上角+s的三个值中最小值，如同前面已经讲过的，

左上角+s是要比较的两字母相同或者相异，用diag表示

左侧+d就是由左侧序列引入的空格，用left表示

顶部+d就是由顶部序列引入的空格。用up表示

假设我们要比对两个字符串：SEND以及AND。

从最右下角开始，它的值为diag，那么就是最末的两个字母对齐，我们得到

D

D

并得到向左上回溯的方向。左上角又得到一个diag，那再得到次末两个字母对齐，

ND

ND

并得到向左上角回溯。左上角得到一个left，那我们得到在测序列需要引入一个空格，

END

-ND

并得到向左侧回溯。左侧得到一个diag，我们得到两字母对齐，

SEND

A-ND

再向左上角得到done，完成比对的traceback。

接下来就是Smith-Waterman算法了。它与Needleman-Wunsch的不同之处有三：

矩阵最初始的两条边（最顶部及最左侧的值）全部以零填充。

之前讲的三个值中取最大值，值不能低于零。否则不计入得分表。

回溯时选择最大值，直到遇到0为止。

其中第二条要注释一下。因为之前一直讲取最小值，这里为什么变成取最大值了呢？主要是评分标准变化了。我们假设一个简单的评分标准：

空格：－1

相同：＋2

不同：－1

我们发现，这里的评分规则和之前的正好相反，值越高表明两字母越相类。

那么SW算法比之NW算法有什么提高呢？NW又称全局算法，SW又称局部算法。因为生物序列有太多的突变，这使得引入突变的部分会增加很多的不确定性（噪音）。而SW算法就避免了突变的部分同未突变的部分做相同的考虑，而是集中精力于得分值比较高的部分，也就是相同的部分。这就更加类似人脑在比对当中的做法，先优化的是小片段，然后再总体拼合，而不是逐一延伸比对。而这一算法也是基于这一算法出现之前统计模型(by Karlin and altschul)。

SW算法比NW算法要更快一些。

如果你想了解更多的算法方面的知识，可从这本文的这些基础知识出发，继续了解blast，blat，bowtie等是如何针对自己的问题而提高比对效率的。