perl实现动态规划全局比对

代码

#!/usr/bin/perl
use warnings;

@data=();
while(<DATA>)
{
	chomp $_;
	@array=split("\t",$_);
	push (@data,[@array]);
}
$seq_A=$data[0][1];
$seq_B=$data[1][1];
$Match=$data[2][1];
$Mismatch=$data[3][1];
$Gap_symbol=$data[4][1];
@seq_A=split("",$seq_A);
@seq_B=split("",$seq_B);

$nrow=length($seq_B);
$ncol=length($seq_A);

@Dynamic=();
$Dynamic[0][0]=0;
foreach $item(1..$ncol){
	$Dynamic[0][$item]=$Dynamic[0][$item-1]+$Gap_symbol;
}
foreach $item(1..$nrow){
	$Dynamic[$item][0]=$Dynamic[$item-1][0]+$Gap_symbol;
}

foreach $item_B(1..$nrow){
	foreach $item_A(1..$ncol){
		$upd=$Dynamic[$item_B-1][$item_A]+$Gap_symbol;
		$leftd=$Dynamic[$item_B][$item_A-1]+$Gap_symbol;
		if($seq_B[$item_B-1] eq $seq_A[$item_A-1]){
			$upleftd=$Dynamic[$item_B-1][$item_A-1]+$Match;
		}else{
			$upleftd=$Dynamic[$item_B-1][$item_A-1]+$Mismatch;
		}
		$Dynamic[$item_B][$item_A]=&Max($upd,$leftd,$upleftd);
	}
}
print "\t0\t";
foreach (@seq_A){print "$_\t";}print "\n";
foreach $item_B(0..$nrow){
	if($item_B>0){print "$seq_B[$item_B-1]\t";}
	else{print "0\t";}
	foreach $item_A(0..$ncol){
		print "$Dynamic[$item_B][$item_A]\t";
	}
	print "\n";
} 

$m=$nrow;
$n=$ncol;
@matched_A=();
@matched_B=();
while($m>0||$n>0){
	&huisu;
}
print "@matched_A\n";
print "@matched_B\n";

sub Max{
	$m=$upd;
	if($m < $leftd){$m=$leftd;}
	if($m < $upleftd){$m=$upleftd;}
	$m;
}

sub huisu{
	$num=$Dynamic[$m][$n]-$Dynamic[$m-1][$n-1];
	if(($num==$Match && $seq_A[$n-1] eq $seq_B[$m-1])||($num==$Mismatch&&$seq_A[$n-1] ne $seq_B[$m-1])){
		unshift @matched_A,$seq_A[$n-1];
		unshift @matched_B,$seq_B[$m-1];
		$m--;
		$n--;
	}elsif(($Dynamic[$m][$n]-$Dynamic[$m][$n-1])==$Gap_symbol || $m==0){
		unshift @matched_A,$seq_A[$n-1];
		unshift @matched_B,'-';
		$n--;
	}else{
		unshift @matched_A,'-';
		unshift @matched_B,$seq_B[$m-1];
		$m--;
	}
}

__DATA__
SequenceA	AACGTACTCA
SequenceB	TCGTACTCT
Match	9
Mismatch	-6
Gap_symbol	-2

声明：写的不好，别喷