adding bzip2 support in ruby

[biopieces.git] / code_perl / Maasha / Align.pm

diff --git a/code_perl/Maasha/Align.pm b/code_perl/Maasha/Align.pm
index b0d270f4cb8a4c758a0020843ae01ead64cf19d7..262ee48255634929b0e7383224952b3800fc88aa 100644 (file)
--- a/code_perl/Maasha/Align.pm
+++ b/code_perl/Maasha/Align.pm
@@ -28,6 +28,7 @@ package Maasha::Align;
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
 
+use warnings;
 use strict;
 use Data::Dumper;
 use IPC::Open2;
@@ -38,8 +39,8 @@ use Maasha::Seq;
 use vars qw ( @ISA @EXPORT );
 
 use constant {
-    HEAD  => 0,
-    SEQ   => 1,
+    SEQ_NAME => 0,
+    SEQ      => 1,
 };
 
 @ISA = qw( Exporter );
@@ -68,7 +69,7 @@ sub align
     $muscle_args  = "-quiet";
     $muscle_args .= $args if $args;
 
-    @aligned_entries = &align_muscle( $entries, $muscle_args );
+    @aligned_entries = align_muscle( $entries, $muscle_args );
 
     return wantarray ? @aligned_entries : \@aligned_entries;
 }
@@ -94,11 +95,11 @@ sub align_muscle
 
     $pid = open2( $fh_out, $fh_in, $cmd );
 
-    map { &Maasha::Fasta::put_entry( $_, $fh_in ) } @{ $entries };
+    map { Maasha::Fasta::put_entry( $_, $fh_in ) } @{ $entries };
 
     close $fh_in;
 
-    while ( $entry = &Maasha::Fasta::get_entry( $fh_out ) ) {
+    while ( $entry = Maasha::Fasta::get_entry( $fh_out ) ) {
         push @aligned_entries, $entry;
     }
 
@@ -126,13 +127,13 @@ sub align_print_pairwise
 
     my ( @entries, $ruler1, $ruler2, $pins );
 
-    $ruler1 = &align_ruler( $entry1, 1 );
-    $ruler2 = &align_ruler( $entry2, 1 );
-    $pins   = &align_pins( $entry1, $entry2 );
+    $ruler1 = align_ruler( $entry1, 1 );
+    $ruler2 = align_ruler( $entry2, 1 );
+    $pins   = align_pins( $entry1, $entry2 );
 
     push @entries, $ruler1, $entry1, $pins, $entry2, $ruler2;
 
-    &align_print( \@entries, $fh, $wrap );
+    align_print( \@entries, $fh, $wrap );
 }
 
 
@@ -153,13 +154,13 @@ sub align_print_multi
 
     my ( @entries, $ruler, $consensus );
 
-    $ruler     = &align_ruler( $entries->[ 0 ] );
-    $consensus = &align_consensus( $entries ) if not $no_cons;
+    $ruler     = align_ruler( $entries->[ 0 ] );
+    $consensus = align_consensus( $entries ) if not $no_cons;
 
     unshift @{ $entries }, $ruler if not $no_ruler;
-    push    @{ $entries }, $consensus;
+    push    @{ $entries }, $consensus if not $no_cons;
 
-    &align_print( $entries, $fh, $wrap );
+    align_print( $entries, $fh, $wrap );
 }
 
 
@@ -180,20 +181,20 @@ sub align_print
 
     $max = 0;
 
-    map { $max = length $_->[ HEAD ] if length $_->[ HEAD ] > $max } @{ $entries };
+    map { $max = length $_->[ SEQ_NAME ] if length $_->[ SEQ_NAME ] > $max } @{ $entries };
 
-    $blocks = &align_wrap( $entries, $wrap );
+    $blocks = align_wrap( $entries, $wrap );
 
     foreach $block ( @{ $blocks } )
     {
         foreach $entry ( @{ $block } )
         {
-            $entry->[ HEAD ] =~ s/stats|ruler|consensus//;                                                                          
+            $entry->[ SEQ_NAME ] =~ s/stats|ruler|consensus//;                                                                          
 
             if ( $fh ) {
-                print $fh $entry->[ HEAD ], " " x ( $max + 3 - length $entry->[ HEAD ] ), $entry->[ SEQ ], "\n";                                 
+                print $fh $entry->[ SEQ_NAME ], " " x ( $max + 3 - length $entry->[ SEQ_NAME ] ), $entry->[ SEQ ], "\n";                                 
             } else {
-                print $entry->[ HEAD ], " " x ( $max + 3 - length $entry->[ HEAD ] ), $entry->[ SEQ ], "\n";                                 
+                print $entry->[ SEQ_NAME ], " " x ( $max + 3 - length $entry->[ SEQ_NAME ] ), $entry->[ SEQ ], "\n";                                 
             }
         }
     }
@@ -225,7 +226,7 @@ sub align_wrap
 
         for ( $c = 0; $c < @{ $entries }; $c++ )
         {
-            if ( $entries->[ $c ]->[ HEAD ] eq "ruler" )
+            if ( $entries->[ $c ]->[ SEQ_NAME ] eq "ruler" )
             {
                 $ruler = substr $entries->[ $c ]->[ SEQ ], $i, $wrap;
 
@@ -241,7 +242,7 @@ sub align_wrap
             }
             else
             {
-                push @lines, [ $entries->[ $c ]->[ HEAD ], substr $entries->[ $c ]->[ SEQ ], $i, $wrap ];
+                push @lines, [ $entries->[ $c ]->[ SEQ_NAME ], substr $entries->[ $c ]->[ SEQ ], $i, $wrap ];
             }
         }
 
@@ -269,20 +270,20 @@ sub align_pins
 
     my ( $blosum, $i, $char1, $char2, $pins );
 
-    $type ||= &Maasha::Seq::seq_guess_type( $entry1->[ SEQ ] );
+    $type ||= Maasha::Seq::seq_guess_type( $entry1->[ SEQ ] );
 
-    $blosum = &blosum_read() if $type =~ /protein/;
+    $blosum = blosum_read() if $type =~ "PROTEIN";
 
     for ( $i = 0; $i < length $entry1->[ SEQ ]; $i++ )
     {
-        $char1 = substr $entry1->[ SEQ ], $i, 1;
-        $char2 = substr $entry2->[ SEQ ], $i, 1;
+        $char1 = uc substr $entry1->[ SEQ ], $i, 1;
+        $char2 = uc substr $entry2->[ SEQ ], $i, 1;
 
         if ( $blosum and $char1 eq $char2 ) {
             $pins .= $char1;
         } elsif ( $char1 eq $char2 ) {
             $pins .= "|";
-        } elsif ( $blosum and $blosum->{ $char1 }->{ $char2 } > 0 ) {
+        } elsif ( $blosum and defined $blosum->{ $char1 }->{ $char2 } and $blosum->{ $char1 }->{ $char2 } > 0 ) {
             $pins .= "+";
         } else {
             $pins .= " ";
@@ -373,7 +374,7 @@ sub align_consensus
 
     my ( $bit_max, $data, $pos, $char, $score, $entry );
 
-    $type    ||= &Maasha::Seq::seq_guess_type( $entries->[ 0 ]->[ SEQ ] );
+    $type    ||= Maasha::Seq::seq_guess_type( $entries->[ 0 ]->[ SEQ ] );
     $min_sim ||= 50;
 
     if ( $type =~ /protein/ ) {
@@ -382,7 +383,7 @@ sub align_consensus
         $bit_max = 2;
     }
 
-    $data = &Maasha::Seq::seqlogo_calc( $bit_max, $entries );
+    $data = Maasha::Seq::seqlogo_calc( $bit_max, $entries );
 
     foreach $pos ( @{ $data } )
     {
@@ -402,7 +403,7 @@ sub align_consensus
         }
     }
 
-    $entry->[ HEAD ] = "Consensus: $min_sim%";
+    $entry->[ SEQ_NAME ] = "Consensus: $min_sim%";
 
     return wantarray ? @{ $entry } : $entry;
 }
@@ -446,7 +447,7 @@ sub align_sim_global
         $match_tot++ if substr( $entry1->[ SEQ ], $i, 1 ) eq substr( $entry2->[ SEQ ], $i, 1 );
     }
 
-    $min = &Maasha::Calc::min( $len1, $len2 );
+    $min = Maasha::Calc::min( $len1, $len2 );
 
     $sim = sprintf( "%.2f", ( $match_tot / $min ) * 100 );
 
@@ -462,8 +463,8 @@ sub align_tile
     # using pairwise alignments. The result is returned as a list of
     # aligned FASTA entries.
     
-    my ( $ref_entry,   # reference entry as [ HEAD, SEQ ] tuple
-         $q_entries,   # list of [ HEAD, SEQ ] tuples
+    my ( $ref_entry,   # reference entry as [ SEQ_NAME, SEQ ] tuple
+         $q_entries,   # list of [ SEQ_NAME, SEQ ] tuples
          $args,        # argument hash
        ) = @_;
 
@@ -477,27 +478,27 @@ sub align_tile
     {
         $seq1 = $entry->[ SEQ ];
 
-        $type = &Maasha::Seq::seq_guess_type( $seq1 );
+        $type = Maasha::Seq::seq_guess_type( $seq1 );
 
-        if ( $type eq "rna" ) {
-            $seq2 = &Maasha::Seq::rna_revcomp( $seq1 );
-        } elsif ( $type eq "dna" ) {
-            $seq2 = &Maasha::Seq::dna_revcomp( $seq1 );
+        if ( $type eq "RNA" ) {
+            $seq2 = Maasha::Seq::rna_revcomp( $seq1 );
+        } elsif ( $type eq "DNA" ) {
+            $seq2 = Maasha::Seq::dna_revcomp( $seq1 );
         } else {
-            &Maasha::Common::error( qq(Bad sequence type->$type) );
+            Maasha::Common::error( qq(Bad sequence type->$type) );
         }
 
-        $align1 = &Maasha::Align::align_muscle( [ $ref_entry, [ $entry->[ HEAD ] . "_+", $seq1 ] ], "-quiet -maxiters 1" );
-        $align2 = &Maasha::Align::align_muscle( [ $ref_entry, [ $entry->[ HEAD ] . "_-", $seq2 ] ], "-quiet -maxiters 1" );
+        $align1 = Maasha::Align::align_muscle( [ $ref_entry, [ $entry->[ SEQ_NAME ] . "_+", $seq1 ] ], "-quiet -maxiters 1" );
+        $align2 = Maasha::Align::align_muscle( [ $ref_entry, [ $entry->[ SEQ_NAME ] . "_-", $seq2 ] ], "-quiet -maxiters 1" );
 
         if ( $args->{ "supress_indels" } )
         {
-            &align_supress_indels( $align1 );
-            &align_supress_indels( $align2 );
+            align_supress_indels( $align1 );
+            align_supress_indels( $align2 );
         }
 
-        $sim1 = &Maasha::Align::align_sim_global( $align1->[ 0 ], $align1->[ 1 ] );
-        $sim2 = &Maasha::Align::align_sim_global( $align2->[ 0 ], $align2->[ 1 ] );
+        $sim1 = Maasha::Align::align_sim_global( $align1->[ 0 ], $align1->[ 1 ] );
+        $sim2 = Maasha::Align::align_sim_global( $align2->[ 0 ], $align2->[ 1 ] );
 
         if ( $sim1 < $args->{ "identity" } and $sim2 < $args->{ "identity" } )
         {
@@ -509,7 +510,7 @@ sub align_tile
 
             $align1->[ 1 ]->[ SEQ ] =~ s/-{$gaps}$// if $gaps;
             
-            $entry->[ HEAD ] = "$align1->[ 1 ]->[ HEAD ]_$sim1";
+            $entry->[ SEQ_NAME ] = "$align1->[ 1 ]->[ SEQ_NAME ]_$sim1";
             $entry->[ SEQ ]  = $align1->[ 1 ]->[ SEQ ];
 
             push @entries, $entry;
@@ -520,7 +521,7 @@ sub align_tile
 
             $align2->[ 1 ]->[ SEQ ] =~ s/-{$gaps}$// if $gaps;
 
-            $entry->[ HEAD ] = "$align2->[ 1 ]->[ HEAD ]_$sim2";
+            $entry->[ SEQ_NAME ] = "$align2->[ 1 ]->[ SEQ_NAME ]_$sim2";
             $entry->[ SEQ ]  = $align2->[ 1 ]->[ SEQ ];
 
             push @entries, $entry;