None available.
sub alignment
{ &opt_default(Gap=>-5,Match=>10,Miss=>-7,output=>"stdout",filename=>"alignment.cvs");
my @args=opt_get(@_);
my $gb=opt_as_gb(shift @args);
my $seq1=\$gb->{SEQ};
$gb=opt_as_gb(shift @args);
my $seq2=\$gb->{SEQ};
my $Gap=opt_val("Gap");
my $Match=opt_val("Match");
my $Miss=opt_val("Miss");
my $output=opt_val("output");
my $filename=opt_val("filename");
my %align;
my $hit_num;
my $miss_num;
my $gap_num1;
my $gap_num2;
my $H=0;
my $V=1;
my $D=2;
my $NDIR=3;
$$seq1=~tr/ \n\t[0-9]//d;
$$seq2=~tr/ \n\t[0-9]//d;
my $i=length($$seq1);
my $j=length($$seq2);
$align{score}=&find_max_score($i,$j,$seq1,$seq2,$Gap,$Match,$Miss);
while($i>0 | $j>0){
if($score_path_matrix[$i][$j]{direction}==$H){
$align{seq1}='-'.$align{seq1};
$align{seq2}=substr($$seq2,$j-1,1).$align{seq2};
$align{hit}='-'.$align{hit};
$j--;
}
if($score_path_matrix[$i][$j]{direction}==$V){
$align{seq1}=substr($$seq1,$i-1,1).$align{seq1};
$align{seq2}='-'.$align{seq2};
$align{hit}='-'.$align{hit};
$i--;
}
if($score_path_matrix[$i][$j]{direction}==$D){
my $tmp1=substr($$seq1,$i-1,1);
my $tmp2=substr($$seq2,$j-1,1);
$align{seq1}=$tmp1.$align{seq1};
$align{seq2}=$tmp2.$align{seq2};
if($tmp1 eq $tmp2){
$align{hit}='+'.$align{hit};
}
else{
$align{hit}='-'.$align{hit};
}
$j--;
$i--;
}
}
$hit_num=$align{hit}=~tr/+/+/;
$miss_num=$align{hit}=~tr/-/-/;
$gap_num1=$align{seq1}=~tr/-/-/;
$gap_num2=$align{seq2}=~tr/-/-/;
if($output eq "stdout"){
&msg_send("Seq1: ",length($$seq1),"base Seq2: ",length($$seq2),"base\n");
&msg_send("Score: $align{score} Match: $hit_num Miss: $miss_num\n");
&msg_send("Gap in seq1: $gap_num1 Gap in seq2: $gap_num2\n");
&msg_send("Match score: $Match Miss score: $Miss Gap Penalty: $Gap\n");
for(my $q=0;$q*60($align{seq1})||$q*60($align{seq2});$q++){
&msg_send(substr($align{seq1},$q*60,60),"\n");
&msg_send(substr($align{seq2},$q*60,60),"\n");
&msg_send(substr($align{hit},$q*60,60),"\n\n");
}
}
elsif($output =~ /f/){
open(FILE,">$filename");
print FILE $align{seq1},"\n";
print FILE $align{seq2},"\n";
print FILE $align{hit},"\n";
close(FILE);
}
return\% align;} | sub diffseq
{ no strict 'refs';
&opt_default(mode=>"align",Gap=>-5,Match=>10,Miss=>-7,output=>"stdout",filename=>"alignment.cvs");
my @args=opt_get(@_);
my $align;
my $seq1;
my $seq2;
my $mode=opt_val("mode");
if($mode eq "align"){
my $gb=opt_as_gb(shift @args);
$seq1=\$gb->{SEQ};
$gb=opt_as_gb(shift @args);
$seq2=\$gb->{SEQ};
}
else{
$align=shift @args;
}
my $Gap=opt_val("Gap");
my $Match=opt_val("Match");
my $Miss=opt_val("Miss");
my $output=opt_val("output");
my $filename=opt_val("filename");
my $i;
my @diff;
my @diff_pos;
my ($s1,$s2);
my ($tmp_gap,$s1_pos,$s2_pos);
if($mode eq "align"){
$align=alignment($seq1,$seq2,-Gap=>$Gap,-Match=>$Match,-Miss=>$Miss,-output=>"n");
}
while($i!=-1){
$i=index($$align{hit},'-',$i);
push(@diff_pos,$i) if($i!=-1);
$i++ if($i!=-1);
}
foreach(@diff_pos){
$s1=substr($$align{seq1},$_,1);
$s2=substr($$align{seq2},$_,1);
$tmp_gap=substr($$align{seq1},0,$_+1)=~tr/-/-/;
$s1_pos=$_-$tmp_gap+1;
$tmp_gap=substr($$align{seq2},0,$_+1)=~tr/-/-/;
$s2_pos=$_-$tmp_gap+1;
my $pos=$_+1;
push(@diff,"insertion,$pos,seq1,$s1_pos,seq2,$s2_pos,$s1,$s2") if($s1 eq "-");
push(@diff,"deletion,$pos,seq1,$s1_pos,seq2,$s2_pos,$s1,$s2") if($s2 eq "-");
push(@diff,"transversion,$pos,seq1,$s1_pos,seq2,$s2_pos,$s1,$s2") if($s1 ne "-" && $s2 ne "-");
}
if($output eq "stdout"){
foreach(@diff){
if(/(\w+),(\d+),seq1,(\d+),seq2,(\d+),(.),(.)/){
my $tmp=sprintf("%5d",$2);
&msg_send("$1\t$tmp: seq1\(",sprintf("%5d",$3),"\)\-\>seq2\(",sprintf("%5d",$4),"\) $5\-\>$6\n");
}
}
}
elsif($output=~/f/){
open(FILE,">$filename");
foreach(@diff){
print FILE $_,"\n";
}
close(FILE);
}
return\@ diff;} | sub find_max_score
{ my $i=shift;
my $j=shift;
my $seq1=shift;
my $seq2=shift;
my $Gap_Penalty=shift;
my $Match=shift;
my $Miss=shift;
my $p;
my $max_score;
my $max_dir;
my @score_tmp;
my @max_score;
my %tmp_hash;
my $H=0;
my $V=1;
my $D=2;
my $NDIR=3;
if($score_path_matrix[$i][$j]{score}){
return $score_path_matrix[$i][$j]{score};
}
else{
if($i==0){$max_score=0;$max_dir=$H;}
elsif($j==0){$max_score=0;$max_dir=$V;}
else{
$score_tmp[$V]=&find_max_score($i-1,$j,$seq1,$seq2,$Gap_Penalty,$Match,$Miss)+$Gap_Penalty;
$score_tmp[$H]=&find_max_score($i,$j-1,$seq1,$seq2,$Gap_Penalty,$Match,$Miss)+$Gap_Penalty;
if(substr($$seq1,$i-1,1) eq substr($$seq2,$j-1,1)){
$score_tmp[$D]=&find_max_score($i-1,$j-1,$seq1,$seq2,$Gap_Penalty,$Match,$Miss)+$Match;
}
else{
$score_tmp[$D]=&find_max_score($i-1,$j-1,$seq1,$seq2,$Gap_Penalty,$Match,$Miss)+$Miss;
}
%tmp_hash=();
foreach(@score_tmp){
$tmp_hash{$_}=$p;
$p++;
}
@max_score=sort{$b <=> $a}@score_tmp;
$max_dir=$tmp_hash{$max_score[0]};
$max_score=$score_tmp[$max_dir];
}
$score_path_matrix[$i][$j]{score}=$max_score;
$score_path_matrix[$i][$j]{direction}=$max_dir;
return $max_score;
}} |