Needleman-Wunsch 算法perl实现版

第一次写博客心情还是有点小激动的。其实写博客的初衷是为了将自己所学到的知识分享给大家，在这个过程中，也能让自己对整个过程有着更清晰的认识，锻炼自己的分析和归纳能力。

废话不多说，今天第一弹为大家献上生物信息学中最基本也是非常关键的序列比对算法。序列比对分为两序列比对和多序列比对，在两序列比对中主要存在两类算法：穷举式精确算法和基于启发式的近似算法。今天给大家介绍的就是穷举式精确算法中的Needleman-Wunsch算法，此算法是全局比对算法，其要点在于要想得到两个序列的最佳比对，首先求得其子序列的最佳比对，子序列的最佳比对又接着从子序列的子序列的最佳比对中取得，这样一层一层递推下去，最后得到两序列比对的最佳形式。

假定采用“固定空位罚分模型”，打分模板采用S(a,a)=1,S(a,b)=-1,S(-,a)=S(a,-)=-2;两条序列分别为s和t，对应长度为m和n。算法流程如下：

1、构建一个（m+1）×（n+1）的矩阵H，其中每个值H_i,j的含义为序列s中前i个子序列和t中前j个子序列的最佳比对得分。

2、H_i,j的最佳比对得分存在三种状况：

• s_i和t_j刚好比对或替换----子序列为H_i-1,j-1;
• s_i和空位比对------------子序列为H_i-1,j;
• 空位和t_j比对------------子序列为H_i,j-1;

以这种方式逐层往下递推，同时每个H_i,j都需记录下其取得最佳比对是从何种子序列而来。

3、从H_m,n中寻找一条到达H_0,0的路径，此路径即为两序列比对的最佳形式

若想取得所有最佳比对的序列形式，则问题转化为在两个节点之间寻找所有路径，即有向图的遍历问题。对于两序列比对，这样做或许意义不大，但对于自己而言刚好是可以锻炼的一个机会，可大学时学的数据结构早已还给老师，于是在网上复习回顾了一下有向图的遍历算法，决定采用深度优先遍历，取得所有最佳比对的路径。

下面是用perl实现的算法：

use strict;
use warnings;

my @F=([],[]);            ##存放两条比对序列的二维数组；
my @H;                    ##存放比对序列得分的二维数组；
my @S;                    ##存放二维数组H中每个值状态的三维数组；
my @S_number;             ##存放H[i][j]所有可选状态的数目，与三维数组S相对应；
my @S_number_copy;        ##存放S_number的一个copy,用于判断回溯节点的原始可选择路径数目；
my @route;                ##存放从H[m][n]到H[0][0]的路径的边；
my @route_node;           ##存放从H[m][n]到H[0][0]的路径的节点；
my @route_node_number;    ##存放从H[m][n]到H[0][0]的路径的节点的剩余状态数；
my $route_number=0;       ##所有最佳比对的数目；
my $m;my $n;              ##两条序列的长度；

my @query;                ##最终输出序列的形式；
my @target;               ##最终输出序列的形式；
my $t_q=0;                
my $t_t=0;

##将两条待比对序列存进二维数组F；
open QUERY,'<','E:\perl\Algorithm\query.txt' or die "$!";
while(<QUERY>){
    chomp;
    s/\s+//g;
    s/^\s+//g;
    s/(\w+)/\U$1/g;
    my @first=split //;
    push @{$F[0]},@first;
}
close QUERY;
open TARGET,'<','E:\perl\Algorithm\target.txt' or die "$!";
while (<TARGET>){
    chomp;
    s/\s+//g;
    s/^\s+//g;
    s/(\w+)/\U$1/g;
    my @second=split //;
    push @{$F[1]},@second;
}
close TARGET;

##二维数组H和S_number的初始化；
$m=@{$F[0]};
$n=@{$F[1]};
$H[0][0]=0;
$S_number[0][0]=0;
@{$H[0]}=map {$_*(-2)} my @n=(0..$n);
for(my $i=1;$i<=$n;$i++){
    $S_number[0][$i]=1;
}
for(my $i=1;$i<=$m;$i++){
    $H[$i][0]=(-2)*$i;
    $S_number[$i][0]=1;
}
##三维数组S的初始化
@{$S[0][0]}=();
for(my $i=1;$i<=$n;$i++){
    @{$S[0][$i]}='left';
}
for (my $i=1;$i<=$m;$i++){
    @{$S[$i][0]}='up';
}
##计算H[i][j]子程序；
sub max{
    my %status;     
    my $a=$_[0];my $b=$_[1];
    if($F[0][$a-1] eq $F[1][$b-1]){
        $status{left_up}=$H[$a-1][$b-1]+1;
    }
    else{
        $status{left_up}=$H[$a-1][$b-1]-1;
    }
    $status{up}=$H[$a-1][$b]-2;
    $status{left}=$H[$a][$b-1]-2;
    my $max=$status{left_up};
    my @score_status;
##利用高水线，取得三种状态下的最大值；
    foreach (keys %status){
        if($status{$_} > $max){
            $max=$status{$_};
        }
    }
    push @score_status,$max;
    foreach (keys %status){
        if($status{$_} == $max){
            push @score_status,$_;
        }
    }
    @score_status;
}
##计算每个H[i][j]的值并将每个值对应的状态存入三维数组S中；
for(my $i=1;$i<=$m;$i++){
    for(my $j=1;$j<=$n;$j++){
        my @temp=&max($i,$j);
        $H[$i][$j]=shift @temp;
        my $status_number=@temp;
        $S_number[$i][$j]=$status_number;
        @{$S[$i][$j]}=@temp;
    }
}
##将二维数组复制一份copy为S_number_copy;
for(my $i=0;$i<=$m;$i++){
    for(my $j=0;$j<=$n;$j++){
        $S_number_copy[$i][$j]=$S_number[$i][$j];
    }
}
##输出每个H[i][j]的值以及对应的状态；
for (my $i=0;$i<=$m;$i++){
    print "@{$H[$i]}\n";
}
for(my $i=0;$i<=$m;$i++){
    for(my $j=0;$j<=$n;$j++){
        print "{"."@{$S[$i][$j]}"."}";
    }
    print "\n";
}
##依据三维数组S，寻找路径上的下一个位置
sub next{
    my $status=$_[0];
    my @position;
        if ($status eq 'left'){
        my $i=$_[1];my $j=$_[2]-1;
        @position=($i,$j);
    }
        if ($status eq 'up'){
        my $i=$_[1]-1;my $j=$_[2];
        @position=($i,$j);
    }
        if ($status eq 'left_up'){
        my $i=$_[1]-1;my $j=$_[2]-1;
        @position=($i,$j);        
    }
    return @position;
}

##深度优先遍历路径，每个路径存放入二维数组route中；
ROUTE:{
    while ($m>0 || $n>0){
     my $status;
     my $status_number;
     my @next_position;
     $status=$S[$m][$n][0];
     my $temp=shift @{$S[$m][$n]};
     push @{$S[$m][$n]},$temp;
     $S_number[$m][$n]=$S_number[$m][$n]-1;
     $status_number=$S_number[$m][$n];
     push @{$route[$route_number]},$status;
     push @{$route_node[$route_number]},($m,$n);
     push @{$route_node_number[$route_number]},$status_number;
     @next_position=&next($status,$m,$n);
     $m=$next_position[0];$n=$next_position[1];
    }
    $route_number++;
    @{$route[$route_number-1]}=reverse @{$route[$route_number-1]};
    @{$route_node_number[$route_number-1]}=reverse @{$route_node_number[$route_number-1]};
    @{$route_node[$route_number-1]}=reverse @{$route_node[$route_number-1]};
    my $count=0;
    foreach(@{$route_node_number[$route_number-1]}){
         $m=$route_node[$route_number-1][$count*2+1];$n=$route_node[$route_number-1][$count*2];     
         my $number=@{$route[$route_number-1]};
        my $last=$number-1;                                                                        
         if ($_>0){
        ##若回溯的节点是顶点
            if ($count==$last){
                redo ROUTE;
            }else{
                my @temp1=@{$route[$route_number-1]}[$count+1..$last];
                push @{$route[$route_number]},reverse @temp1;
                
                my @temp2=@{$route_node[$route_number-1]}[$count*2+2..$last*2+1];
                push @{$route_node[$route_number]},reverse @temp2;
                
                my @temp3=@{$route_node_number[$route_number-1]}[$count+1..$last];
                push @{$route_node_number[$route_number]},reverse @temp3;
                
                redo ROUTE;
            }
         }
         ##若回溯的的节点之前已被回溯且可选择路径布置一条，重新将其可回溯状态设为原始可选择路径数目的n-1次；
         if($_<0 && $S_number_copy[$m][$n]>1){
                 my @temp1=@{$route[$route_number-1]}[$count+1..$last];
                push @{$route[$route_number]},reverse @temp1;
                
                my @temp2=@{$route_node[$route_number-1]}[$count*2+2..$last*2+1];
                push @{$route_node[$route_number]},reverse @temp2;
                
                my @temp3=@{$route_node_number[$route_number-1]}[$count+1..$last];
                push @{$route_node_number[$route_number]},reverse @temp3;
                
                $S_number[$m][$n]=$S_number_copy[$m][$n]-1;
                redo ROUTE;
         }
        $count++;
    }
}

##输出两个序列的比对
$route_number=@route;
print "The total route number is:$route_number\n";
for(my $i=0;$i<$route_number;$i++){
    print "@{$route[$i]}\n";
    foreach(@{$route[$i]}){
        if ($_ eq 'left_up'){
            push @query,$F[0][$t_q];
            push @target,$F[1][$t_t];
            $t_q++;
            $t_t++;
            
    }
        if ($_ eq 'left'){
            push @query,"-";
            push @target,$F[1][$t_t];
            $t_t++;
    }
        if ($_ eq 'up'){
            push @query,$F[0][$t_q];
            push @target,"-";
            $t_q++;
    }
 }
    print "@query\n";
    print "@target\n";
    $t_q=0;$t_t=0;
    @query=();@target=();
}

下面是两条测试的序列：

s=CTGAAAA;t=CAGAA;

输出结果如下：

但是大家千万不要用这个程序比对两条较长的序列，一是因为算法本身的空间和时间复杂度都很高，二是因为深度优先遍历的缘故，三是程序本身并没有进行优化。比对较长序列时肯定会导致电脑死机（除非是大型计算机）。如何改进这些问题有待进一步的学习.........

好，今天就到此为止，接下来将为大家献上更多用perl实现的生物信息学算法。