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->{ 'SCORE' }    = sprintf( "%.2f", Maasha::Fastq::phred_str_mean( $entry->[ QUAL ] ) );
 122
 123     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG' }  != 0;
 124     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG2' } != 0;
 125
 126     $record->{ 'S_END' } = $record->{ 'S_BEG' } + length( $record->{ 'SEQ' } ) - 1;
 127
 128     $record->{ 'READ_PAIRED' }    = $record->{ 'FLAG' } & 0x0001;
 129     $record->{ 'READ_MAPPED' }    = $record->{ 'FLAG' } & 0x0002;
 130     $record->{ 'Q_SEQ_UNMAPPED' } = $record->{ 'FLAG' } & 0x0004;
 131     $record->{ 'MATE_UNMAPPED' }  = $record->{ 'FLAG' } & 0x0008;
 132
 133     if ( $record->{ 'FLAG' } & 0x0010 ) {
 134         $record->{ 'STRAND' } = '+';
 135     } else {
 136         $record->{ 'STRAND' } = '-';
 137     }
 138
 139     $record->{ 'MATE_STRAND' }    = $record->{ 'FLAG' } & 0x0020;
 140     $record->{ '1ST_READ' }       = $record->{ 'FLAG' } & 0x0040;
 141     $record->{ '2ND_READ' }       = $record->{ 'FLAG' } & 0x0080;
 142     $record->{ 'ALIGN_PRIMARY' }  = $record->{ 'FLAG' } & 0x0100;
 143     $record->{ 'READ_FAIL' }      = $record->{ 'FLAG' } & 0x0200;
 144     $record->{ 'READ_BAD' }       = $record->{ 'FLAG' } & 0x0400;
 145
 146
 147     # extra fields - (hate hate hate SAM!!!)
 148
 149     for ( $i = QUAL + 1; $i < @{ $entry }; $i++ )
 150     {
 151         ( $tag, $vtype, $value ) = split /:/, $entry->[ $i ];
 152
 153         if ( $tag eq "MD" or $tag eq "NM" ) {
 154             $record->{ $tag } = $value;
 155         }
 156     }
 157
 158     if ( defined $record->{ "NM" } and $record->{ "NM" } > 0 ) {
 159         $record->{ "ALIGN" } = sam2align( $record->{ 'SEQ' }, $record->{ 'CIGAR' }, $record->{ 'NM' }, $record->{ 'MD' } );
 160     } else {
 161         $record->{ "ALIGN" } = ".";
 162     }
 163
 164     return wantarray ? %{ $record } : $record;
 165 }
 166
 167
 168 sub sam2align
 169 {
 170     my ( $seq,
 171          $cigar,
 172          $nm,
 173          $md,
 174        ) = @_;
 175
 176     # Returns a string.
 177
 178     my ( $offset, $len, $op, $i, $nt, $nt2, @align, $aln_str );
 179
 180     $offset = 0;
 181
 182     while ( length( $cigar ) > 0 )
 183     {
 184         if ( $cigar =~ s/^(\d+)(\w)// )
 185         {
 186             $len = $1;
 187             $op  = $2;
 188
 189             if ( $op eq 'I' )
 190             {
 191                 for ( $i = 0; $i < $len; $i++ )
 192                 {
 193                     $nt = substr $seq, $offset + $i, 1;
 194
 195                     push @align, [ $offset + $i, "->$nt" ];
 196
 197                     $nm--;
 198                 }
 199             }
 200
 201             $offset += $len;
 202         }
 203     }
 204
 205     $offset = 0;
 206
 207     while ( $nm > 0 )
 208     {
 209         if ( $md =~ s/^(\d+)// )
 210         {
 211             $offset += $1;
 212         }
 213         elsif ( $md =~ s/^(\w)// )
 214         {
 215             $nt = $1;
 216
 217             $nt2 = substr $seq, $offset, 1;
 218
 219             push @align, [ $offset, "$nt>$nt2" ];
 220
 221             $offset++;
 222
 223             $nm--;
 224         }
 225         elsif( $md =~ s/^\^// )
 226         {
 227             $nt = substr $seq, $offset, 1;
 228
 229             push @align, [ $offset, "$nt>-" ];
 230
 231             $offset++;
 232
 233             $nm--;
 234         }
 235         else
 236         {
 237             Maasha::Common::error( qq(Bad SAM field -> MD: "$md") );
 238         }
 239     }
 240
 241
 242     $aln_str = "";
 243
 244     @align = sort { $a->[ 0 ] <=> $b->[ 0 ] } @align;
 245
 246     map { $aln_str .= "$_->[ 0 ]:$_->[ 1 ]," } @align;
 247
 248     $aln_str =~ s/,$//;
 249
 250     return $aln_str;
 251 }
 252
 253
 254 sub biopiece2sam
 255 {
 256
 257 }
 258
 259
 260
 261 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 262
 263 1;