File Coverage

Bio/SearchIO/waba.pm
Criterion Covered Total %
statement 137 167 82.0
branch 40 64 62.5
condition 8 15 53.3
subroutine 17 20 85.0
pod 10 11 90.9
total 212 277 76.5


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::SearchIO::waba
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::waba - SearchIO parser for Jim Kent WABA program
17             alignment output
18              
19             =head1 SYNOPSIS
20              
21             # do not use this object directly, rather through Bio::SearchIO
22              
23             use Bio::SearchIO;
24             my $in = Bio::SearchIO->new(-format => 'waba',
25             -file => 'output.wab');
26             while( my $result = $in->next_result ) {
27             while( my $hit = $result->next_hit ) {
28             while( my $hsp = $result->next_hsp ) {
29              
30             }
31             }
32             }
33              
34             =head1 DESCRIPTION
35              
36             This parser will process the waba output (NOT the human readable format).
37              
38             =head1 FEEDBACK
39              
40             =head2 Mailing Lists
41              
42             User feedback is an integral part of the evolution of this and other
43             Bioperl modules. Send your comments and suggestions preferably to
44             the Bioperl mailing list. Your participation is much appreciated.
45              
46             bioperl-l@bioperl.org - General discussion
47             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
48              
49             =head2 Support
50              
51             Please direct usage questions or support issues to the mailing list:
52              
53             I
54              
55             rather than to the module maintainer directly. Many experienced and
56             reponsive experts will be able look at the problem and quickly
57             address it. Please include a thorough description of the problem
58             with code and data examples if at all possible.
59              
60             =head2 Reporting Bugs
61              
62             Report bugs to the Bioperl bug tracking system to help us keep track
63             of the bugs and their resolution. Bug reports can be submitted via the
64             web:
65              
66             https://github.com/bioperl/bioperl-live/issues
67              
68             =head1 AUTHOR - Jason Stajich
69              
70             Email jason-at-bioperl.org
71              
72             =head1 APPENDIX
73              
74             The rest of the documentation details each of the object methods.
75             Internal methods are usually preceded with a _
76              
77             =cut
78              
79              
80             # Let the code begin...
81              
82              
83             package Bio::SearchIO::waba;
84 1     1   5 use vars qw(%MODEMAP %MAPPING @STATES);
  1         2  
  1         51  
85 1     1   5 use strict;
  1         1  
  1         19  
86              
87             # Object preamble - inherits from Bio::Root::Root
88              
89 1     1   262 use Bio::Search::Result::ResultFactory;
  1         2  
  1         21  
90 1     1   225 use Bio::Search::HSP::HSPFactory;
  1         2  
  1         22  
91              
92 1     1   243 use POSIX;
  1         4104  
  1         4  
93              
94             BEGIN {
95             # mapping of NCBI Blast terms to Bioperl hash keys
96 1     1   2415 %MODEMAP = ('WABAOutput' => 'result',
97             'Hit' => 'hit',
98             'Hsp' => 'hsp'
99             );
100 1         2 @STATES = qw(Hsp_qseq Hsp_hseq Hsp_stateseq);
101 1         27 %MAPPING =
102             (
103             'Hsp_query-from'=> 'HSP-query_start',
104             'Hsp_query-to' => 'HSP-query_end',
105             'Hsp_hit-from' => 'HSP-hit_start',
106             'Hsp_hit-to' => 'HSP-hit_end',
107             'Hsp_qseq' => 'HSP-query_seq',
108             'Hsp_hseq' => 'HSP-hit_seq',
109             'Hsp_midline' => 'HSP-homology_seq',
110             'Hsp_stateseq' => 'HSP-hmmstate_seq',
111             'Hsp_align-len' => 'HSP-hsp_length',
112            
113             'Hit_id' => 'HIT-name',
114             'Hit_accession' => 'HIT-accession',
115              
116             'WABAOutput_program' => 'RESULT-algorithm_name',
117             'WABAOutput_version' => 'RESULT-algorithm_version',
118             'WABAOutput_query-def'=> 'RESULT-query_name',
119             'WABAOutput_query-db' => 'RESULT-query_database',
120             'WABAOutput_db' => 'RESULT-database_name',
121             );
122             }
123              
124              
125 1     1   6 use base qw(Bio::SearchIO);
  1         2  
  1         1312  
126              
127             =head2 new
128              
129             Title : new
130             Usage : my $obj = Bio::SearchIO::waba->new();
131             Function: Builds a new Bio::SearchIO::waba object
132             Returns : Bio::SearchIO::waba
133             Args : see Bio::SearchIO
134              
135             =cut
136              
137             sub _initialize {
138 1     1   3 my ($self,@args) = @_;
139 1         4 $self->SUPER::_initialize(@args);
140 1         5 $self->_eventHandler->register_factory('result', Bio::Search::Result::ResultFactory->new(-type => 'Bio::Search::Result::WABAResult'));
141              
142 1         4 $self->_eventHandler->register_factory('hsp', Bio::Search::HSP::HSPFactory->new(-type => 'Bio::Search::HSP::WABAHSP'));
143             }
144              
145              
146             =head2 next_result
147              
148             Title : next_result
149             Usage : my $hit = $searchio->next_result;
150             Function: Returns the next Result from a search
151             Returns : Bio::Search::Result::ResultI object
152             Args : none
153              
154             =cut
155              
156             sub next_result{
157 3     3 1 533 my ($self) = @_;
158 3         12 local $/ = "\n";
159 3         4 local $_;
160            
161 3         6 my ($curquery,$curhit);
162 3         4 my $state = -1;
163 3         11 $self->start_document();
164 3         4 my @hit_signifs;
165 3         18 while( defined ($_ = $self->_readline )) {
166            
167 17 100       31 if( $state == -1 ) {
    50          
168 5         58 my ($qid, $qhspid,$qpercent, $junk,
169             $alnlen,$qdb,$qacc,$qstart,$qend,$qstrand,
170             $hitdb,$hacc,$hstart,$hend,
171             $hstrand) =
172             ( /^(\S+)\.(\S+)\s+align\s+ # get the queryid
173             (\d+(\.\d+)?)\%\s+ # get the percentage
174             of\s+(\d+)\s+ # get the length of the alignment
175             (\S+)\s+ # this is the query database
176             (\S+):(\-?\d+)\-(\-?\d+) # The accession:start-end for query
177             \s+([\-\+]) # query strand
178             \s+(\S+)\. # hit db
179             (\S+):(\-?\d+)\-(\-?\d+) # The accession:start-end for hit
180             \s+([\-\+])\s*$ # hit strand
181             /ox );
182            
183             # Curses. Jim's code is 0 based, the following is to readjust
184 5 50       14 if( $hstart < 0 ) { $hstart *= -1}
  0         0  
185 5 50       10 if( $hend < 0 ) { $hend *= -1}
  0         0  
186 5 50       8 if( $qstart < 0 ) { $qstart *= -1}
  0         0  
187 5 50       8 if( $qend < 0 ) { $qend *= -1}
  0         0  
188 5         6 $hstart++; $hend++; $qstart++; $qend++;
  5         6  
  5         5  
  5         4  
189 5 50       8 if( ! defined $alnlen ) {
190 0         0 $self->warn("Unable to parse the rest of the WABA alignment info for: '$_'");
191 0         0 last;
192             }
193 5         8 $self->{'_reporttype'} = 'WABA'; # hardcoded - only
194             # one type of WABA AFAIK
195 5 100 100     18 if( defined $curquery &&
196             $curquery ne $qid ) {
197 1         4 $self->end_element({'Name' => 'Hit'});
198 1         8 $self->_pushback($_);
199 1         3 $self->end_element({'Name' => 'WABAOutput'});
200 1         4 return $self->end_document();
201             }
202            
203 4 50 66     14 if( defined $curhit &&
    100          
204             $curhit ne $hacc) {
205             # slight duplication here -- keep these in SYNC
206 0         0 $self->end_element({'Name' => 'Hit'});
207 0         0 $self->start_element({'Name' => 'Hit'});
208 0         0 $self->element({'Name' => 'Hit_id',
209             'Data' => $hacc});
210 0         0 $self->element({'Name' => 'Hit_accession',
211             'Data' => $hacc});
212              
213             } elsif ( ! defined $curquery ) {
214 2         10 $self->start_element({'Name' => 'WABAOutput'});
215 2         5 $self->{'_result_count'}++;
216 2         14 $self->element({'Name' => 'WABAOutput_query-def',
217             'Data' => $qid });
218 2         8 $self->element({'Name' => 'WABAOutput_program',
219             'Data' => 'WABA'});
220 2         6 $self->element({'Name' => 'WABAOutput_query-db',
221             'Data' => $qdb});
222 2         6 $self->element({'Name' => 'WABAOutput_db',
223             'Data' => $hitdb});
224            
225             # slight duplication here -- keep these N'SYNC ;-)
226 2         6 $self->start_element({'Name' => 'Hit'});
227 2         6 $self->element({'Name' => 'Hit_id',
228             'Data' => $hacc});
229 2         6 $self->element({'Name' => 'Hit_accession',
230             'Data' => $hacc});
231             }
232              
233            
234             # strand is inferred by start,end values
235             # in the Result Builder
236 4 50       8 if( $qstrand eq '-' ) {
237 0         0 ($qstart,$qend) = ($qend,$qstart);
238             }
239 4 50       8 if( $hstrand eq '-' ) {
240 0         0 ($hstart,$hend) = ($hend,$hstart);
241             }
242              
243 4         14 $self->start_element({'Name' => 'Hsp'});
244 4         13 $self->element({'Name' => 'Hsp_query-from',
245             'Data' => $qstart});
246 4         10 $self->element({'Name' => 'Hsp_query-to',
247             'Data' => $qend});
248 4         12 $self->element({'Name' => 'Hsp_hit-from',
249             'Data' => $hstart});
250 4         12 $self->element({'Name' => 'Hsp_hit-to',
251             'Data' => $hend});
252 4         11 $self->element({'Name' => 'Hsp_align-len',
253             'Data' => $alnlen});
254            
255 4         7 $curquery = $qid;
256 4         4 $curhit = $hacc;
257 4         12 $state = 0;
258             } elsif( ! defined $curquery ) {
259 0 0       0 $self->warn("skipping because no Hit begin line was recognized\n$_") if( $_ !~ /^\s+$/ );
260 0         0 next;
261             } else {
262 12         19 chomp;
263 12         38 $self->element({'Name' => $STATES[$state++],
264             'Data' => $_});
265 12 100       32 if( $state >= scalar @STATES ) {
266 4         5 $state = -1;
267 4         8 $self->end_element({'Name' => 'Hsp'});
268             }
269             }
270             }
271 2 100       7 if( defined $curquery ) {
272 1         4 $self->end_element({'Name' => 'Hit'});
273 1         4 $self->end_element({'Name' => 'WABAOutput'});
274 1         3 return $self->end_document();
275             }
276 1         6 return;
277             }
278              
279             =head2 start_element
280              
281             Title : start_element
282             Usage : $eventgenerator->start_element
283             Function: Handles a start element event
284             Returns : none
285             Args : hashref with at least 2 keys 'Data' and 'Name'
286              
287              
288             =cut
289              
290             sub start_element{
291 64     64 1 65 my ($self,$data) = @_;
292             # we currently don't care about attributes
293 64         66 my $nm = $data->{'Name'};
294 64 100       98 if( my $type = $MODEMAP{$nm} ) {
295 8         15 $self->_mode($type);
296 8 50       17 if( $self->_eventHandler->will_handle($type) ) {
297 8         26 my $func = sprintf("start_%s",lc $type);
298 8         13 $self->_eventHandler->$func($data->{'Attributes'});
299             }
300 8         13 unshift @{$self->{'_elements'}}, $type;
  8         16  
301             }
302 64 100       96 if($nm eq 'WABAOutput') {
303 2         4 $self->{'_values'} = {};
304 2         3 $self->{'_result'}= undef;
305 2         4 $self->{'_mode'} = '';
306             }
307              
308             }
309              
310             =head2 end_element
311              
312             Title : start_element
313             Usage : $eventgenerator->end_element
314             Function: Handles an end element event
315             Returns : none
316             Args : hashref with at least 2 keys 'Data' and 'Name'
317              
318              
319             =cut
320              
321             sub end_element {
322 64     64 1 65 my ($self,$data) = @_;
323 64         64 my $nm = $data->{'Name'};
324 64         51 my $rc;
325             # Hsp are sort of weird, in that they end when another
326             # object begins so have to detect this in end_element for now
327 64 100       77 if( $nm eq 'Hsp' ) {
328 4         8 foreach ( qw(Hsp_qseq Hsp_midline Hsp_hseq) ) {
329             $self->element({'Name' => $_,
330 12         25 'Data' => $self->{'_last_hspdata'}->{$_}});
331             }
332 4         19 $self->{'_last_hspdata'} = {}
333             }
334              
335 64 100       121 if( my $type = $MODEMAP{$nm} ) {
    50          
336 8 50       19 if( $self->_eventHandler->will_handle($type) ) {
337 8         21 my $func = sprintf("end_%s",lc $type);
338             $rc = $self->_eventHandler->$func($self->{'_reporttype'},
339 8         15 $self->{'_values'});
340             }
341 8         11 shift @{$self->{'_elements'}};
  8         13  
342              
343             } elsif( $MAPPING{$nm} ) {
344 56 50       71 if ( ref($MAPPING{$nm}) =~ /hash/i ) {
345 0         0 my $key = (keys %{$MAPPING{$nm}})[0];
  0         0  
346 0         0 $self->{'_values'}->{$key}->{$MAPPING{$nm}->{$key}} = $self->{'_last_data'};
347             } else {
348 56         97 $self->{'_values'}->{$MAPPING{$nm}} = $self->{'_last_data'};
349             }
350             } else {
351 0         0 $self->warn( "unknown nm $nm ignoring\n");
352             }
353 64         67 $self->{'_last_data'} = ''; # remove read data if we are at
354             # end of an element
355 64 100       80 $self->{'_result'} = $rc if( $nm eq 'WABAOutput' );
356 64         87 return $rc;
357              
358             }
359              
360             =head2 element
361              
362             Title : element
363             Usage : $eventhandler->element({'Name' => $name, 'Data' => $str});
364             Function: Convience method that calls start_element, characters, end_element
365             Returns : none
366             Args : Hash ref with the keys 'Name' and 'Data'
367              
368              
369             =cut
370              
371             sub element{
372 56     56 1 64 my ($self,$data) = @_;
373 56         81 $self->start_element($data);
374 56         80 $self->characters($data);
375 56         81 $self->end_element($data);
376             }
377              
378              
379             =head2 characters
380              
381             Title : characters
382             Usage : $eventgenerator->characters($str)
383             Function: Send a character events
384             Returns : none
385             Args : string
386              
387              
388             =cut
389              
390             sub characters{
391 56     56 1 64 my ($self,$data) = @_;
392              
393 56 100       78 return unless ( defined $data->{'Data'} );
394 52 50       298 if( $data->{'Data'} =~ /^\s+$/ ) {
395 0 0       0 return unless $data->{'Name'} =~ /Hsp\_(midline|qseq|hseq)/;
396             }
397              
398 52 100 100     68 if( $self->in_element('hsp') &&
399             $data->{'Name'} =~ /Hsp\_(qseq|hseq|midline)/ ) {
400            
401 16         113 $self->{'_last_hspdata'}->{$data->{'Name'}} .= $data->{'Data'};
402             }
403            
404 52         82 $self->{'_last_data'} = $data->{'Data'};
405             }
406              
407             =head2 _mode
408              
409             Title : _mode
410             Usage : $obj->_mode($newval)
411             Function:
412             Example :
413             Returns : value of _mode
414             Args : newvalue (optional)
415              
416              
417             =cut
418              
419             sub _mode{
420 8     8   11 my ($self,$value) = @_;
421 8 50       14 if( defined $value) {
422 8         8 $self->{'_mode'} = $value;
423             }
424 8         12 return $self->{'_mode'};
425             }
426              
427             =head2 within_element
428              
429             Title : within_element
430             Usage : if( $eventgenerator->within_element($element) ) {}
431             Function: Test if we are within a particular element
432             This is different than 'in' because within can be tested
433             for a whole block.
434             Returns : boolean
435             Args : string element name
436              
437              
438             =cut
439              
440             sub within_element{
441 0     0 1 0 my ($self,$name) = @_;
442             return 0 if ( ! defined $name &&
443             ! defined $self->{'_elements'} ||
444 0 0 0     0 scalar @{$self->{'_elements'}} == 0) ;
  0   0     0  
445 0         0 foreach ( @{$self->{'_elements'}} ) {
  0         0  
446 0 0       0 if( $_ eq $name ) {
447 0         0 return 1;
448             }
449             }
450 0         0 return 0;
451             }
452              
453             =head2 in_element
454              
455             Title : in_element
456             Usage : if( $eventgenerator->in_element($element) ) {}
457             Function: Test if we are in a particular element
458             This is different than 'in' because within can be tested
459             for a whole block.
460             Returns : boolean
461             Args : string element name
462              
463              
464             =cut
465              
466             sub in_element{
467 52     52 1 57 my ($self,$name) = @_;
468 52 50       71 return 0 if ! defined $self->{'_elements'}->[0];
469 52         186 return ( $self->{'_elements'}->[0] eq $name)
470             }
471              
472              
473             =head2 start_document
474              
475             Title : start_document
476             Usage : $eventgenerator->start_document
477             Function: Handles a start document event
478             Returns : none
479             Args : none
480              
481              
482             =cut
483              
484             sub start_document{
485 3     3 1 6 my ($self) = @_;
486 3         7 $self->{'_lasttype'} = '';
487 3         18 $self->{'_values'} = {};
488 3         13 $self->{'_result'}= undef;
489 3         7 $self->{'_mode'} = '';
490 3         7 $self->{'_elements'} = [];
491             }
492              
493              
494             =head2 end_document
495              
496             Title : end_document
497             Usage : $eventgenerator->end_document
498             Function: Handles an end document event
499             Returns : Bio::Search::Result::ResultI object
500             Args : none
501              
502              
503             =cut
504              
505             sub end_document{
506 2     2 1 4 my ($self,@args) = @_;
507 2         10 return $self->{'_result'};
508             }
509              
510             =head2 result_count
511              
512             Title : result_count
513             Usage : my $count = $searchio->result_count
514             Function: Returns the number of results we have processed
515             Returns : integer
516             Args : none
517              
518              
519             =cut
520              
521             sub result_count {
522 0     0 1   my $self = shift;
523 0           return $self->{'_result_count'};
524             }
525              
526 0     0 0   sub report_count { shift->result_count }
527              
528             1;