史密斯-沃特曼算法

史密斯-沃特曼算法（Smith-Waterman algorithm）是一种进行局部序列比对（相对于全局比对）的算法，用于找出两个核苷酸序列或蛋白质序列之间的相似区域。该算法的目的不是进行全序列的比对，而是找出两个序列中具有高相似度的片段。该算法由坦普尔·史密斯和迈克尔·沃特曼于1981年提出。史密斯-沃特曼算法是尼德曼-翁施算法的一个变体，二者都是动态规划算法。这一算法的优势在于可以在给定的打分方法下找出两个序列的最优的局部比对（打分方法使用了置换矩阵和空位罚分）。该算法和尼德曼-翁施算法的主要区别在于该算法不存在负分（负分被替换为零），因此局部比对成为可能。回溯从分数最高的矩阵元素开始，直到遇到分数为零的元素停止。分数最高的局部比对结果在此过程中产生。在实际运用中，人们通常使用该算法的优化版本。生物学上一个基本问题是，一个基因或者蛋白是否和别的基因或者蛋白存在联系。蛋白质在序列层次的关联暗示了其同源性，同时暗示了它们可能有类似的功能。通过分析多个DNA或者蛋白质序列，我们可能会发现一些保守序列和区域。这种基因序列或者蛋白质序列之间关联性的分析就是通过序列比对来完成的。如今，人们完成了许多不同生物体的基因组分析，获得了大量的基因和蛋白质序列数据。序列比对能够向我们揭示某一生物体中蛋白质之间的关联以及蛋白质在不同生物体中的关联，从而进一步理解生命和进化。序列比对起源于1970年。Saul B. Needleman和Christian D. Wunsch于1970提出了一个启发式的同源算法，称为尼德曼-翁施算法。该算法使用迭代的方法计算出一个矩阵，从而达到全局比对的效果。一些其他的启发式分析方法也在70年代被提出。Sankoff，Reichert，Beyer等人提出了一些严谨的数学的算法用来分析序列。然而这些方法通常难以揭示序列的生物学意义。 Sellers提出了一套测量序列距离的方法。沃特曼等于1976年在该方法的基础上，增加了插入和删除任意长度片段的方法。该方法解决了如何用最小的基因突变的次数将一个序列转变成另一个序列的问题。之后，坦普尔·史密斯和迈克尔·沃特曼于1981年提出了史密斯-沃特曼算法，解决了局部比对问题。设要比对的两序列为 A = a 1 a 2 . . . a n {displaystyle A=a_{1}a_{2}...a_{n}} 和 B = b 1 b 2 . . . b m {displaystyle B=b_{1}b_{2}...b_{m}} ，其中 n {displaystyle n} 和 m {displaystyle m} 分别为序列 A {displaystyle A} 和 B {displaystyle B} 的长度。史密斯-沃特曼算法通过字母的匹配，删除和插入操作，把两条序列进行排列。实际上插入和删除都是引入空位的操作（在不同序列上）。序列1上某一位置的插入相当于在序列2上对应位置的删除。在实际计算中，只需要考虑何时引入空位即可。该算法主要分两步，计算得分矩阵和寻找最优比对序列。可以简单描述为：史密斯-沃特曼算法为局部比对算法，用于找出两序列中最相似的区域。尼德曼-翁施算法为全局比对算法，用于完整匹配两个序列。这两种算法的目不同，因此研究者在选择算法时要根据实际需求来决定。两种算法都使用了置换矩阵，空位罚分，得分矩阵，以及回溯的方法，具有一定的相似度。然而，史密斯-沃特曼算法与尼德曼-翁施算法的三个主要区别在于：其中，打分阶段分数不为负是非常重要的一点区别。它使得对序列局部的比对成为可能。当任何位置分数低于0，即表示此前的序列不具有相似性；而此时将此元素分数设为0，即达到了“重设”的效果，使得从此位置开始的比对不受之前位置的影响。初始化阶段的区别使得史密斯-沃特曼算法中从一序列任意位置开始匹配另一序列不受罚分，而尼德曼-翁施算法因为要匹配完整的序列，所以两端的空位也需要罚分。每一碱基对或残基对匹配或错配的分数通常用一矩阵表示，称为置换矩阵。基本的置换矩阵为匹配则加分，错配则扣分。以核苷酸序列为例，若匹配为+1，错配为-1，则置换矩阵为：该置换矩阵可以表示为： s ( a i , b j ) = { 1 , a i = b j − 1 , a i ≠ b j {displaystyle s(a_{i},b_{j})={begin{cases}1,quad a_{i}=b_{j}\-1,quad a_{i}neq b_{j}end{cases}}}不同的碱基对或残基对可以有不同的分数以适应于特定的场合。氨基酸序列比对的置换矩阵一般更加复杂。通常性质相似的残基对具有正分数，反之，不相似的具有负分数。此外，残基对的相似或不相似程度决定了分数的大小。参见PAM（英语：Point accepted mutation），BLOSUM（英语：BLOSUM）。当碱基对或残基对之间匹配会导致更低分数时，需要空位的引入，即让碱基或残基与空位匹配。空位罚分决定了插入或者删除的分值。最基本的空位罚分方式为每次插入或者删除的得分相同。然而在生物学上，两序列之间一般存在着具有不同相似度的区域。另外，单个基因突变事件可能导致一长串空位的插入。因此，一个连续的较长的空位优于多个分散的小的空位。虽然多个分散的小的空位可以产生更多匹配，但一个连续的较长的空位代表这个区域只在一个序列中出现，使用后者可以避免为了得到高分而过度匹配这段序列。实现该方法只需要引入空位起始罚分和空位延长罚分的概念。空位起始罚分通常高于空位延长罚分。例如EMBOSS Water的默认参数为空位起始罚分=10，空位延长罚分=0.5。实际应用中空位罚分方法有多种（参见空位罚分（英语：Gap_penalty）），此处详细介绍这两种最常见的方法。设 W k {displaystyle W_{k}} 为空位罚分函数，其中 k {displaystyle k} 为空位长度：W k = k W 1 {displaystyle W_{k}=kW_{1}}该模型空位的罚分正比于空位的长度。 W 1 {displaystyle W_{1}} 为单个空位的罚分。使用该模型时，原史密斯-沃特曼算法可简化为：H i j = max { H i − 1 , j − 1 + s ( a i , b j ) , H i − 1 , j − W 1 , H i , j − 1 − W 1 , 0 {displaystyle H_{ij}=max {begin{cases}H_{i-1,j-1}+s(a_{i},b_{j}),\H_{i-1,j}-W_{1},\H_{i,j-1}-W_{1},\0end{cases}}}此时该算法步数为 O ( m n ) {displaystyle O(mn)} ， m {displaystyle m} 和 n {displaystyle n} 分别为两序列的长度。在对某一位置考虑空位罚分时，只需考虑相邻的上方的和左侧的位置，而不需考虑整列和整行。W k = u ( k − 1 ) + v ( u > 0 , v > 0 ) {displaystyle W_{k}=u(k-1)+vquad (u>0,v>0)}该模型考虑空位起始罚分和空位延长罚分，其中 v {displaystyle v} 为开始空位的罚分， u {displaystyle u} 为每次延长空位的罚分。例如，一个长度为2的空位的罚分为 u + v {displaystyle u+v} 。未经优化的该模型的算法步数为 O ( m 2 n ) {displaystyle O(m^{2}n)} 。Gotoh将该算法的步数优化至 O ( m n ) {displaystyle O(mn)} ，然而，优化后的方法只能找出一个局部比对结果，并且并不保证能成功获得结果。Altschul进一步优化了Gotoh的算法，在保持复杂度不变的情况下解决了Gotoh算法的不足。Myers和Miller指出，Gotoh和Altschul的算法忽视了Hirschberg在1975年提出的在线性空间内进行序列比对的方法，并在改进的算法中运用了该方法。Myers和Miller的算法仅使用 O ( n ) {displaystyle O(n)} 空间， n {displaystyle n} 为较短的序列的长度。以两序列TACGGGCCCGCTAC和TAGCCCTATCGGTCA进行比对为例，分别考虑两种模型。当使用空位权值恒定模型时，即不考虑空位起始和延长的区别，得到如下结果（置换矩阵为DNAfull，空位起始和延长均为1.0，其他参数为默认）：当使用空位延伸罚分模型时，设空位起始罚分为5.0，空位延长罚分为1.0时，得到如下结果：可以看出，空位延伸罚分模型虽然总分略低，但考虑空位起始罚分和空位延长罚分能够有效避免多个小的空位。得分矩阵的作用是对两序列中两两位置进行打分以逐步记录最优比对。动态编程的思想体现于此。最终的最优局部比对结果是由逐步扩展当前位置最优比对结果来实现的（即通过作出一系列的对如何比对当前位置可以得到更高分数的决定来得到最终结果）。矩阵的长度和宽度分别为两序列长度+1。额外的首行和首列是为了让序列的末端可以和空位匹配。首行和首列均设为0。以CTCAA和GGTCA为例，初始得分矩阵为：以序列TGTTACGG和GGTTGACTA为例。使用如下参数：初始化得分矩阵，并进行打分。如下图所示。该图展示了该得分矩阵前三个元素的打分过程。黄色代表当前正在考虑的两个碱基。红色代表该位置得分的来源。打分完成后的得分矩阵如左下图所示。其中蓝色代表最高分元素。注意某些元素有不止一种得分来源，这样在回溯过程中会产生多个路径。如果有多个相同的最高分，则应当考虑每个最高分的回溯。回溯过程见右下图。局部最相似序列在从最高分回溯的过程中自右向左产生。该例的结果为：