G::Seq

Align

G::Seq::Align - Perl extension for blah blah blah

  blah blah blah

Stub documentation for G::Seq::Align was created by h2xs. It looks like the
author of the extension was negligent enough to leave the stub
unedited.

Blah blah blah.

Methods description

Methods code

    my $self = shift;

    &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;

    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;

    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;
    }

    my $pkg = shift;
    my $filename = shift;
    my $option = shift;
    my $this;

    return $this;

General documentation

A. U. Thor, a.u.thor@a.galaxy.far.far.away

perl(1).