Needleman-Wunsch 算法perl实现版

jjfjjj

　　第一次写博客心情还是有点小激动的。其实写博客的初衷是为了将自己所学到的知识分享给大家，在这个过程中，也能让自己对整个过程有着更清晰的认识，锻炼自己的分析和归纳能力。
　　废话不多说，今天第一弹为大家献上生物信息学中最基本也是非常关键的序列比对算法。序列比对分为两序列比对和多序列比对，在两序列比对中主要存在两类算法：穷举式精确算法和基于启发式的近似算法。今天给大家介绍的就是穷举式精确算法中的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，其中每个值Hi,j的含义为序列s中前i个子序列和t中前j个子序列的最佳比对得分。
　　2、Hi,j的最佳比对得分存在三种状况：

• si和tj刚好比对或替换----子序列为Hi-1,j-1;
  
• si和空位比对------------子序列为Hi-1,j;
  
• 空位和tj比对------------子序列为Hi,j-1;
  

　　以这种方式逐层往下递推，同时每个Hi,j都需记录下其取得最佳比对是从何种子序列而来。
　　3、从Hm,n中寻找一条到达H0,0的路径，此路径即为两序列比对的最佳形式
　　若想取得所有最佳比对的序列形式，则问题转化为在两个节点之间寻找所有路径，即有向图的遍历问题。对于两序列比对，这样做或许意义不大，但对于自己而言刚好是可以锻炼的一个机会，可大学时学的数据结构早已还给老师，于是在网上复习回顾了一下有向图的遍历算法，决定采用深度优先遍历，取得所有最佳比对的路径。
　　下面是用perl实现的算法：



  1 use strict;
  2 use warnings;
  3
  4 my @F=([],[]);            ##存放两条比对序列的二维数组；
  5 my @H;                    ##存放比对序列得分的二维数组；
  6 my @S;                    ##存放二维数组H中每个值状态的三维数组；
  7 my @S_number;             ##存放H[j]所有可选状态的数目，与三维数组S相对应；
  8 my @S_number_copy;        ##存放S_number的一个copy,用于判断回溯节点的原始可选择路径数目；
  9 my @route;                ##存放从H[m][n]到H[0][0]的路径的边；
10 my @route_node;           ##存放从H[m][n]到H[0][0]的路径的节点；
11 my @route_node_number;    ##存放从H[m][n]到H[0][0]的路径的节点的剩余状态数；
12 my $route_number=0;       ##所有最佳比对的数目；
13 my $m;my $n;              ##两条序列的长度；
14
15 my @query;                ##最终输出序列的形式；
16 my @target;               ##最终输出序列的形式；
17 my $t_q=0;               
18 my $t_t=0;
19
20 ##将两条待比对序列存进二维数组F；
21 open QUERY,'<','E:\perl\Algorithm\query.txt' or die "$!";
22 while(<QUERY>){
23     chomp;
24     s/\s+//g;
25     s/^\s+//g;
26     s/(\w+)/\U$1/g;
27     my @first=split //;
28     push @{$F[0]},@first;
29 }
30 close QUERY;
31 open TARGET,'<','E:\perl\Algorithm\target.txt' or die "$!";
32 while (<TARGET>){
33     chomp;
34     s/\s+//g;
35     s/^\s+//g;
36     s/(\w+)/\U$1/g;
37     my @second=split //;
38     push @{$F[1]},@second;
39 }
40 close TARGET;
41
42 ##二维数组H和S_number的初始化；
43 $m=@{$F[0]};
44 $n=@{$F[1]};
45 $H[0][0]=0;
46 $S_number[0][0]=0;
47 @{$H[0]}=map {$_*(-2)} my @n=(0..$n);
48 for(my $i=1;$i<=$n;$i++){
49     $S_number[0][$i]=1;
50 }
51 for(my $i=1;$i<=$m;$i++){
52     $H[$i][0]=(-2)*$i;
53     $S_number[$i][0]=1;
54 }
55 ##三维数组S的初始化
56 @{$S[0][0]}=();
57 for(my $i=1;$i<=$n;$i++){
58     @{$S[0][$i]}='left';
59 }
60 for (my $i=1;$i<=$m;$i++){
61     @{$S[$i][0]}='up';
62 }
63 ##计算H[j]子程序；
64 sub max{
65     my %status;     
66     my $a=$_[0];my $b=$_[1];
67     if($F[0][$a-1] eq $F[1][$b-1]){
68         $status{left_up}=$H[$a-1][$b-1]+1;
69     }
70     else{
71         $status{left_up}=$H[$a-1][$b-1]-1;
72     }
73     $status{up}=$H[$a-1][$b]-2;
74     $status{left}=$H[$a][$b-1]-2;
75     my $max=$status{left_up};
76     my @score_status;
77 ##利用高水线，取得三种状态下的最大值；
78     foreach (keys %status){
79         if($status{$_} > $max){
80             $max=$status{$_};
81         }
82     }
83     push @score_status,$max;
84     foreach (keys %status){
85         if($status{$_} == $max){
86             push @score_status,$_;
87         }
88     }
89     @score_status;
90 }
91 ##计算每个H[j]的值并将每个值对应的状态存入三维数组S中；
92 for(my $i=1;$i<=$m;$i++){
93     for(my $j=1;$j<=$n;$j++){
94         my @temp=&max($i,$j);
95         $H[$i][$j]=shift @temp;
96         my $status_number=@temp;
97         $S_number[$i][$j]=$status_number;
98         @{$S[$i][$j]}=@temp;
99     }
100 }
101 ##将二维数组复制一份copy为S_number_copy;
102 for(my $i=0;$i<=$m;$i++){
103     for(my $j=0;$j<=$n;$j++){
104         $S_number_copy[$i][$j]=$S_number[$i][$j];
105     }
106 }
107 ##输出每个H[j]的值以及对应的状态；
108 for (my $i=0;$i<=$m;$i++){
109     print "@{$H[$i]}\n";
110 }
111 for(my $i=0;$i<=$m;$i++){
112     for(my $j=0;$j<=$n;$j++){
113         print "{"."@{$S[$i][$j]}"."}";
114     }
115     print "\n";
116 }
117 ##依据三维数组S，寻找路径上的下一个位置
118 sub next{
119     my $status=$_[0];
120     my @position;
121         if ($status eq 'left'){
122         my $i=$_[1];my $j=$_[2]-1;
123         @position=($i,$j);
124     }
125         if ($status eq 'up'){
126         my $i=$_[1]-1;my $j=$_[2];
127         @position=($i,$j);
128     }
129         if ($status eq 'left_up'){
130         my $i=$_[1]-1;my $j=$_[2]-1;
131         @position=($i,$j);        
132     }
133     return @position;
134 }
135
136 ##深度优先遍历路径，每个路径存放入二维数组route中；
137 ROUTE:{
138     while ($m>0 || $n>0){
139      my $status;
140      my $status_number;
141      my @next_position;
142      $status=$S[$m][$n][0];
143      my $temp=shift @{$S[$m][$n]};
144      push @{$S[$m][$n]},$temp;
145      $S_number[$m][$n]=$S_number[$m][$n]-1;
146      $status_number=$S_number[$m][$n];
147      push @{$route[$route_number]},$status;
148      push @{$route_node[$route_number]},($m,$n);
149      push @{$route_node_number[$route_number]},$status_number;
150      @next_position=&next($status,$m,$n);
151      $m=$next_position[0];$n=$next_position[1];
152     }
153     $route_number++;
154     @{$route[$route_number-1]}=reverse @{$route[$route_number-1]};
155     @{$route_node_number[$route_number-1]}=reverse @{$route_node_number[$route_number-1]};
156     @{$route_node[$route_number-1]}=reverse @{$route_node[$route_number-1]};
157     my $count=0;
158     foreach(@{$route_node_number[$route_number-1]}){
159          $m=$route_node[$route_number-1][$count*2+1];$n=$route_node[$route_number-1][$count*2];     
160          my $number=@{$route[$route_number-1]};
161         my $last=$number-1;                                                                        
162          if ($_>0){
163         ##若回溯的节点是顶点
164             if ($count==$last){
165                 redo ROUTE;
166             }else{
167                 my @temp1=@{$route[$route_number-1]}[$count+1..$last];
168                 push @{$route[$route_number]},reverse @temp1;
169                 
170                 my @temp2=@{$route_node[$route_number-1]}[$count*2+2..$last*2+1];
171                 push @{$route_node[$route_number]},reverse @temp2;
172                 
173                 my @temp3=@{$route_node_number[$route_number-1]}[$count+1..$last];
174                 push @{$route_node_number[$route_number]},reverse @temp3;
175                 
176                 redo ROUTE;
177             }
178          }
179          ##若回溯的的节点之前已被回溯且可选择路径布置一条，重新将其可回溯状态设为原始可选择路径数目的n-1次；
180          if($_<0 && $S_number_copy[$m][$n]>1){
181                  my @temp1=@{$route[$route_number-1]}[$count+1..$last];
182                 push @{$route[$route_number]},reverse @temp1;
183                 
184                 my @temp2=@{$route_node[$route_number-1]}[$count*2+2..$last*2+1];
185                 push @{$route_node[$route_number]},reverse @temp2;
186                 
187                 my @temp3=@{$route_node_number[$route_number-1]}[$count+1..$last];
188                 push @{$route_node_number[$route_number]},reverse @temp3;
189                 
190                 $S_number[$m][$n]=$S_number_copy[$m][$n]-1;
191                 redo ROUTE;
192          }
193         $count++;
194     }
195 }
196
197 ##输出两个序列的比对
198 $route_number=@route;
199 print "The total route number is:$route_number\n";
200 for(my $i=0;$i<$route_number;$i++){
201     print "@{$route[$i]}\n";
202     foreach(@{$route[$i]}){
203         if ($_ eq 'left_up'){
204             push @query,$F[0][$t_q];
205             push @target,$F[1][$t_t];
206             $t_q++;
207             $t_t++;
208            
209     }
210         if ($_ eq 'left'){
211             push @query,"-";
212             push @target,$F[1][$t_t];
213             $t_t++;
214     }
215         if ($_ eq 'up'){
216             push @query,$F[0][$t_q];
217             push @target,"-";
218             $t_q++;
219     }
220  }
221     print "@query\n";
222     print "@target\n";
223     $t_q=0;$t_t=0;
224     @query=();@target=();
225 }
　　下面是两条测试的序列：
　　s=CTGAAAA;t=CAGAA;
　　输出结果如下：

　　但是大家千万不要用这个程序比对两条较长的序列，一是因为算法本身的空间和时间复杂度都很高，二是因为深度优先遍历的缘故，三是程序本身并没有进行优化。比对较长序列时肯定会导致电脑死机（除非是大型计算机）。如何改进这些问题有待进一步的学习.........
　　好，今天就到此为止，接下来将为大家献上更多用perl实现的生物信息学算法。