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 Data::Dumper;
  36
  37 use vars qw ( @ISA @EXPORT );
  38
  39 @ISA = qw( Exporter );
  40
  41 use constant {
  42     QNAME => 0,
  43     FLAG  => 1,
  44     RNAME => 2,
  45     POS   => 3,
  46     MAPQ  => 4,
  47     CIGAR => 5,
  48     MRNM  => 6,
  49     MPOS  => 7,
  50     ISIZE => 8,
  51     SEQ   => 9,
  52     QUAL  => 10,
  53     TAG   => 11,
  54     VTYPE => 12,
  55     VALUE => 13,
  56 };
  57
  58
  59 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  60
  61
  62 sub get_entry
  63 {
  64     # Martin A. Hansen, September 2009.
  65
  66     # Parses a SAM entry from a given file handle.
  67
  68     my ( $fh,   # file handle
  69        ) = @_;
  70
  71     # Returns a list.
  72
  73     my ( $line, @fields );
  74
  75     while ( $line = <$fh> )
  76     {
  77         chomp $line;
  78
  79         next if substr( $line, 0, 1 ) eq '@';
  80
  81         @fields = split "\t", $line;
  82
  83         return wantarray ? @fields : \@fields;
  84     }
  85
  86     return;
  87 }
  88
  89
  90 sub sam2biopiece
  91 {
  92     # Martin A. Hansen, September 2009.
  93
  94     # Converts a SAM entry to a Biopiece record.
  95
  96     my ( $entry,   # SAM entry
  97        ) = @_;
  98
  99     # Returns a hashref.
 100
 101     my ( $record );
 102
 103     $record->{ 'REC_TYPE' } = 'SAM';
 104     $record->{ 'Q_ID' }     = $entry->[ QNAME ];
 105     $record->{ 'FLAG' }     = $entry->[ FLAG ];
 106     $record->{ 'S_ID' }     = $entry->[ RNAME ];
 107     $record->{ 'S_BEG' }    = $entry->[ POS ];
 108     $record->{ 'MAPQ' }     = $entry->[ MAPQ ];
 109     $record->{ 'CIGAR' }    = $entry->[ CIGAR ];
 110     $record->{ 'MRNM' }     = $entry->[ MRNM ];
 111     $record->{ 'S_BEG2' }   = $entry->[ MPOS ];
 112     $record->{ 'ISIZE' }    = $entry->[ ISIZE ];
 113     $record->{ 'SEQ' }      = $entry->[ SEQ ];
 114     $record->{ 'SCORES' }   = $entry->[ QUAL ];
 115
 116     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG' }  != 0;
 117     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG2' } != 0;
 118
 119     $record->{ 'S_END' } = $record->{ 'S_BEG' } + length( $record->{ 'SEQ' } ) - 1;
 120
 121     if ( $record->{ 'FLAG' } & 0x0010 ) {
 122         $record->{ 'STRAND' } = '+';
 123     } else {
 124         $record->{ 'STRAND' } = '-';
 125     }
 126
 127     return wantarray ? %{ $record } : $record;
 128 }
 129
 130
 131
 132
 133 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 134
 135 1;