Blast 比对原理

my $pattern='ADFG';
my $text='FADFGH';

my $TLEN=length $text;
my $PLEN=length $pattern;

my $D=[];
for (my $t=0;$t<=$TLEN;++$t){
$D->[$t][0]=0;
}

for (my $p=0;$p<=$PLEN;++$p){
$D->[0][$p]=$p;

}

for (my $t=1;$t<=$TLEN;++$t){
for (my $p=1;$p<=$PLEN;++$p){

$D->[$t][$p]=min3(
substr($text,$t-1,1)eq substr($pattern,$p-1,1)?$D->[$t-1][$p-1]:$D->[$t-1][$p-1]+1,
$D->[$t-1][$p]+1,
$D->[$t][$p-1]+1

)

}

}

for (my $p=0;$p<=$PLEN;++$p){
for (my $t=0;$t<=$TLEN;++$t){
print $D->[$t][$p]," ";
}
print "\n";
}

打印结果

image.png

my @matches=();
my $bestscore=1000000000000;
for(my $t=1;$t<=$PLEN;++$t){

if($D->[$t][$PLEN]<$bestscore){
$bestscore=$D->[$t][$PLEN];
@matches=($t);

}
elsif($D->[$t][$PLEN]==$bestscore){
push @matches,$t;

}
}
print "\nThe best match for the pattern  $pattern\n";

print "has an edit distance of $bestscore\n";

print "and appears in the text editing location";

print "s" if(@matches >1);

print " @matches\n";

子函数

sub min3{

my($i,$j,$k)=@_;
my($temp);

$temp=($i<$j ? $i:$j);
$temp<$k ? $temp:$k;

}

打印结果

image.png