]> git.donarmstrong.com Git - biopieces.git/blob - code_perl/Maasha/SAM.pm
fixed SAM bug
[biopieces.git] / 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     if ( $record->{ 'READ_PAIRED' } ) {
147         $record->{ 'BLOCK_COUNT' } = 2;
148     } else {
149         $record->{ 'BLOCK_COUNT' } = 1;
150     }
151
152     # extra fields - (hate hate hate SAM!!!)
153     
154     for ( $i = QUAL + 1; $i < @{ $entry }; $i++ )
155     {
156         ( $tag, $vtype, $value ) = split /:/, $entry->[ $i ];
157
158         if ( $tag eq "MD" or $tag eq "NM" ) {
159             $record->{ $tag } = $value;       
160         }
161     }
162
163     if ( defined $record->{ "NM" } and $record->{ "NM" } > 0 ) {
164         $record->{ "ALIGN" } = sam2align( $record->{ 'SEQ' }, $record->{ 'CIGAR' }, $record->{ 'NM' }, $record->{ 'MD' } );
165     } else {
166         $record->{ "ALIGN" } = ".";
167     }
168
169     return wantarray ? %{ $record } : $record;
170 }
171
172
173 sub sam2align
174 {
175     my ( $seq,
176          $cigar,
177          $nm,
178          $md,
179        ) = @_;
180
181     # Returns a string.
182
183     my ( $offset, $len, $op, $i, $nt, $nt2, @nts, $insertions, @align, $aln_str );
184
185     $offset     = 0;
186     $insertions = 0;
187
188     while ( length( $cigar ) > 0 )
189     {
190         if ( $cigar =~ s/^(\d+)(\w)// )
191         {
192             $len = $1;
193             $op  = $2;
194
195             if ( $op eq 'I' )
196             {
197                 for ( $i = 0; $i < $len; $i++ )
198                 {
199                     $nt = substr $seq, $offset + $i, 1;
200
201                     push @align, [ $offset + $i, "->$nt" ];
202
203                     $nm--;
204                 }
205             }
206
207             $offset += $len;
208         }
209     }
210     
211     $offset = 0;
212
213     while ( $nm > 0 )
214     {
215         if ( $md =~ s/^(\d+)// )   # match
216         {
217             $offset += $1;
218         }
219         elsif( $md =~ s/^\^([ATCGN]+)// )  # insertions
220         {
221             @nts = split //, $1;
222
223             foreach $nt ( @nts )
224             {
225                 $nt2 = substr $seq, $offset, 1;
226
227                 push @align, [ $offset, "$nt2>-" ];
228         
229                 $offset++;
230
231                 $insertions++;
232
233                 $nm--;
234             }
235         }
236         elsif ( $md =~ s/^([ATCGN]+)// ) # mismatch
237         {
238             @nts = split //, $1;
239
240             foreach $nt ( @nts )
241             {
242                 $nt2 = substr $seq, $offset - $insertions, 1;
243
244                 push @align, [ $offset, "$nt>$nt2" ];
245
246                 $offset++;
247
248                 $nm--;
249             }
250         }
251         else
252         {
253             Maasha::Common::error( qq(Bad SAM field -> MD: "$md") );
254         }
255     }
256
257     $aln_str = "";
258
259     @align = sort { $a->[ 0 ] <=> $b->[ 0 ] } @align;
260
261     map { $aln_str .= "$_->[ 0 ]:$_->[ 1 ]," } @align;
262
263     $aln_str =~ s/,$//;
264
265     return $aln_str;
266 }
267
268
269 sub biopiece2sam
270 {
271
272 }
273
274
275
276 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
277
278 1;