code_perl/Maasha/SAM.pm

   1 package Maasha::SAM;
   2
   3 # Copyright (C) 2009 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines to handle SAM and BAM format.
  26 # http://samtools.sourceforge.net/
  27
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31
  32 use warnings;
  33 use strict;
  34 use Maasha::Common;
  35 use Maasha::Fastq;
  36 use Data::Dumper;
  37
  38 use vars qw( @ISA @EXPORT );
  39
  40 @ISA = qw( Exporter );
  41
  42 use constant {
  43     QNAME => 0,
  44     FLAG  => 1,
  45     RNAME => 2,
  46     POS   => 3,
  47     MAPQ  => 4,
  48     CIGAR => 5,
  49     MRNM  => 6,
  50     MPOS  => 7,
  51     ISIZE => 8,
  52     SEQ   => 9,
  53     QUAL  => 10,
  54 };
  55
  56
  57 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  58
  59
  60 sub sam_entry_get
  61 {
  62     # Martin A. Hansen, September 2009.
  63
  64     # Parses a SAM entry from a given file handle.
  65
  66     my ( $fh,   # file handle
  67        ) = @_;
  68
  69     # Returns a list.
  70
  71     my ( $line, @fields );
  72
  73     while ( $line = <$fh> )
  74     {
  75         chomp $line;
  76
  77         next if substr( $line, 0, 1 ) eq '@';
  78
  79         @fields = split "\t", $line;
  80
  81         return wantarray ? @fields : \@fields;
  82     }
  83
  84     return;
  85 }
  86
  87
  88 sub sam_entry_put
  89 {
  90
  91 }
  92
  93
  94 sub sam2biopiece
  95 {
  96     # Martin A. Hansen, September 2009.
  97
  98     # Converts a SAM entry to a Biopiece record.
  99
 100     my ( $entry,   # SAM entry
 101        ) = @_;
 102
 103     # Returns a hashref.
 104
 105     my ( $record, $i, $tag, $vtype, $value );
 106
 107     return if $entry->[ RNAME ] eq '*';
 108
 109     $record->{ 'REC_TYPE' } = 'SAM';
 110     $record->{ 'Q_ID' }     = $entry->[ QNAME ];
 111     $record->{ 'FLAG' }     = $entry->[ FLAG ];
 112     $record->{ 'S_ID' }     = $entry->[ RNAME ];
 113     $record->{ 'S_BEG' }    = $entry->[ POS ];
 114     $record->{ 'MAPQ' }     = $entry->[ MAPQ ];
 115     $record->{ 'CIGAR' }    = $entry->[ CIGAR ];
 116     $record->{ 'MRNM' }     = $entry->[ MRNM ];
 117     $record->{ 'S_BEG2' }   = $entry->[ MPOS ];
 118     $record->{ 'ISIZE' }    = $entry->[ ISIZE ];
 119     $record->{ 'SEQ' }      = $entry->[ SEQ ];
 120     $record->{ 'SCORES' }   = $entry->[ QUAL ];
 121     $record->{ 'SCORES' }   =~ s/(.)/chr( ( ord( $1 ) - 33 ) + 64 )/ge; # convert phred-33 to phred-64 scores
 122     $record->{ 'SCORE' }    = sprintf( "%.2f", Maasha::Fastq::phred_str_mean( $entry->[ QUAL ] ) );
 123
 124     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG' }  != 0;
 125     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG2' } != 0;
 126
 127     $record->{ 'S_END' } = $record->{ 'S_BEG' } + length( $record->{ 'SEQ' } ) - 1;
 128
 129     $record->{ 'READ_PAIRED' }    = $record->{ 'FLAG' } & 0x0001;
 130     $record->{ 'READ_MAPPED' }    = $record->{ 'FLAG' } & 0x0002;
 131     $record->{ 'Q_SEQ_UNMAPPED' } = $record->{ 'FLAG' } & 0x0004;
 132     $record->{ 'MATE_UNMAPPED' }  = $record->{ 'FLAG' } & 0x0008;
 133
 134     if ( $record->{ 'FLAG' } & 0x0010 ) {   # 0 for forward, 1 for reverse
 135         $record->{ 'STRAND' } = '-';
 136     } else {
 137         $record->{ 'STRAND' } = '+';
 138     }
 139
 140     $record->{ 'MATE_STRAND' }    = $record->{ 'FLAG' } & 0x0020;
 141     $record->{ '1ST_READ' }       = $record->{ 'FLAG' } & 0x0040;
 142     $record->{ '2ND_READ' }       = $record->{ 'FLAG' } & 0x0080;
 143     $record->{ 'ALIGN_PRIMARY' }  = $record->{ 'FLAG' } & 0x0100;
 144     $record->{ 'READ_FAIL' }      = $record->{ 'FLAG' } & 0x0200;
 145     $record->{ 'READ_BAD' }       = $record->{ 'FLAG' } & 0x0400;
 146
 147     if ( $record->{ 'READ_PAIRED' } ) {
 148         $record->{ 'BLOCK_COUNT' } = 2;
 149     } else {
 150         $record->{ 'BLOCK_COUNT' } = 1;
 151     }
 152
 153     # extra fields - (hate hate hate SAM!!!)
 154
 155     for ( $i = QUAL + 1; $i < @{ $entry }; $i++ )
 156     {
 157         ( $tag, $vtype, $value ) = split /:/, $entry->[ $i ];
 158
 159         if ( $tag eq "MD" or $tag eq "NM" ) {
 160             $record->{ $tag } = $value;
 161         }
 162     }
 163
 164     if ( defined $record->{ "NM" } and $record->{ "NM" } > 0 ) {
 165         $record->{ "ALIGN" } = sam2align( $record->{ 'SEQ' }, $record->{ 'CIGAR' }, $record->{ 'NM' }, $record->{ 'MD' } );
 166     } else {
 167         $record->{ "ALIGN" } = ".";
 168     }
 169
 170     return wantarray ? %{ $record } : $record;
 171 }
 172
 173
 174 sub sam2align
 175 {
 176     my ( $seq,
 177          $cigar,
 178          $nm,
 179          $md,
 180        ) = @_;
 181
 182     # Returns a string.
 183
 184     my ( $offset, $len, $op, $i, $nt, $nt2, @nts, $insertions, @align, $aln_str );
 185
 186     $offset     = 0;
 187     $insertions = 0;
 188
 189     while ( length( $cigar ) > 0 )
 190     {
 191         if ( $cigar =~ s/^(\d+)(\w)// )
 192         {
 193             $len = $1;
 194             $op  = $2;
 195
 196             if ( $op eq 'I' )
 197             {
 198                 for ( $i = 0; $i < $len; $i++ )
 199                 {
 200                     $nt = substr $seq, $offset + $i, 1;
 201
 202                     push @align, [ $offset + $i, "->$nt" ];
 203
 204                     $nm--;
 205                 }
 206             }
 207
 208             $offset += $len;
 209         }
 210     }
 211
 212     $offset = 0;
 213
 214     while ( $nm > 0 )
 215     {
 216         if ( $md =~ s/^(\d+)// )   # match
 217         {
 218             $offset += $1;
 219         }
 220         elsif( $md =~ s/^\^([ATCGN]+)// )  # insertions
 221         {
 222             @nts = split //, $1;
 223
 224             foreach $nt ( @nts )
 225             {
 226                 $nt2 = substr $seq, $offset, 1;
 227
 228                 push @align, [ $offset, "$nt2>-" ];
 229
 230                 $offset++;
 231
 232                 $insertions++;
 233
 234                 $nm--;
 235             }
 236         }
 237         elsif ( $md =~ s/^([ATCGN]+)// ) # mismatch
 238         {
 239             @nts = split //, $1;
 240
 241             foreach $nt ( @nts )
 242             {
 243                 $nt2 = substr $seq, $offset - $insertions, 1;
 244
 245                 push @align, [ $offset, "$nt>$nt2" ];
 246
 247                 $offset++;
 248
 249                 $nm--;
 250             }
 251         }
 252         else
 253         {
 254             Maasha::Common::error( qq(Bad SAM field -> MD: "$md") );
 255         }
 256     }
 257
 258     $aln_str = "";
 259
 260     @align = sort { $a->[ 0 ] <=> $b->[ 0 ] } @align;
 261
 262     map { $aln_str .= "$_->[ 0 ]:$_->[ 1 ]," } @align;
 263
 264     $aln_str =~ s/,$//;
 265
 266     return $aln_str;
 267 }
 268
 269
 270 sub biopiece2sam
 271 {
 272
 273 }
 274
 275
 276
 277 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 278
 279 1;