序列比对那点事儿
本来这应该是一本书，那样的话的确需要花一点心思，就写成一篇短文吧。
从字符比对开始说起吧。第一个问题最简单，如何判断两个字符串是相等的。
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
",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>
" 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
";
print "$align2
";
|
下面就是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 << "
Dynamic programming matrix: " << "
";
print_matrix( F, seq_1, seq_2 );
cout << "
Traceback matrix: " << "
";
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 << "
";
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 << "
";
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>
" 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
";
print "$align2
";
|
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等是如何针对自己的问题而提高比对效率的。