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 sam_entry_get
  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 sam_entry_put
  91 {
  92
  93 }
  94
  95
  96 sub sam2biopiece
  97 {
  98     # Martin A. Hansen, September 2009.
  99
 100     # Converts a SAM entry to a Biopiece record.
 101
 102     my ( $entry,   # SAM entry
 103        ) = @_;
 104
 105     # Returns a hashref.
 106
 107     my ( $record );
 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->{ 'TAG' }      = $entry->[ TAG ];
 122     $record->{ 'VTYPE' }    = $entry->[ VTYPE ];
 123     $record->{ 'VALUE' }    = $entry->[ VALUE ];
 124
 125     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG' }  != 0;
 126     $record->{ 'S_BEG' } -= 1 if $record->{ 'S_BEG2' } != 0;
 127
 128     $record->{ 'S_END' } = $record->{ 'S_BEG' } + length( $record->{ 'SEQ' } ) - 1;
 129
 130     $record->{ 'READ_PAIRED' }    = $record->{ 'FLAG' } & 0x0001;
 131     $record->{ 'READ_MAPPED' }    = $record->{ 'FLAG' } & 0x0002;
 132     $record->{ 'Q_SEQ_UNMAPPED' } = $record->{ 'FLAG' } & 0x0004;
 133     $record->{ 'MATE_UNMAPPED' }  = $record->{ 'FLAG' } & 0x0008;
 134
 135     if ( $record->{ 'FLAG' } & 0x0010 ) {
 136         $record->{ 'STRAND' } = '+';
 137     } else {
 138         $record->{ 'STRAND' } = '-';
 139     }
 140
 141     $record->{ 'MATE_STRAND' }    = $record->{ 'FLAG' } & 0x0020;
 142     $record->{ '1ST_READ' }       = $record->{ 'FLAG' } & 0x0040;
 143     $record->{ '2ND_READ' }       = $record->{ 'FLAG' } & 0x0080;
 144     $record->{ 'ALIGN_PRIMARY' }  = $record->{ 'FLAG' } & 0x0100;
 145     $record->{ 'READ_FAIL' }      = $record->{ 'FLAG' } & 0x0200;
 146     $record->{ 'READ_BAD' }       = $record->{ 'FLAG' } & 0x0400;
 147
 148     return wantarray ? %{ $record } : $record;
 149 }
 150
 151
 152 sub biopiece2sam
 153 {
 154
 155 }
 156
 157
 158
 159 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 160
 161 1;