globalblast2.pl

实现一个生物信息学上面的blast的比对,有很好的借鉴作用

#!/usr/bin/perl
use strict;
#use FileHandle;
#my $fh2 = new FileHandle(">f:\\data1.txt");
open  IN_FILE,  "c:\\perl\\cmyc.fasta" or die("culd not open the file");
open  OUT_FILE, ">c:\\perl\\123.txt" or die("culd not open the file");
my @line =<IN_FILE>;
my $string=join("",@line);
my @result = split  />/,$string;
my @seqname;
my @seq;
my $i;
my $j;
my @len;
my @seqall;
my $long;
for ($j=0; $j<=$#result ;$j++)
{
@seqall=split /\n/, $result[$j];
$seqname[$j]=$seqall[0];
 for($i=1;$i<=$#seqall+1;$i++)  
{ 
     $seq[$j].= $seqall[$i];
}
$len[$j]=length($seq[$j]);
}
my $m;
my $n;
#取两个不同序列比对
for($m=1;$m<=$#result;$m++)
{for($n=$m+1;$n<=$#result;$n++)

{ 
my @matrix;
my @score;
my @array;
$matrix[0][0]=0;
#初始化数组
#if($len[$m]>=$len[$n])
#{$long=$len[$m];}
#else {$long=$len[$n];}
#for($i=0;$i<$long;$i++)
#{ $array[$i]=' ';
#}
#初始化矩阵
for ($i=1;$i<$len[$m]+1 ;$i++) 
{
	$matrix[$i][0]=$matrix[$i-1][0]-1;
	$score[$i][0]=-1;
}

for ($i=1;$i<$len[$n]+1 ;$i++) 
{
	$matrix[0][$i]=$matrix[0][$i-1]-1;
	$score[0][$i]=-1;
}
 for ($i=1;$i<=$len[$m] ;$i++)
{
 for ($j=1;$j<=$len[$n] ;$j++)
 { 
	 my $s1=substr($seq[$m],$i-1,1);
	 my $s2=substr($seq[$n],$j-1,1);
	 if ($s1 eq $ s2) 
	 {
		$score[$i][$j]=2;
	 }
		else  
	 {
		$score[$i][$j]=-1;
	 }
			
	 if ($matrix[$i-1][$j-1]+$score[$i][$j]>$matrix[$i-1][$j]-1) 
	 {
			if ($matrix[$i-1][$j-1]+$score[$i][$j]>$matrix[$i][$j-1]-1) 
			{
				$matrix[$i][$j]=$matrix[$i-1][$j-1]+$score[$i][$j];
			}
			else
		    {
				$matrix[$i][$j]=$matrix[$i][$j-1]-1;
			}
	 }
	 else
	{
		    if ($matrix[$i-1][$j]-1>$matrix[$i][$j-1]-1) 
			{
				$matrix[$i][$j]=$matrix[$i-1][$j]-1;
		    }
			else
		    {
				$matrix[$i][$j]=$matrix[$i][$j-1]-1;
			}
    }
 }
}
my $prestr;
my $bacstr;
my $gap="-";
#my @array;
my $seqn='';
$i=$len[$m];
$j=$len[$n];
my $x;
my $score=$matrix[$i][$j];
my $identity;
my $max;
for ($i=$len[$m],$j=$len[$n];$i>0 and $j>0 ;)
{
	
	
	if ($matrix[$i][$j] eq $matrix[$i-1][$j-1]+$score[$i][$j] ) 
	{     
                $seqn='|'.$seqn   ;   
                $i--;
		            $j--;
		            $score+=$matrix[$i-1][$j-1];
		            $x++;
		
	}
	if ($matrix[$i][$j] eq $matrix[$i][$j-1]-1 ) 
	{
		$j--;
		$prestr=substr($seq[$m],0,$i);
		$bacstr=substr($seq[$m],$i,$len[1]-$i);
		$len[$m]++;
		$seq[$m]=$prestr.$gap;
		$seq[$m].=$bacstr;
		$seqn=' '.$seqn;
	  $score+=$matrix[$i][$j-1];
	}
		
	if ($matrix[$i][$j] eq $matrix[$i-1][$j]-1) 
	{
			 $i--;
			 $prestr=substr($seq[$n],0,$j);
			 $bacstr=substr($seq[$n],$j,$len[2]-$j);
			 $len[$n]++;
			 $seq[$n]=$prestr.$gap;
			 $seq[$n].=$bacstr;
			 $seqn=' '.$seqn;
			 $score+=$matrix[$i-1][$j];
	}
	if($len[$m]>$len[$n])
	{$max=$len[$m];
		$identity=$x/$max;}
	else {
		$max=$len[$n];
		$identity=$x/$max;}
		
}
if ($i eq 0 and $j ne 0)
{
	$seq[$m]=$gap.$seq[$m];
}

if ($j eq 0 and $i ne 0)
{
	$seq[$n]=$gap.$seq[$n];
}
my $i;
my $lefta= 0;
my $leftb= 0;
print OUT_FILE "score=",$score ;
print OUT_FILE "  identity=",$identity ;
print OUT_FILE "\n" ;
for($i=0;$i<=$len[$m];$i+=60)
	  {     
                my $j;
                my $str;		
                #my $strnogap;	
		
		my $righta;
		my $rightb;
		$str= substr($seq[$m],$i,60);
		#$strnogap= $str;
		#$strnogap=~ s/\-+//g;
		$righta= $lefta+length($str)-1;
		printf OUT_FILE "Query: %-5d%s %-5d\n",$lefta,$str,$righta;
		#printf OUT_FILE "            ";
		#for($j=$i;$j<=$i+59;$j++)
                #{print OUT_FILE $array[$j] ;
        #}         
                $str= substr($seqn,$i,60);
                printf OUT_FILE "            %s\n",$str;
                #print OUT_FILE "\n";
		$str= substr($seq[$n],$i,60);
		#$strnogap= $str;
		#$strnogap=~ s/\-+//g;
		$rightb= $leftb+length($str)-1;
		printf OUT_FILE "Sbjct: %-5d%s %-5d\n",$leftb,$str,$rightb;
		print OUT_FILE "\n";
		$lefta= $righta+1;
		$leftb= $rightb+1;
	}
	}
close(IN_FILE);
close(OUT_FILE);
}