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，其中每个值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所有可选状态的数目，与三维数组S相对应；
8 my @S_number_copy; ##存放S_number的一个copy,用于判断回溯节点的原始可选择路径数目；
9 my @route; ##存放从H到H的路径的边；
10 my @route_node; ##存放从H到H的路径的节点；
11 my @route_node_number; ##存放从H到H的路径的节点的剩余状态数；
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},@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},@second;
39 }
40 close TARGET;
41
42 ##二维数组H和S_number的初始化；
43 $m=@{$F};
44 $n=@{$F};
45 $H=0;
46 $S_number=0;
47 @{$H}=map {$_*(-2)} my @n=(0..$n);
48 for(my $i=1;$i<=$n;$i++){
49 $S_number[$i]=1;
50 }
51 for(my $i=1;$i<=$m;$i++){
52 $H[$i]=(-2)*$i;
53 $S_number[$i]=1;
54 }
55 ##三维数组S的初始化
56 @{$S}=();
57 for(my $i=1;$i<=$n;$i++){
58 @{$S[$i]}='left';
59 }
60 for (my $i=1;$i<=$m;$i++){
61 @{$S[$i]}='up';
62 }
63 ##计算H子程序；
64 sub max{
65 my %status;
66 my $a=$_;my $b=$_;
67 if($F[$a-1] eq $F[$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的值并将每个值对应的状态存入三维数组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的值以及对应的状态；
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=$_;
120 my @position;
121 if ($status eq 'left'){
122 my $i=$_;my $j=$_-1;
123 @position=($i,$j);
124 }
125 if ($status eq 'up'){
126 my $i=$_-1;my $j=$_;
127 @position=($i,$j);
128 }
129 if ($status eq 'left_up'){
130 my $i=$_-1;my $j=$_-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];
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;$n=$next_position;
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[$t_q];
205 push @target,$F[$t_t];
206 $t_q++;
207 $t_t++;
208
209 }
210 if ($_ eq 'left'){
211 push @query,"-";
212 push @target,$F[$t_t];
213 $t_t++;
214 }
215 if ($_ eq 'up'){
216 push @query,$F[$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实现的生物信息学算法。
页:
[1]