File Coverage

Bio/SearchIO/megablast.pm
Criterion Covered Total %
statement 97 129 75.1
branch 30 50 60.0
condition 4 12 33.3
subroutine 13 18 72.2
pod 11 12 91.6
total 155 221 70.1


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::SearchIO::megablast
3             #
4             # Please direct questions and support issues to
5             #
6             # Cared for by Jason Stajich
7             #
8             # Copyright Jason Stajich
9             #
10             # You may distribute this module under the same terms as perl itself
11              
12             # POD documentation - main docs before the code
13              
14             =head1 NAME
15              
16             Bio::SearchIO::megablast - a driver module for Bio::SearchIO to parse
17             megablast reports (format 0)
18              
19             =head1 SYNOPSIS
20              
21             # do not use this module directly
22              
23             use Bio::SearchIO;
24             # for default format output from megablast
25             my $in = Bio::SearchIO->new(-file => 'file.mbl',
26             -format => 'megablast',
27             -report_format => 0);
28              
29             while( my $r = $in->next_result ) {
30             while( my $hit = $r->next_hit ) {
31             while( my $hsp = $hit->next_hsp ) {
32             }
33             }
34             }
35              
36             =head1 DESCRIPTION
37              
38             Beware!
39              
40             Because of the way megablast report format 0 is coded, realize that score
41             means # gap characters + # mismatches for a HSP.
42              
43             The docs from NCBI regarding FORMAT 0
44             # 0: Produce one-line output for each alignment, in the form
45             #
46             # 'subject-id'=='[+-]query-id' (s_off q_off s_end q_end) score
47             #
48             # Here subject(query)-id is a gi number, an accession or some other type of
49             # identifier found in the FASTA definition line of the respective sequence.
50             #
51             # + or - corresponds to same or different strand alignment.
52             #
53             # Score for non-affine gapping parameters means the total number of
54             # differences (mismatches + gap characters). For affine case it is the
55             # actual (raw) score of the alignment.
56              
57             FORMAT 1 parsing has not been implemented
58             FORMAT 2 parsing should work with the SearchIO 'blast' parser
59              
60             =head1 FEEDBACK
61              
62             =head2 Mailing Lists
63              
64             User feedback is an integral part of the evolution of this and other
65             Bioperl modules. Send your comments and suggestions preferably to
66             the Bioperl mailing list. Your participation is much appreciated.
67              
68             bioperl-l@bioperl.org - General discussion
69             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
70              
71             =head2 Support
72              
73             Please direct usage questions or support issues to the mailing list:
74              
75             I
76              
77             rather than to the module maintainer directly. Many experienced and
78             reponsive experts will be able look at the problem and quickly
79             address it. Please include a thorough description of the problem
80             with code and data examples if at all possible.
81              
82             =head2 Reporting Bugs
83              
84             Report bugs to the Bioperl bug tracking system to help us keep track
85             of the bugs and their resolution. Bug reports can be submitted via the
86             web:
87              
88             https://github.com/bioperl/bioperl-live/issues
89              
90             =head1 AUTHOR - Jason Stajich
91              
92             Email jason-at-bioperl.org
93              
94             =head1 APPENDIX
95              
96             The rest of the documentation details each of the object methods.
97             Internal methods are usually preceded with a _
98              
99             =cut
100              
101              
102             # Let the code begin...
103              
104              
105             package Bio::SearchIO::megablast;
106 1     1   4 use strict;
  1         1  
  1         29  
107 1     1   3 use vars qw(%MAPPING %MODEMAP $DEFAULT_BLAST_WRITER_CLASS);
  1         1  
  1         116  
108              
109 1     1   4 use base qw(Bio::SearchIO);
  1         1  
  1         116  
110              
111             BEGIN {
112             # mapping of MegaBlast terms to Bioperl hash keys
113 1     1   4 %MODEMAP = ('MegaBlastOutput' => 'result',
114             'Hit' => 'hit',
115             'Hsp' => 'hsp'
116             );
117              
118             # This should really be done more intelligently, like with
119             # XSLT
120              
121 1         5 %MAPPING =
122             (
123             'Hsp_query-from' => 'HSP-query_start',
124             'Hsp_query-to' => 'HSP-query_end',
125             'Hsp_hit-from' => 'HSP-hit_start',
126             'Hsp_hit-to' => 'HSP-hit_end',
127             'Hit_score' => 'HIT-score',
128             'Hsp_score' => 'HSP-score',
129            
130             'Hsp_identity' => 'HSP-identical',
131             'Hsp_positive' => 'HSP-conserved',
132              
133             'Hit_id' => 'HIT-name',
134            
135             'MegaBlastOutput_program' => 'RESULT-algorithm_name',
136             'MegaBlastOutput_query-def'=> 'RESULT-query_name',
137             );
138              
139              
140 1         875 $DEFAULT_BLAST_WRITER_CLASS = 'Bio::SearchIO::Writer::HitTableWriter';
141             }
142              
143              
144             =head2 new
145              
146             Title : new
147             Usage : my $obj = Bio::SearchIO::blast->new();
148             Function: Builds a new Bio::SearchIO::blast object
149             Returns : Bio::SearchIO::blast
150             Args : -fh/-file => filehandle/filename to BLAST file
151             -format => 'blast'
152              
153             =cut
154              
155             sub _initialize {
156 1     1   1 my ($self,@args) = @_;
157 1         5 $self->SUPER::_initialize(@args);
158 1         4 my ($fmt) = $self->_rearrange([qw(REPORT_FORMAT)], @args);
159              
160 1 50       4 $self->throw("Must provide a value for -report_format when initializing a megablast parser") unless defined $fmt ;
161 1         5 $self->report_format($fmt);
162 1         2 return 1;
163             }
164              
165             =head2 next_result
166              
167             Title : next_result
168             Usage : my $hit = $searchio->next_result;
169             Function: Returns the next Result from a search
170             Returns : Bio::Search::Result::ResultI object
171             Args : none
172              
173             =cut
174              
175             sub next_result{
176 1     1 1 7 my ($self) = @_;
177            
178 1         5 local $/ = "\n";
179 1         1 local $_;
180              
181 1         3 my $fmt = $self->report_format;
182 1         2 my ($lastquery,$lasthit);
183 1         15 while( defined($_ = $self->_readline) ) {
184 21 50       33 if( $fmt == 0 ) {
185 21 50       107 if( /^\'(\S+)\'\=\=\'(\+|\-)(\S+)\'\s+
186             \((\d+)\s+(\d+)\s+(\d+)\s+(\d+)\)\s+
187             (\d+)/ox )
188             {
189 21         66 my ($hit,$strand,$query,
190             $h_start,$q_start,$h_end,$q_end,
191             $score) = ($1,$2,$3,$4,$5,$6,$7,$8);
192 21 100       44 if( ! defined $lastquery ) {
    50          
193 1         5 $self->start_element({'Name' => 'MegaBlastOutput'});
194 1         7 $self->element({'Name' => 'MegaBlastOutput_program',
195             'Data' => 'MEGABLAST'});
196 1         3 $self->element({'Name' => 'MegaBlastOutput_query-def',
197             'Data' => $query});
198             } elsif( $lastquery ne $query ) {
199 0         0 $self->_pushback($_);
200 0 0       0 $self->end_element({'Name' => 'Hit'}) if( defined $lasthit);
201 0         0 $self->end_element({ 'Name' => 'MegaBlastOutput'});
202 0         0 $lasthit = undef;
203 0         0 $lastquery = undef;
204 0         0 return $self->end_document();
205             }
206              
207 21 100 100     62 if( ! defined $lasthit || $lasthit ne $hit ) {
208 4 100       13 $self->end_element({'Name' => 'Hit'}) if( defined $lasthit);
209 4         12 $self->start_element({'Name' => 'Hit'});
210 4         14 $self->element({'Name' => 'Hit_id',
211             'Data' => $hit});
212             }
213 21         43 $self->start_element({'Name' => 'Hsp'});
214 21         54 $self->element({'Name' => 'Hsp_score',
215             'Data' => $score});
216              
217             # flip flop start/end if strand is < 0
218             # since strandedness is inferred from the query
219             # because of the way it is coded all queries will
220             # be on the forward strand and hits will be either
221             # +/-
222              
223             # also the NCBI docs state:
224             # 0: Produce one-line output for each alignment, in the form
225             #
226             # 'subject-id'=='[+-]query-id' (s_off q_off s_end q_end) score
227             #
228             # Here subject(query)-id is a gi number, an accession or some other type of
229             # identifier found in the FASTA definition line of the respective sequence.
230             #
231             # + or - corresponds to same or different strand alignment.
232             #
233             # Score for non-affine gapping parameters means the total number of
234             # differences (mismatches + gap characters). For affine case it is the
235             # actual (raw) score of the alignment.
236            
237             # and yet when rev strand hits are made I see
238             # (MBL 2.2.4)
239             # 'Contig634'=='-503384' (1 7941 321 7620) 19
240             # so the query is on the rev strand and the
241             # subject is on the fwd strand
242             # so I am flip-flopping everything when I see a '-'
243 21 100       44 if( $strand eq '-' ) {
244 18         24 ($h_start,$h_end) = ( $h_end,$h_start);
245 18         23 ($q_start,$q_end) = ( $q_end,$q_start);
246             }
247 21         40 $self->element({'Name' => 'Hsp_hit-from',
248             'Data' => $h_start});
249 21         48 $self->element({'Name' => 'Hsp_hit-to',
250             'Data' => $h_end});
251 21         47 $self->element({'Name' => 'Hsp_query-from',
252             'Data' => $q_start});
253 21         47 $self->element({'Name' => 'Hsp_query-to',
254             'Data' => $q_end});
255             # might not be quite right -- need to know length of the HSP
256 21         43 my $numid = (abs($q_end - $q_start) - $score);
257              
258 21         39 $self->element({'Name' => 'Hsp_identity',
259             'Data' => $numid});
260 21         38 $self->element({'Name' => 'Hsp_positive',
261             'Data' => $numid});
262              
263 21         38 $self->end_element({'Name' => 'Hsp'});
264 21         26 $lasthit = $hit;
265 21         55 $lastquery = $query;
266             } else {
267 0         0 $self->debug("Unknown line in fmt0 parsing: $_");
268             }
269             }
270             }
271 1 50 33     6 if( defined $lastquery && $fmt == 0 ) {
272 1 50       18 $self->end_element({'Name' => 'Hit'}) if( defined $lasthit);
273 1         4 $self->end_element({ 'Name' => 'MegaBlastOutput'});
274 1         4 return $self->end_document();
275             }
276 0         0 return;
277             }
278              
279             =head2 report_format
280              
281             Title : report_format
282             Usage : $obj->report_format($newval)
283             Function: Get/Set the report_format value
284             Returns : value of report_format (a scalar)
285             Args : on set, new value (a scalar or undef, optional)
286              
287              
288             =cut
289              
290             sub report_format{
291 2     2 1 2 my $self = shift;
292 2 100       6 return $self->{'_report_format'} = shift if @_;
293 1         2 return $self->{'_report_format'};
294             }
295              
296              
297             =head2 start_element
298              
299             Title : start_element
300             Usage : $eventgenerator->start_element
301             Function: Handles a start element event
302             Returns : none
303             Args : hashref with at least 2 keys 'Data' and 'Name'
304              
305              
306             =cut
307              
308             sub start_element{
309 179     179 1 104 my ($self,$data) = @_;
310             # we currently do not care about attributes
311 179         142 my $nm = $data->{'Name'};
312 179 100       222 if( my $type = $MODEMAP{$nm} ) {
313 26         36 $self->_mode($type);
314 26 50       46 if( $self->_eventHandler->will_handle($type) ) {
315 26         65 my $func = sprintf("start_%s",lc $type);
316 26         36 $self->_eventHandler->$func($data->{'Attributes'});
317             }
318 26         32 unshift @{$self->{'_elements'}}, $type;
  26         43  
319             }
320              
321 179 100       252 if($nm eq 'MegaBlastOutput') {
322 1         2 $self->{'_values'} = {};
323 1         2 $self->{'_result'}= undef;
324 1         2 $self->{'_mode'} = '';
325             }
326              
327             }
328              
329             =head2 end_element
330              
331             Title : start_element
332             Usage : $eventgenerator->end_element
333             Function: Handles an end element event
334             Returns : none
335             Args : hashref with at least 2 keys 'Data' and 'Name'
336              
337              
338             =cut
339              
340             sub end_element {
341 179     179 1 103 my ($self,$data) = @_;
342 179         126 my $nm = $data->{'Name'};
343 179         103 my $rc;
344              
345 179 100       315 if( my $type = $MODEMAP{$nm} ) {
    50          
346 26 50       38 if( $self->_eventHandler->will_handle($type) ) {
347 26         48 my $func = sprintf("end_%s",lc $type);
348             $rc = $self->_eventHandler->$func($self->{'_reporttype'},
349 26         35 $self->{'_values'});
350             }
351 26         22 shift @{$self->{'_elements'}};
  26         30  
352              
353             } elsif( $MAPPING{$nm} ) {
354 153 50       166 if ( ref($MAPPING{$nm}) =~ /hash/i ) {
355 0         0 my $key = (keys %{$MAPPING{$nm}})[0];
  0         0  
356 0         0 $self->{'_values'}->{$key}->{$MAPPING{$nm}->{$key}} = $self->{'_last_data'};
357             } else {
358 153         181 $self->{'_values'}->{$MAPPING{$nm}} = $self->{'_last_data'};
359             }
360             } else {
361 0         0 $self->warn( "unknown nm $nm ignoring\n");
362             }
363 179         133 $self->{'_last_data'} = ''; # remove read data if we are at
364             # end of an element
365 179 100       215 $self->{'_result'} = $rc if( $nm eq 'MegaBlastOutput' );
366 179         155 return $rc;
367              
368             }
369              
370             =head2 element
371              
372             Title : element
373             Usage : $eventhandler->element({'Name' => $name, 'Data' => $str});
374             Function: Convience method that calls start_element, characters, end_element
375             Returns : none
376             Args : Hash ref with the keys 'Name' and 'Data'
377              
378              
379             =cut
380              
381             sub element{
382 153     153 1 111 my ($self,$data) = @_;
383 153         142 $self->start_element($data);
384 153         155 $self->characters($data);
385 153         147 $self->end_element($data);
386             }
387              
388              
389             =head2 characters
390              
391             Title : characters
392             Usage : $eventgenerator->characters($str)
393             Function: Send a character events
394             Returns : none
395             Args : string
396              
397              
398             =cut
399              
400             sub characters{
401 153     153 1 94 my ($self,$data) = @_;
402 153 50       177 return unless defined $data->{'Data'};
403 153         139 $self->{'_last_data'} = $data->{'Data'};
404             }
405              
406             =head2 _mode
407              
408             Title : _mode
409             Usage : $obj->_mode($newval)
410             Function:
411             Example :
412             Returns : value of _mode
413             Args : newvalue (optional)
414              
415              
416             =cut
417              
418             sub _mode{
419 26     26   16 my ($self,$value) = @_;
420 26 50       35 if( defined $value) {
421 26         25 $self->{'_mode'} = $value;
422             }
423 26         22 return $self->{'_mode'};
424             }
425              
426             =head2 within_element
427              
428             Title : within_element
429             Usage : if( $eventgenerator->within_element($element) ) {}
430             Function: Test if we are within a particular element
431             This is different than 'in' because within can be tested
432             for a whole block.
433             Returns : boolean
434             Args : string element name
435              
436              
437             =cut
438              
439             sub within_element{
440 0     0 1 0 my ($self,$name) = @_;
441             return 0 if ( ! defined $name &&
442             ! defined $self->{'_elements'} ||
443 0 0 0     0 scalar @{$self->{'_elements'}} == 0) ;
  0   0     0  
444 0         0 foreach ( @{$self->{'_elements'}} ) {
  0         0  
445 0 0       0 if( $_ eq $name ) {
446 0         0 return 1;
447             }
448             }
449 0         0 return 0;
450             }
451              
452             =head2 in_element
453              
454             Title : in_element
455             Usage : if( $eventgenerator->in_element($element) ) {}
456             Function: Test if we are in a particular element
457             This is different than 'in' because within can be tested
458             for a whole block.
459             Returns : boolean
460             Args : string element name
461              
462              
463             =cut
464              
465             sub in_element{
466 0     0 1 0 my ($self,$name) = @_;
467 0 0       0 return 0 if ! defined $self->{'_elements'}->[0];
468 0         0 return ( $self->{'_elements'}->[0] eq $name)
469             }
470              
471              
472             =head2 start_document
473              
474             Title : start_document
475             Usage : $eventgenerator->start_document
476             Function: Handles a start document event
477             Returns : none
478             Args : none
479              
480              
481             =cut
482              
483             sub start_document{
484 0     0 1 0 my ($self) = @_;
485 0         0 $self->{'_lasttype'} = '';
486 0         0 $self->{'_values'} = {};
487 0         0 $self->{'_result'}= undef;
488 0         0 $self->{'_mode'} = '';
489 0         0 $self->{'_elements'} = [];
490             }
491              
492              
493             =head2 end_document
494              
495             Title : end_document
496             Usage : $eventgenerator->end_document
497             Function: Handles an end document event
498             Returns : Bio::Search::Result::ResultI object
499             Args : none
500              
501              
502             =cut
503              
504             sub end_document{
505 1     1 1 2 my ($self,@args) = @_;
506 1         6 return $self->{'_result'};
507             }
508              
509             =head2 result_count
510              
511             Title : result_count
512             Usage : my $count = $searchio->result_count
513             Function: Returns the number of results we have processed
514             Returns : integer
515             Args : none
516              
517              
518             =cut
519              
520             sub result_count {
521 0     0 1   my $self = shift;
522 0           return $self->{'_result_count'};
523             }
524              
525 0     0 0   sub report_count { shift->result_count }
526              
527             1;