Li_阴宅

这个屌丝很懒，什么也没留下！

热门标签

序列比对那点事儿_构建载体序列比对中间多了一个碱基

作者：Li_阴宅 | 2024-08-20 13:24:00

踩

构建载体序列比对中间多了一个碱基

从字符比对开始说起吧。第一个问题最简单，如何判断两个字符串是相等的。

int strcmp(const char *s1, const char *s2)
{
  int ret = 0;
  while (!(ret = *(unsigned char *) s1 - *(unsigned char *) s2) && *s2) ++s1, ++s2;
 
  if (ret < 0)
    ret = -1;
  else if (ret > 0)
    ret = 1 ;
 
  return ret;
}

这种比较方式是逐一比较字符，如果遇到某个字符不一样了，就停下来，返回结果，否则就认为两个字符串是一样的。

第二个问题就是如果两个字符串一样长，它们之间有几个字符不一样？解答这个问题，有一个专门的名词，叫做汉明距离(hamming distance)。举个例子来说，
agctagct 与 acctaggt 的汉明距离就是2，因为它们相应的位置上有两个字符不一样。

/*
 *  hammimgDistance.c
 *  
 *
 *  Created by 糗世界 on 4/11/12.
 *  Copyright 2012 qiuworld.com. All rights reserved.
 *
 */
 
#include <stdio.h>
 
int hamdist(const char *s1, const char *s2)
{
	int dist = 0;
 
	while(*s1 && *s2)
	{
		if (!(*(unsigned char *)s1==*(unsigned char *)s2)) {
			++dist;
		}
		s1++,s2++;
	}
 
	return dist;
}
 
int main(int argc, const char* argv[]){
	const char *s1=argv[1];
	const char *s2=argv[2];
	printf("Hamming Distance of '%s' and '%s' is %d\n",s1,s2,hamdist(s1,s2));
	return 0;
}

第三个问题，两个字符串，其中一个要经过多少次编辑才能变成另一个？这里的编辑包括变换，删除以及插入。比如kitten要变成sitten就需要把首字符k变成s，就是说需要编辑一次。把sittin变成sitting就需要在末尾加上一个字符g，也是需要编辑一次。这个问题很早就被一个俄罗斯人解决了，他的名字叫Vladimir Levenshtein，所以这个算法就被称为levenshtein距离。比较汉明距离和levenshtein距离，汉明距离只允许替换，而levenshtein距离还允许删除和插入。

levenshtein的具体算法是：

获取字符串Ａ的长度为n，获取字符串Ｂ的长度为m，如果m=0，返回n；如果n=0，返回m。
初始化一个二维数组D[m+1,n+1]，将第一行填充为0至n，将第一列填充为0至m。注意这里的填充不是0。
循环比对A的每一个字符与B的每一个字符。对于字符A[i](i取值1至n)与B[j](j取值1至m)，如果A[i]==B[j]，那么cost为0，否则cost为1。这时比较二维数组D中的三个值：D[i-1,j]+1（左侧值＋1）, D[i,j-1]+1（顶部值＋1）, d[i-1,j-1]+cost（左上角＋cost）；这三个值分别代表删除，插入或者替换。取三个值中的最小值填充D[i,j]
返回D数组最右下角D[m＋1,n＋1]的值就是levenshtein距离。

理解上面的步骤似乎有点难度的样子，我们看一眼一个实例。假设我们需要计算两个字符串atctg及agct的levenshtein距离。第一步，获得两个字符串的长度，分别为5和4.第二步，初始化一个二维数组D[5,6]，将第一行填充为0至5，将第一列填充为0至4。

 	 	a	t	c	t	g
	0	1	2	3	4	5
a	1
g	2
c	3
t	4

第三步，循环，逐列填充。填充第2列。我们看第二行第二列，因为a==a，所以cost=0,左上角值＋cost＝0,左侧值＋1＝2，顶部值＋1＝2。三者中最小值是左上角值＋cost。所以第二行第二列填0。相同步骤，依次向下填充。

 	 	a	t	c	t	g
	0	1	2	3	4	5
a	1	0
g	2	1
c	3	2
t	4	3

填充第3列。

 	 	a	t	c	t	g
	0	1	2	3	4	5
a	1	0	1
g	2	1	1
c	3	2	2
t	4	3	2

填充第4列。

 	 	a	t	c	t	g
	0	1	2	3	4	5
a	1	0	1	2
g	2	1	1	2
c	3	2	2	1
t	4	3	2	2

填充第5列。

 	 	a	t	c	t	g
	0	1	2	3	4	5
a	1	0	1	2	3
g	2	1	1	2	3
c	3	2	2	1	2
t	4	3	2	2	1

填充第6列。

 	 	a	t	c	t	g
	0	1	2	3	4	5
a	1	0	1	2	3	4
g	2	1	1	2	3	3
c	3	2	2	1	2	3
t	4	3	2	2	1	2

得到最右下角的值为2，于是返回levenshtein距离为2.就是说，如果想把agct变换成atctg，需要2步，1，把g换成t，在最后加上g。

/*
 *  levenshteinDistance.cpp
 *  
 *
 *  Created by 糗世界 on 4/11/12.
 *  Copyright 2012 qiuworld.com All rights reserved.
 *
 */
 
#include <iostream>
#include <algorithm>
#include <string>
using namespace std;
 
int minimum(int a, int b, int c){
	int m = min(a,b);
	return min(m,c);
}
 
int levedist(string s1, string s2)
{
	const int n=s1.length()+1;
	const int m=s2.length()+1;
	if(n == 1) return m-1;
	if(m == 1) return n-1;
 
	int **mat= new int*[m];
	for (int j=0; j<m; j++) {
		mat[j] = new int[n];
	}
 
	for (int i =0 ; i < m; ++i) {
		for (int j=0; j < n; ++j) {
			mat[i][j] = 0;
			mat[0][j] = j;
		}
		mat[i][0] = i;
	}
 
	for (int i=1; i<m; ++i) {
		for (int j=1; j<n; ++j) {
			int cost = (s2[i-1]==s1[j-1]) ? 0 : 1;
			mat[i][j] = minimum(mat[i-1][j]+1, mat[i][j-1]+1,mat[i-1][j-1]+cost);
		}
	}
	return mat[m-1][n-1];
}
 
int main(int argc, const char* argv[]){
	string s1=argv[1];
	string s2=argv[2];
	cout << "Levenshtein Distance of '" << s1 << "' and '" << s2 << "' is "<< levedist(s1,s2) << endl;
	return 0;
}

第四个问题，问题3中的编辑距离是不是字符串Ａ变成字符串Ｂ的最短编辑次数？如果说只包含替换，插入和删除三种操作的话，的确是的。在1964年，Frederick J. Damerau在Vladimir I. Levenshtein的编辑距离的基础上提出了一个新的概念，就是交换相邻两字符的顺序算一步操作。这个概念的提出是因为人们在拼写字母为基础的文字时，有可能把字母的顺序搞错，在计算机输入方面，这个问题变得尤其明显。于是得到了levenshtein distance的一种变型Damerau-Levenshtein distance。

举个例子，我们希望把CA变成ABC。如果使用levenshtein距离的话是3，CA->A->AB->ABC。如果使用Damerau-Levenshtein距离的话就只有两步CA->AC->ABC。很明显，如果可以只交换相邻两字符的话，Damerau-Levenshtein距离是很有用的。

我们先来看一下算法。Damerau-Levenshtein距离算法基本和Levenshtein距离算法一致，只是多了一步，在mat[i][j] = minimum(mat[i-1][j]+1, mat[i][j-1]+1,mat[i-1][j-1]+cost);判断完之后，再判断一下s1[i]==s2[j-1]&&s1[i-1]＝＝s2[j]，如果是真的话，那么就是相邻交换，mat[i][j]就可以从mat[i-1][j]+1, mat[i][j-1]+1,mat[i-1][j-1]+cost,mat[i-2][j-2]+cost四个值当中选择最小值。下面的代码与这里写的会不一样，是因为它对问题做了一下优化，算得更快些，但其思路确实是开始于这里介绍的。

/*
 *  DamerauLevenshteinDistance.cpp
 *  
 *
 *  Created by 糗世界 on 4/11/12.
 *  Copyright 2012 qiuworld.com. All rights reserved.
 *
 */
 
#include <iostream>
#include <algorithm>
#include <map>
#include <string>
using namespace std;
 
int damelevedist(string s1, string s2)
{
	const int n=s1.length();
	const int m=s2.length();
	if(n == 0) return m;
	if(m == 0) return n;
	int INF=m+n;
	string s3=s1+s2;
 
	int **mat= new int*[m+2];
	for (int j=0; j<=m+1; j++) {
		mat[j] = new int[n+2];
	}
 
	mat[0][0] = INF;
 
	for (int i=0; i<=m; i++) {
		mat[i+1][1] = i;
		mat[i+1][0] = INF;
	}
	for (int j=0; j<=n; j++) {
		mat[1][j+1] = j;
		mat[0][j+1] = INF;
	}
 
	map<char,int> sd;
	for (int key=0; key < INF; key++) {
		char letter = s3[key];
		sd[letter] = 0;
	}
 
	for (int i=1; i<=m; i++) {
		int db = 0;
		for (int j=1; j<=n; j++) {
			char i1 = sd[s2[i-1]];
			int j1 = db;
 
			if (s2[i-1]==s1[j-1]) {
				mat[i+1][j+1] = mat[i][j];
				db = j;
			}else {
				mat[i+1][j+1] = min(mat[i][j],min(mat[i+1][j],mat[i][j+1])) + 1;
			}
 
			mat[i+1][j+1] = min(mat[i+1][j+1], mat[i1][j1] + (i-i1-1) + 1 + (j -j1 -1));
		}
		sd[s1[i-1]] = i;
	}
 
	return mat[m+1][n+1];
}
 
int main(int argc, const char* argv[]){
	string s1=argv[1];
	string s2=argv[2];
	cout << "Damerau-Levenshtein Distance of '" << s1 << "' and '" << s2 << "' is "<< damelevedist(s1,s2) << endl;
	return 0;
}

我们来比较一下levenshtein distance及damerau-levenshtein分数距阵的差别。比如我们想比较的两个字符串分别为”acttg”和”atctg”。

qiuworld.com$ ./damelevedist acttg atctg
 	 	a	c	t	t	g
 	10	10	10	10	10	10	
a	10	0	1	2	3	4	
t	10	1	0	1	2	3	
c	10	2	1	1	1	2	
t	10	3	2	1	1	2	
g	10	4	3	2	1	1	
Damerau-Levenshtein Distance of 'acttg' and 'atctg' is 1
qiuworld.com$ ./levedist acttg atctg
 	 	a	c	t	t	g
 	0	1	2	3	4	5	
a	1	0	1	2	3	4	
t	2	1	1	1	2	3	
c	3	2	1	2	2	3	
t	4	3	2	1	2	3	
g	5	4	3	2	2	2	
Levenshtein Distance of 'acttg' and 'atctg' is 2

好了，前面讲了那么多字符比对的事情，现在应该开始回到正题了，关于序列比对。首先我们来完成一个任务，手工比对一下ATGGCGT和ATGAGT

ATGGCGT
*** !**
ATG-AGT

是不是很容易？但是我们是如何得出这样的结果的呢？可能的思路应该是从第一个字母开始，查找两个序列一致的部分，然后分割成小段，把一致的部分都对齐，然后把不一致的部分补齐。
但如果需要使用计算机来完成这一过程，似乎有些不太合理，因为计算量的问题。好吧，简单的，计算机是如何来解决这一问题的呢？

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

两序列比对的结果必须遵守它们原有的序列顺序。
可以引入空格。空格的意思就是在其中一个序列当中插入一个空的字符，表示该位置相对的另一序列的碱基对应一个空碱基。
空格对空格是没有意义的，也是不被允许的。

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

A T G G C G T
A T G - A G T
1 1 1 0-1 1 1
score: +1+1+1+0-1+1+1=4
A T G G C G T
A - T G A G T
1 0-1 1-1 1 1
score: +1+0-1+1-1+1+1=2

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

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

	C	T	A	G
C 	 1	-1	-1	-1
T	-1	 1	-1	-1
A	-1	-1	 1	-1
G	-1	-1	-1	 1

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

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

	C	T	A	G
C 	 2	 1	-1	-1
T	 1	 2	-1	-1
A	-1	-1	 2	 1
G	-1	-1	 1	 2

当然还有更复杂的，就是依据统计结果用马尔可夫链来生成评分表，这个马尔可夫链可以是一维的，就是说我们只考虑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。

traceback示例

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

D
D

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

ND
ND

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

END
-ND

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

SEND
A-ND

再向左上角得到done，完成比对的traceback。
下面就是Trace-back with Needleman-Wunsch algorithm的perl程序示例：

# Needleman-Wunsch  Algorithm 
 
# usage statement
die "usage: $0 <sequence 1> <sequence 2>\n" unless @ARGV == 2;
 
# get sequences from command line
my ($seq1, $seq2) = @ARGV;
 
# scoring scheme
my $MATCH    =  1; # +1 for letters that match
my $MISMATCH = -1; # -1 for letters that mismatch
my $GAP      = -1; # -1 for any gap
 
# initialization
my @matrix;
$matrix[0][0]{score}   = 0;
$matrix[0][0]{pointer} = "none";
for(my $j = 1; $j <= length($seq1); $j++) {
    $matrix[0][$j]{score}   = $GAP * $j;
    $matrix[0][$j]{pointer} = "left";
}
for (my $i = 1; $i <= length($seq2); $i++) {
    $matrix[$i][0]{score}   = $GAP * $i;
    $matrix[$i][0]{pointer} = "up";
}
 
# fill
for(my $i = 1; $i <= length($seq2); $i++) {
    for(my $j = 1; $j <= length($seq1); $j++) {
        my ($diagonal_score, $left_score, $up_score);
 
        # calculate match score
        my $letter1 = substr($seq1, $j-1, 1);
        my $letter2 = substr($seq2, $i-1, 1);                            
        if ($letter1 eq $letter2) {
            $diagonal_score = $matrix[$i-1][$j-1]{score} + $MATCH;
        }
        else {
            $diagonal_score = $matrix[$i-1][$j-1]{score} + $MISMATCH;
        }
 
        # calculate gap scores
        $up_score   = $matrix[$i-1][$j]{score} + $GAP;
        $left_score = $matrix[$i][$j-1]{score} + $GAP;
 
        # choose best score
        if ($diagonal_score >= $up_score) {
            if ($diagonal_score >= $left_score) {
                $matrix[$i][$j]{score}   = $diagonal_score;
                $matrix[$i][$j]{pointer} = "diagonal";
            }
        else {
                $matrix[$i][$j]{score}   = $left_score;
                $matrix[$i][$j]{pointer} = "left";
            }
        } else {
            if ($up_score >= $left_score) {
                $matrix[$i][$j]{score}   = $up_score;
                $matrix[$i][$j]{pointer} = "up";
            }
            else {
                $matrix[$i][$j]{score}   = $left_score;
                $matrix[$i][$j]{pointer} = "left";
            }
        }
    }
}
 
# trace-back
 
my $align1 = "";
my $align2 = "";
 
# start at last cell of matrix
my $j = length($seq1);
my $i = length($seq2);
 
while (1) {
    last if $matrix[$i][$j]{pointer} eq "none"; # ends at first cell of matrix
 
    if ($matrix[$i][$j]{pointer} eq "diagonal") {
        $align1 .= substr($seq1, $j-1, 1);
        $align2 .= substr($seq2, $i-1, 1);
        $i--;
        $j--;
    }
    elsif ($matrix[$i][$j]{pointer} eq "left") {
        $align1 .= substr($seq1, $j-1, 1);
        $align2 .= "-";
        $j--;
    }
    elsif ($matrix[$i][$j]{pointer} eq "up") {
        $align1 .= "-";
        $align2 .= substr($seq2, $i-1, 1);
        $i--;
    }    
}
 
$align1 = reverse $align1;
$align2 = reverse $align2;
print "$align1\n";
print "$align2\n";

下面就是Trace-back with Needleman-Wunsch algorithm的C++程序示例：
nw.h

/*
 * nw.h for program nw.
 *
 */
 
#include <iostream>
#include <string>
#include <algorithm>
#include <stdlib.h>
#include <stdio.h>
 
#define DEBUG 0
 
using namespace std;
 
 
extern int  nw( 
                 string, string, 
                 string&, string&,
                 bool
              );
 
extern int  nw_align( 
                      int **, char **,
                      string, string, 
                      string&, string&,
                      int 
                    );
 
extern void  dpm_init        ( int **, char **, int, int, int );
extern void  print_al        ( string&, string& );
extern void  print_matrix    ( int ** const, string, string );
extern void  print_traceback ( char ** const, string, string );
extern int   max             ( int, int, int, char * );

nw.cpp

/*-------------------------------------------
 * 
 *        nw.c++ for program nw
 * 
 -------------------------------------------*/
 
#include "nw.h"
 
using namespace std;
 
 
int nw(                                                          
        string       seq_1,          /*  Needleman-Wunsch   */
        string       seq_2,          /*  algorithm for      */
        string&      seq_1_al,       /*  global alignment   */
        string&      seq_2_al,       /*  of nt sequence.    */
        bool         prm
      )
{
        int  d = 2 ;                 /* gap penalty */
 
        int  L1 = seq_1.length();
        int  L2 = seq_2.length();
 
        // Dynamic programming matrix
        int ** F = new int * [ L2+1 ];
        for( int i = 0; i <= L2; i++ )  F[ i ] = new int [ L1 ];
 
        // Traceback matrix
        char ** traceback = new char * [ L2+1 ];
        for( int i = 0; i <= L2; i++ )  traceback[ i ] = new char [ L1 ];
 
        // Initialize traceback and F matrix (fill in first row and column)
        dpm_init( F, traceback, L1, L2, d );
 
        // Create alignment
        nw_align( F, traceback, seq_1, seq_2, seq_1_al, seq_2_al, d );
 
        #if DEBUG
            int  L_al = seq_1_al.length();
            cout << "Length after alignment: " << L_al << endl;
        #endif
 
        if( prm )
        {
                cout << "\nDynamic programming matrix: " << "\n\n";
                print_matrix( F, seq_1, seq_2 );
 
                cout << "\nTraceback matrix: " << "\n\n";
                print_traceback( traceback, seq_1, seq_2 );
 
                cout << endl;
        }
 
        for( int i = 0; i <= L2; i++ )  delete F[ i ];
        delete[] F;
        for( int i = 0; i <= L2; i++ )  delete traceback[ i ];
        delete[] traceback;
 
        return  0 ;
}
 
 
void  dpm_init( int ** F, char ** traceback, int L1, int L2, int d )
{
        F[ 0 ][ 0 ] =  0 ;
        traceback[ 0 ][ 0 ] = 'n' ;
 
        int i=0, j=0;
 
        for( j = 1; j <= L1; j++ )
        {
                F[ 0 ][ j ] =  -j * d ;
                traceback[ 0 ][ j ] =  '-' ;
        }
        for( i = 1; i <= L2; i++ )
        {
                F[ i ][ 0 ] =  -i * d ;
                traceback[ i ][ 0 ] =  '|' ;
        }
}
 
 
int nw_align(                  // Needleman-Wunsch algorithm
              int **     F,
              char **    traceback,
              string     seq_1,
              string     seq_2,
              string&    seq_1_al,
              string&    seq_2_al,
              int        d         // Gap penalty
            )
{
        int        k = 0, x = 0, y = 0;
        int        fU, fD, fL ;
        char       ptr, nuc ;
        int        i = 0, j = 0;
 
        const int  a =  2;   // Match
        const int  b = -1;   // Mismatch
 
        const int  s[ 4 ][ 4 ] = { { a, b, b, b },    /* substitution matrix */
                                   { b, a, b, b },
                                   { b, b, a, b },
                                   { b, b, b, a } } ;
 
        int  L1 = seq_1.length();
        int  L2 = seq_2.length();
 
        for( i = 1; i <= L2; i++ )
        {
                for( j = 1; j <= L1; j++ )
                {
                        nuc = seq_1[ j-1 ] ;
 
                        switch( nuc )
                        {
                                case 'A':  x = 0 ;  break ;
                                case 'C':  x = 1 ;  break ;
                                case 'G':  x = 2 ;  break ;
                                case 'T':  x = 3 ;
                        }
 
                        nuc = seq_2[ i-1 ] ;
 
                        switch( nuc )
                        {
                                case 'A':  y = 0 ;  break ;
                                case 'C':  y = 1 ;  break ;
                                case 'G':  y = 2 ;  break ;
                                case 'T':  y = 3 ;
                        }
 
                        fU = F[ i-1 ][ j ] - d ;
                        fD = F[ i-1 ][ j-1 ] + s[ x ][ y ] ;
                        fL = F[ i ][ j-1 ] - d ;
 
                        F[ i ][ j ] = max( fU, fD, fL, &ptr ) ;
 
                        traceback[ i ][ j ] =  ptr ;
                }
        }
        i-- ; j-- ;
 
        while( i > 0 || j > 0 )
        {
                switch( traceback[ i ][ j ] )
                {
                        case '|' :      seq_1_al += '-' ; 
                                        seq_2_al += seq_2[ i-1 ] ; 
                                        i-- ;
                                        break ;
 
                        case '\\':      seq_1_al += seq_1[ j-1 ] ; 
                                        seq_2_al += seq_2[ i-1 ] ; 
                                        i-- ;  j-- ;
                                        break ;
 
                        case '-' :      seq_1_al += seq_1[ j-1 ] ; 
                                        seq_2_al += '-' ; 
                                        j-- ;
                }
                k++ ;
        }
 
        reverse( seq_1_al.begin(), seq_1_al.end() );
        reverse( seq_2_al.begin(), seq_2_al.end() );
 
        return  0 ;
}
 
 
int  max( int f1, int f2, int f3, char * ptr )
{
        int  max = 0 ;
 
        if( f1 >= f2 && f1 >= f3 )  
        {
                max = f1 ;
                *ptr = '|' ;
        }
        else if( f2 > f3 )              
        {
                max = f2 ;
                *ptr = '\\' ;
        }
        else
        {
                max = f3 ;
                *ptr = '-' ;
        }
 
        return  max ;   
}
 
 
void  print_matrix( int ** F, string seq_1, string seq_2 )
{
        int  L1 = seq_1.length();
        int  L2 = seq_2.length();
 
        cout << "        ";
        for( int j = 0; j < L1; j++ )
        {
                cout << seq_1[ j ] << "   ";
        }
        cout << "\n  ";
 
        for( int i = 0; i <= L2; i++ )
        {
                if( i > 0 )
                {
                        cout << seq_2[ i-1 ] << " ";
                }
                for( int j = 0; j <= L1; j++ )
                {
                        cout.width( 3 );
                        cout << F[ i ][ j ] << " ";
                }
                cout << endl;
        }
}
 
 
void  print_traceback( char ** traceback, string seq_1, string seq_2 )
{
        int  L1 = seq_1.length();
        int  L2 = seq_2.length();
 
        cout << "    ";
        for( int j = 0; j < L1; j++ )
        {
                cout << seq_1[ j ] << " ";
        }
        cout << "\n  ";
 
        for( int i = 0; i <= L2; i++ )
        {
                if( i > 0 )
                {
                        cout << seq_2[ i-1 ] << " ";
                }
                for( int j = 0; j <= L1; j++ )
                {
                        cout << traceback[ i ][ j ] << " ";
                }
                cout << endl;
        }
}
 
 
void  print_al( string& seq_1_al, string& seq_2_al )
{
        cout << seq_1_al << endl;
        cout << seq_2_al << endl;
}

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

矩阵最初始的两条边（最顶部及最左侧的值）全部以零填充。
之前讲的三个值中取最大值，值不能低于零。否则不计入得分表。
回溯时选择最大值，直到遇到0为止。

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

空格：－1
相同：＋2
不同：－1

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

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

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

我们来看一下SW算法的perl代码吧。

# Smith-Waterman  Algorithm
 
# usage statement
die "usage: $0 <sequence 1> <sequence 2>\n" unless @ARGV == 2;
 
# get sequences from command line
my ($seq1, $seq2) = @ARGV;
 
# scoring scheme
my $MATCH    =  1; # +1 for letters that match
my $MISMATCH = -1; # -1 for letters that mismatch
my $GAP      = -1; # -1 for any gap
 
# initialization
my @matrix;
$matrix[0][0]{score}   = 0;
$matrix[0][0]{pointer} = "none";
for(my $j = 1; $j <= length($seq1); $j++) {
    $matrix[0][$j]{score}   = 0;
    $matrix[0][$j]{pointer} = "none";
}
for (my $i = 1; $i <= length($seq2); $i++) {
    $matrix[$i][0]{score}   = 0;
    $matrix[$i][0]{pointer} = "none";
}
 
# fill
my $max_i     = 0;
my $max_j     = 0;
my $max_score = 0;
 
for(my $i = 1; $i <= length($seq2); $i++) {
    for(my $j = 1; $j <= length($seq1); $j++) {
        my ($diagonal_score, $left_score, $up_score);
 
        # calculate match score
        my $letter1 = substr($seq1, $j-1, 1);
        my $letter2 = substr($seq2, $i-1, 1);       
        if ($letter1 eq $letter2) {
            $diagonal_score = $matrix[$i-1][$j-1]{score} + $MATCH;
        }
        else {
            $diagonal_score = $matrix[$i-1][$j-1]{score} + $MISMATCH;
        }
 
        # calculate gap scores
        $up_score   = $matrix[$i-1][$j]{score} + $GAP;
        $left_score = $matrix[$i][$j-1]{score} + $GAP;
 
        if ($diagonal_score <= 0 and $up_score <= 0 and $left_score <= 0) {
            $matrix[$i][$j]{score}   = 0;
            $matrix[$i][$j]{pointer} = "none";
            next; # terminate this iteration of the loop
        }
 
        # choose best score
        if ($diagonal_score >= $up_score) {
            if ($diagonal_score >= $left_score) {
                $matrix[$i][$j]{score}   = $diagonal_score;
                $matrix[$i][$j]{pointer} = "diagonal";
            }
            else {
                $matrix[$i][$j]{score}   = $left_score;
                $matrix[$i][$j]{pointer} = "left";
            }
        } else {
            if ($up_score >= $left_score) {
                $matrix[$i][$j]{score}   = $up_score;
                $matrix[$i][$j]{pointer} = "up";
            }
            else {
                $matrix[$i][$j]{score}   = $left_score;
                $matrix[$i][$j]{pointer} = "left";
            }
        }
 
        # set maximum score
        if ($matrix[$i][$j]{score} > $max_score) {
            $max_i     = $i;
            $max_j     = $j;
            $max_score = $matrix[$i][$j]{score};
        }
    }
}
 
# trace-back
 
my $align1 = "";
my $align2 = "";
 
my $j = $max_j;
my $i = $max_i;
 
while (1) {
    last if $matrix[$i][$j]{pointer} eq "none";
 
    if ($matrix[$i][$j]{pointer} eq "diagonal") {
        $align1 .= substr($seq1, $j-1, 1);
        $align2 .= substr($seq2, $i-1, 1);
        $i--; $j--;
    }
    elsif ($matrix[$i][$j]{pointer} eq "left") {
        $align1 .= substr($seq1, $j-1, 1);
        $align2 .= "-";
        $j--;
    }
    elsif ($matrix[$i][$j]{pointer} eq "up") {
        $align1 .= "-";
        $align2 .= substr($seq2, $i-1, 1);
        $i--;
    }   
}
 
$align1 = reverse $align1;
$align2 = reverse $align2;
print "$align1\n";
print "$align2\n";

c++代码

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <string>
#include <cmath>
#include <sys/time.h>
using namespace std;
 
double similarity_score(char a,char b);
double find_array_max(double array[],int length);
void insert_at(char arr[], int n, int idx, char val);
void checkfile(int open, char filename[]);
string read_sequence(ifstream& f);
 
int ind;
double mu,delta;
 
 
int main(int argc, char** argv){
 
  // read info from arguments
  if(argc!=6){
    cout<<"Give me the propen number of input arguments:"<<endl<<"1 : mu"<<endl;
    cout<<"2 : delta"<<endl<<"3 : filename sequence A"<<endl<<"4 : filename sequence B"<<endl;
    cout<<"5 : maximal length N of sequences"<<endl;exit(1);
  }
  mu    = atof(argv[1]);
  delta = atof(argv[2]);
  /
  // give it the filenames
  char *nameof_seq_a = argv[3];
  char *nameof_seq_b = argv[4];
  int N_max = atoi(argv[5]);
  string seq_a,seq_b; 
 
  // read the sequences into two vectors:
  ifstream stream_seq_b;                      // first define the input-streams for seq_a and seq_b
  stream_seq_b.open(nameof_seq_b);            // the same for seq_b
  checkfile(! stream_seq_b,nameof_seq_b);
  seq_b = read_sequence(stream_seq_b);
  ifstream stream_seq_a;
  stream_seq_a.open(nameof_seq_a);            // open the file for input
  checkfile(! stream_seq_a,nameof_seq_a);     // check, whether the file was opened successfully
  seq_a = read_sequence(stream_seq_a);  
 
 
  // string s_a=seq_a,s_b=seq_b;
  int N_a = seq_a.length();                     // get the actual lengths of the sequences
  int N_b = seq_b.length();
 
  
 
  // initialize H
  double H[N_a+1][N_b+1];     
  for(int i=0;i<=N_a;i++){
    for(int j=0;j<=N_b;j++){
      H[i][j]=0.;
    }
  } 
 
  double temp[4];
  int I_i[N_a+1][N_b+1],I_j[N_a+1][N_b+1];     // Index matrices to remember the 'path' for backtracking
 
  // here comes the actual algorithm
 
  for(int i=1;i<=N_a;i++){
    for(int j=1;j<=N_b;j++){
      temp[0] = H[i-1][j-1]+similarity_score(seq_a[i-1],seq_b[j-1]); 
      temp[1] = H[i-1][j]-delta;                  
      temp[2] = H[i][j-1]-delta;                 
      temp[3] = 0.;
      H[i][j] = find_array_max(temp,4);
      switch(ind){
      case 0:                                  // score in (i,j) stems from a match/mismatch
   	I_i[i][j] = i-1;
	I_j[i][j] = j-1;
	break;
      case 1:                                  // score in (i,j) stems from a deletion in sequence A
     	I_i[i][j] = i-1;
	I_j[i][j] = j;
	break;
      case 2:                                  // score in (i,j) stems from a deletion in sequence B
      	I_i[i][j] = i;
	I_j[i][j] = j-1;
	break;
      case 3:                                  // (i,j) is the beginning of a subsequence
      	I_i[i][j] = i;
	I_j[i][j] = j;	
	break;
      }
    }
  }
 
  // Print the matrix H to the console
  cout<<"**********************************************"<<endl;
  cout<<"The scoring matrix is given by  "<<endl<<endl;
  for(int i=1;i<=N_a;i++){
    for(int j=1;j<=N_b;j++){
      cout<<H[i][j]<<" ";
    }
    cout<<endl;
    }
 
  // search H for the maximal score
  double H_max = 0.;
  int i_max=0,j_max=0;
  for(int i=1;i<=N_a;i++){
    for(int j=1;j<=N_b;j++){
      if(H[i][j]>H_max){
	H_max = H[i][j];
	i_max = i;
	j_max = j;
      }
    }
  }
 
  //cout<<H_max<<endl;
 
     // Backtracking from H_max
  int current_i=i_max,current_j=j_max;
  int next_i=I_i[current_i][current_j];
  int next_j=I_j[current_i][current_j];
  int tick=0;
  char consensus_a[N_a+N_b+2],consensus_b[N_a+N_b+2];
 
  while(((current_i!=next_i) || (current_j!=next_j)) && (next_j!=0) && (next_i!=0)){
 
    if(next_i==current_i)  consensus_a[tick] = '-';                  // deletion in A
    else                   consensus_a[tick] = seq_a[current_i-1];   // match/mismatch in A
 
    if(next_j==current_j)  consensus_b[tick] = '-';                  // deletion in B
    else                   consensus_b[tick] = seq_b[current_j-1];   // match/mismatch in B
 
    current_i = next_i;
    current_j = next_j;
    next_i = I_i[current_i][current_j];
    next_j = I_j[current_i][current_j];
    tick++;
    }
 
 // Output of the consensus motif to the console
  cout<<endl<<"***********************************************"<<endl;
  cout<<"The alignment of the sequences"<<endl<<endl;
  for(int i=0;i<N_a;i++){cout<<seq_a[i];}; cout<<"  and"<<endl;
  for(int i=0;i<N_b;i++){cout<<seq_b[i];}; cout<<endl<<endl;
  cout<<"is for the parameters  mu = "<<mu<<" and delta = "<<delta<<" given by"<<endl<<endl;  
  for(int i=tick-1;i>=0;i--) cout<<consensus_a[i]; 
  cout<<endl;
  for(int j=tick-1;j>=0;j--) cout<<consensus_b[j];
  cout<<endl;
 
 
} // END of main
 
 
 
/
// auxiliary functions used by main:
/
 
 
void checkfile(int open, char filename[]){
 
  if (open){cout << "Error: Can't open the file "<<filename<<endl;exit(1);}
  else cout<<"Opened file "<<filename<<endl;
}
 
/
 
double similarity_score(char a,char b){
 
  double result;
  if(a==b){
      result=1.;
    }
  else{
      result=-mu;
    }
  return result;
}
 
/
 
double find_array_max(double array[],int length){
 
  double max = array[0];            // start with max = first element
  ind = 0;
 
  for(int i = 1; i<length; i++){
      if(array[i] > max){
	max = array[i];
	ind = i; 
      }
  }
  return max;                    // return highest value in array
}
 
/
 
string read_sequence(ifstream& f)
{
  // overflows.
  string seq;
  char line[5000];
  while( f.good() )
    {
      f.getline(line,5000);
      // 		cout << "Line:" << line << endl;
      if( line[0] == 0 || line[0]=='#' )
	continue;
      for(int i = 0; line[i] != 0; ++i)
	{
	  int c = toupper(line[i]);
	  if( c != 'A' && c != 'G' && c != 'C' && c != 'T' )
	    continue;
	  //cout << char(c);
	  seq.push_back(char(c));
	}
    }
  return seq;
}

结语：当我们开始使用不同的算法的时候，我们已经开始脱离生物知识而集中精力于计算机知识了。而在计算方面，就会需要考虑比如算法的复杂性啊，动态编程啊，之类的，那也许是一个从生物出发，走入计算机，再进入数学的过程。如果你想了解更多的算法方面的知识，可从这本文的这些基础知识出发，继续了解blast，blat，bowtie等是如何针对自己的问题而提高比对效率的。

http://blog.qiuworld.com:8080/archives/2637

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Li_阴宅/article/detail/1007396