序列每天从0开始_从零开始生物信息学(3):序列比对-Smith–Waterman算法

前言

1. 在上一部分我们已经讲解了生物信息学序列比对中一个非常经典的算法Needleman-Wunsch算法，没看过的童鞋可见此处这是一个基于全局匹配的动态规划算法。这一部分主要讲下基于Needleman-Wunsch算法扩展得到的Smith-Waterman算法，这是一个基于局部匹配的动态规划算法，后来常用的序列比对算法例如FASTA和大名鼎鼎的BLAST也是基于这个算法改进的。

Smith-Waterman算法

Smith-Waterman算法是由Temple F. Smith and Michael S. Waterman 在1981年提出来的，比Needleman-Wunsch算法晚了11年。

Smith-Waterman算法的整个计算流程和Needleman-Wunsch算法很接近，在几处细节上有改动。

为了方便大家理解，还是举个例子来讲解下Smith-Waterman算法原理。假如现在我们有两个序列：

TGTTACGG
GGTTGACTA

同样的，我们定义得分准则如下：

如果两个位点的碱基匹配（MATCH），得3分，出现不匹配(MISMATCH)或者为-3分。并且如果出现空缺（GAP）惩罚分数做一个线性增长惩罚：Wk = 2k (k指的是出现GAP的次数，初始值为2，也就是每次出现一个GAP惩罚的分数依次为-2，-4，-6，递增)。

惩罚得分矩阵

和Needleman-Wunsch算法相似，我们也需要构造一个得分矩阵用于Smith-Waterman算法的比对序列回溯，矩阵的行列同样是是两个序列的碱基排列，第一行和第一列为置0，这是和Needleman-Wunsch算法不同的地方。

然后同样地从左上到右下的顺序计算每个位点的得分，每个位点的得分是与它位置的上面，左边和左上角三个位置的得分相关。最后取三者的最大值，具体计算如下：

其中s(ai, bj)函数就是我们上面定义的匹配或者不匹配两种情况，W是横方向或者纵方向的GAP惩罚。最后计算得分方式如下：

三个方向的得分=该方向上一位点得分+移动过程得分

总得分=max(三个方向的得分， 0)

为啥要和0比较呢，可以认为就是这里的得分为负了，起到一个重置的作用，以此循环从上到下，从左到右得到整个矩阵的得分。

例如，第三个图的对应位置的左上角为0分，但是对应横方向和纵方向的G是匹配的，所以斜方向得分为0+3，横竖两个W都是0-2=-2，所以总的得分就是max(3，-2，-2, 0)=3.其他的依次类推，最后得到整个惩罚矩阵：

红色箭头是每个位点最高值的来源，然后我们可以发现整个数组有个最大值（蓝色标记），我们需要从最大值处回溯得到整个比对序列。同样的，从矩阵的最大值开始，比较每个当前位点的左上，上方和左方三个方向，如果最大值出现在上面，则横向这条序列引入一个GAP ("-")，纵向这条序列取该处碱基；如果最大值出现在左边，则纵向这条序列引入一个GAP ("-")，横向这条序列取该处碱基; 如果最大值出现在左上角，则不引入GAP，纵向和横向均取该处碱基。值得注意的是，回溯的过程在回溯的值为0就停止，最后将整个序列翻转，就可以得到最终的序列比对结果。

最终的局部序列比对结果如下：

G T T - A C
| | |   | |
G T T G A C

Smith-Waterman VS Needleman-Wunsch

首先，我们需要知道我们需要搞清楚这两个算法的核心区别，Smith-Waterman算法针对的序列的局部匹配性。

对于两段序列而言，即使他们的不相关，或者不具备相似性，他们依旧可以拥有相似的局部序列分布，并且这个局部序列矩阵得分会高于全局的矩阵得分（从矩阵的最大值不是最右下角的值就可以看出），如果我们能在他们之间比对出一些局部片段，具有高度的相似性，一定程度上，我们可以得到两个序列具备同源性，也就是可能拥有相同的祖先，对于一些外显子序列，有可能这些片段会转录出一些相似的RNA及蛋白质。还有就是这种比对方法可会揭示一些匹配的序列段，而本来这些序列段是被一些完全不相关的残基所淹没。

因此，局部比对的相似性往往比全局比对相似性更具有价值。

全局比对VS局部比对

局部匹配与惩罚矩阵

其实，之前一直困扰我的还有为什么每个位点的得分都需要和0比较。后来我想了下，应该是0表示两个序列相似性的一个临界值，如果到当前位点，两个序列的相似性都为0或者以下，说明GAP和MISMATCH的情况是一定程度是多余MATCH的，那么间接的理解就是之前的序列相关性较低，可以认为是不相关，所以当是负数时，最大值直接置0，也就是重新开始序列相似度计算。这也一定程度上说明为什么直接回溯到0就停止了。

最后总结下两种算法计算方式的差异：

这次就到这里，有什么问题，也欢迎和我联系！

欢迎大家关注我的知乎专栏：从零开始生物信息学

相同内容也可以关注我的微信公众号: 壹读基因：