adding bzip2 support in ruby

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

diff --git a/code_perl/Maasha/SAM.pm b/code_perl/Maasha/SAM.pm
index e794d75f69b066b8c94ca4419505166efb5315f2..ff8cb66e61dc8682db2aea75bec1b84fe2bbbdf7 100644 (file)
--- a/code_perl/Maasha/SAM.pm
+++ b/code_perl/Maasha/SAM.pm
@@ -118,6 +118,7 @@ sub sam2biopiece
     $record->{ 'ISIZE' }    = $entry->[ ISIZE ];
     $record->{ 'SEQ' }      = $entry->[ SEQ ];
     $record->{ 'SCORES' }   = $entry->[ QUAL ];
+    $record->{ 'SCORES' }   =~ s/(.)/chr( ( ord( $1 ) - 33 ) + 64 )/ge; # convert phred-33 to phred-64 scores
     $record->{ 'SCORE' }    = sprintf( "%.2f", Maasha::Fastq::phred_str_mean( $entry->[ QUAL ] ) );
 
     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG' }  != 0;
@@ -130,10 +131,10 @@ sub sam2biopiece
     $record->{ 'Q_SEQ_UNMAPPED' } = $record->{ 'FLAG' } & 0x0004;
     $record->{ 'MATE_UNMAPPED' }  = $record->{ 'FLAG' } & 0x0008;
 
-    if ( $record->{ 'FLAG' } & 0x0010 ) {
-        $record->{ 'STRAND' } = '+';
-    } else {
+    if ( $record->{ 'FLAG' } & 0x0010 ) {   # 0 for forward, 1 for reverse
         $record->{ 'STRAND' } = '-';
+    } else {
+        $record->{ 'STRAND' } = '+';
     }
 
     $record->{ 'MATE_STRAND' }    = $record->{ 'FLAG' } & 0x0020;
@@ -143,6 +144,11 @@ sub sam2biopiece
     $record->{ 'READ_FAIL' }      = $record->{ 'FLAG' } & 0x0200;
     $record->{ 'READ_BAD' }       = $record->{ 'FLAG' } & 0x0400;
 
+    if ( $record->{ 'READ_PAIRED' } ) {
+        $record->{ 'BLOCK_COUNT' } = 2;
+    } else {
+        $record->{ 'BLOCK_COUNT' } = 1;
+    }
 
     # extra fields - (hate hate hate SAM!!!)
     
@@ -175,9 +181,10 @@ sub sam2align
 
     # Returns a string.
 
-    my ( $offset, $len, $op, $i, $nt, $nt2, @align, $aln_str );
+    my ( $offset, $len, $op, $i, $nt, $nt2, @nts, $insertions, @align, $aln_str );
 
-    $offset = 0;
+    $offset     = 0;
+    $insertions = 0;
 
     while ( length( $cigar ) > 0 )
     {
@@ -206,31 +213,41 @@ sub sam2align
 
     while ( $nm > 0 )
     {
-        if ( $md =~ s/^(\d+)// )
+        if ( $md =~ s/^(\d+)// )   # match
         {
             $offset += $1;
         }
-        elsif ( $md =~ s/^(\w)// )
+        elsif( $md =~ s/^\^([ATCGN]+)// )  # insertions
         {
-            $nt = $1;
+            @nts = split //, $1;
 
-            $nt2 = substr $seq, $offset, 1;
+            foreach $nt ( @nts )
+            {
+                $nt2 = substr $seq, $offset, 1;
 
-            push @align, [ $offset, "$nt>$nt2" ];
+                push @align, [ $offset, "$nt2>-" ];
+        
+                $offset++;
 
-            $offset++;
+                $insertions++;
 
-            $nm--;
+                $nm--;
+            }
         }
-        elsif( $md =~ s/^\^// )
+        elsif ( $md =~ s/^([ATCGN]+)// ) # mismatch
         {
-            $nt = substr $seq, $offset, 1;
+            @nts = split //, $1;
 
-            push @align, [ $offset, "$nt>-" ];
-        
-            $offset++;
+            foreach $nt ( @nts )
+            {
+                $nt2 = substr $seq, $offset - $insertions, 1;
 
-            $nm--;
+                push @align, [ $offset, "$nt>$nt2" ];
+
+                $offset++;
+
+                $nm--;
+            }
         }
         else
         {
@@ -238,7 +255,6 @@ sub sam2align
         }
     }
 
-
     $aln_str = "";
 
     @align = sort { $a->[ 0 ] <=> $b->[ 0 ] } @align;