File Coverage

Bio/SearchIO/hmmer2.pm
Criterion Covered Total %
statement 300 326 92.0
branch 155 188 82.4
condition 65 78 83.3
subroutine 14 15 93.3
pod 10 10 100.0
total 544 617 88.1


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::SearchIO::hmmer2
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::hmmer2 - A parser for HMMER output (hmmpfam, hmmsearch)
17              
18             =head1 SYNOPSIS
19              
20             # do not use this class directly it is available through Bio::SearchIO
21             use Bio::SearchIO;
22             my $in = Bio::SearchIO->new(-format => 'hmmer2',
23             -file => 't/data/L77119.hmmer');
24             while( my $result = $in->next_result ) {
25             # this is a Bio::Search::Result::HMMERResult object
26             print $result->query_name(), " for HMM ", $result->hmm_name(), "\n";
27             while( my $hit = $result->next_hit ) {
28             print $hit->name(), "\n";
29             while( my $hsp = $hit->next_hsp ) {
30             print "length is ", $hsp->length(), "\n";
31             }
32             }
33             }
34              
35             =head1 DESCRIPTION
36              
37             This object implements a parser for HMMER output.
38              
39             =head1 FEEDBACK
40              
41             =head2 Mailing Lists
42              
43             User feedback is an integral part of the evolution of this and other
44             Bioperl modules. Send your comments and suggestions preferably to
45             the Bioperl mailing list. Your participation is much appreciated.
46              
47             bioperl-l@bioperl.org - General discussion
48             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
49              
50             =head2 Support
51              
52             Please direct usage questions or support issues to the mailing list:
53              
54             I
55              
56             rather than to the module maintainer directly. Many experienced and
57             reponsive experts will be able look at the problem and quickly
58             address it. Please include a thorough description of the problem
59             with code and data examples if at all possible.
60              
61             =head2 Reporting Bugs
62              
63             Report bugs to the Bioperl bug tracking system to help us keep track
64             of the bugs and their resolution. Bug reports can be submitted via the
65             web:
66              
67             https://github.com/bioperl/bioperl-live/issues
68              
69             =head1 AUTHOR - Jason Stajich
70              
71             Email jason-at-bioperl.org
72              
73             =head1 APPENDIX
74              
75             The rest of the documentation details each of the object methods.
76             Internal methods are usually preceded with a _
77              
78             =cut
79              
80             # Let the code begin...
81              
82             package Bio::SearchIO::hmmer2;
83              
84 1     1   5 use strict;
  1         2  
  1         25  
85              
86 1     1   5 use Bio::Factory::ObjectFactory;
  1         2  
  1         23  
87              
88 1         40 use vars qw(%MAPPING %MODEMAP
89 1     1   5 );
  1         1  
90              
91 1     1   4 use base qw(Bio::SearchIO::hmmer);
  1         2  
  1         132  
92              
93             BEGIN {
94              
95             # mapping of HMMER items to Bioperl hash keys
96 1     1   8 %MODEMAP = (
97             'HMMER_Output' => 'result',
98             'Hit' => 'hit',
99             'Hsp' => 'hsp'
100             );
101              
102 1         3133 %MAPPING = (
103             'Hsp_bit-score' => 'HSP-bits',
104             'Hsp_score' => 'HSP-score',
105             'Hsp_evalue' => 'HSP-evalue',
106             'Hsp_query-from' => 'HSP-query_start',
107             'Hsp_query-to' => 'HSP-query_end',
108             'Hsp_hit-from' => 'HSP-hit_start',
109             'Hsp_hit-to' => 'HSP-hit_end',
110             'Hsp_positive' => 'HSP-conserved',
111             'Hsp_identity' => 'HSP-identical',
112             'Hsp_gaps' => 'HSP-hsp_gaps',
113             'Hsp_hitgaps' => 'HSP-hit_gaps',
114             'Hsp_querygaps' => 'HSP-query_gaps',
115             'Hsp_qseq' => 'HSP-query_seq',
116             'Hsp_csline' => 'HSP-cs_seq',
117             'Hsp_hseq' => 'HSP-hit_seq',
118             'Hsp_midline' => 'HSP-homology_seq',
119             'Hsp_align-len' => 'HSP-hsp_length',
120             'Hsp_query-frame' => 'HSP-query_frame',
121             'Hsp_hit-frame' => 'HSP-hit_frame',
122              
123             'Hit_id' => 'HIT-name',
124             'Hit_len' => 'HIT-length',
125             'Hit_accession' => 'HIT-accession',
126             'Hit_desc' => 'HIT-description',
127             'Hit_signif' => 'HIT-significance',
128             'Hit_score' => 'HIT-score',
129              
130             'HMMER_program' => 'RESULT-algorithm_name',
131             'HMMER_version' => 'RESULT-algorithm_version',
132             'HMMER_query-def' => 'RESULT-query_name',
133             'HMMER_query-len' => 'RESULT-query_length',
134             'HMMER_query-acc' => 'RESULT-query_accession',
135             'HMMER_querydesc' => 'RESULT-query_description',
136             'HMMER_hmm' => 'RESULT-hmm_name',
137             'HMMER_seqfile' => 'RESULT-sequence_file',
138             'HMMER_db' => 'RESULT-database_name',
139             );
140             }
141              
142             =head2 next_result
143              
144             Title : next_result
145             Usage : my $hit = $searchio->next_result;
146             Function: Returns the next Result from a search
147             Returns : Bio::Search::Result::ResultI object
148             Args : none
149              
150             =cut
151              
152             sub next_result {
153 19     19 1 4056 my ($self) = @_;
154 19         30 my $seentop = 0;
155 19         23 my $reporttype;
156 19         29 my ( $buffer, $last, @hitinfo, @hspinfo, %hspinfo, %hitinfo );
157 19         86 local $/ = "\n";
158              
159 19         57 my $verbose = $self->verbose; # cache for speed?
160 19         66 $self->start_document();
161              
162 19         55 while ( defined( $buffer = $self->_readline ) ) {
163 461         525 my $lineorig = $buffer;
164              
165 461         461 chomp $buffer;
166 461 100 100     2242 if ($buffer =~ /^HMMER\s+(\S+)\s+\((.+)\)/o) {
    100 33        
    100          
    100          
    100          
    100          
    100          
167 16         87 my ( $prog, $version ) = split( /\s+/, $buffer );
168 16 50       42 if ($seentop) {
169 0         0 $self->_pushback($buffer);
170 0         0 $self->end_element( { 'Name' => 'HMMER_Output' } );
171 0         0 return $self->end_document();
172             }
173 16         34 $self->{'_hmmidline'} = $buffer;
174 16         71 $self->start_element( { 'Name' => 'HMMER_Output' } );
175 16         36 $self->{'_result_count'}++;
176 16         22 $seentop = 1;
177 16 50       30 if ( defined $last ) {
178 16         55 ($reporttype) = split( /\s+/, $last );
179 16 100       40 $reporttype = uc($reporttype) if defined $reporttype;
180 16         67 $self->element(
181             {
182             'Name' => 'HMMER_program',
183             'Data' => $reporttype
184             }
185             );
186             }
187             $self->element(
188             {
189 16         58 'Name' => 'HMMER_version',
190             'Data' => $version
191             }
192             );
193             }
194             elsif ($buffer =~ s/^HMM file:\s+//o) {
195 16         42 $self->{'_hmmfileline'} = $lineorig;
196 16         55 $self->element(
197             {
198             'Name' => 'HMMER_hmm',
199             'Data' => $buffer
200             }
201             );
202             }
203             elsif ($buffer =~ s/^Sequence\s+(file|database):\s+//o) {
204 16         35 $self->{'_hmmseqline'} = $lineorig;
205 16 100       50 if ( $1 eq 'database' ) {
206 5         17 $self->element(
207             {
208             'Name' => 'HMMER_db',
209             'Data' => $buffer
210             }
211             );
212             }
213             $self->element(
214             {
215 16         52 'Name' => 'HMMER_seqfile',
216             'Data' => $buffer
217             }
218             );
219             }
220             elsif ($buffer =~ s/^Query(?:\s+(?:sequence|HMM))?(?:\s+\d+)?:\s+//o) {
221 17 100       36 if ( !$seentop ) {
222              
223             # we're in a multi-query report
224 1         3 $self->_pushback($lineorig);
225 1         6 $self->_pushback( $self->{'_hmmseqline'} );
226 1         5 $self->_pushback( $self->{'_hmmfileline'} );
227 1         4 $self->_pushback( $self->{'_hmmidline'} );
228 1         4 next;
229             }
230 16         42 $buffer =~ s/\s+$//;
231 16         61 $self->element(
232             {
233             'Name' => 'HMMER_query-def',
234             'Data' => $buffer
235             }
236             );
237             }
238             elsif ($buffer =~ s/^Accession:\s+//o) {
239 8         16 $buffer =~ s/\s+$//;
240 8         29 $self->element(
241             {
242             'Name' => 'HMMER_query-acc',
243             'Data' => $buffer
244             }
245             );
246             }
247             elsif ($buffer =~ s/^Description:\s+//o) {
248 8         22 $buffer =~ s/\s+$//;
249 8         26 $self->element(
250             {
251             'Name' => 'HMMER_querydesc',
252             'Data' => $buffer
253             }
254             );
255             }
256             elsif ( defined $self->{'_reporttype'}
257             && ( $self->{'_reporttype'} eq 'HMMSEARCH'
258             || $self->{'_reporttype'} eq 'HMMPFAM' )
259             ) {
260             # PROCESS RESULTS HERE
261 365 100 100     1159 if ($buffer =~ /^Scores for (?:complete sequences|sequence family)/o) {
    100          
    100          
    100          
262 16         36 while ( defined( $buffer = $self->_readline ) ) {
263 3173 100       5859 last if ($buffer =~ /^\s+$/);
264 3157 100 100     9200 next if ( $buffer =~ /^Model\s+Description/o
      100        
265             || $buffer =~ /^Sequence\s+Description/o
266             || $buffer =~ /^\-\-\-/o );
267              
268 3125         3093 chomp $buffer;
269 3125         13679 my @line = split( /\s+/, $buffer );
270 3125         5257 my ( $name, $domaintotal, $evalue, $score ) =
271             ( shift @line, pop @line, pop @line, pop @line );
272 3125         5961 my $desc = join( ' ', @line );
273 3125         6260 push @hitinfo, [ $name, $desc, $score, $evalue, $domaintotal ];
274 3125         8412 $hitinfo{$name} = $#hitinfo;
275             }
276             }
277             elsif ($buffer =~ /^Parsed for domains:/o) {
278 16         27 @hspinfo = ();
279              
280 16         42 while ( defined( $buffer = $self->_readline ) ) {
281 5035 100       9360 last if ($buffer =~ /^\s+$/);
282 5019 50       6090 if ($buffer =~ m!^//!) {
283 0         0 $self->_pushback($buffer);
284 0         0 last;
285             }
286 5019 100 100     10791 next if ( $buffer =~ /^(?:Model|Sequence)\s+Domain/ || $buffer =~ /^\-\-\-/ );
287              
288 4987         5106 chomp $buffer;
289 4987 50       31307 if (
290             my ( $name, $domainct, $domaintotal,
291             $seq_start, $seq_end, $seq_cov,
292             $hmm_start, $hmm_end, $hmm_cov,
293             $score, $evalue )
294             = ( $buffer =~
295             m!^(\S+)\s+ # domain name
296             (\d+)/(\d+)\s+ # domain num out of num
297             (\d+)\s+(\d+)\s+ # seq start, end
298             (\S+)\s+ # seq coverage
299             (\d+)\s+(\d+)\s+ # hmm start, end
300             (\S+)\s+ # hmm coverage
301             (\S+)\s+ # score
302             (\S+) # evalue
303             \s*$!ox
304             )
305             ) {
306 4987         8261 my $hindex = $hitinfo{$name};
307 4987 50       6267 if ( !defined $hindex ) {
308 0         0 push @hitinfo,
309             [ $name, '', $score, $evalue, $domaintotal ];
310 0         0 $hitinfo{$name} = $#hitinfo;
311 0         0 $hindex = $#hitinfo;
312             }
313              
314 4987         5264 my $info = $hitinfo[$hindex];
315 4987 50       5342 if ( !defined $info ) {
316 0 0       0 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    0          
317 0         0 $self->warn(
318             "Incomplete Sequence information, can't find $name hitinfo says $hitinfo{$name}"
319             );
320             }
321             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
322 0         0 $self->warn(
323             "Incomplete Domain information, can't find $name hitinfo says $hitinfo{$name}"
324             );
325             }
326 0         0 next;
327             }
328              
329             # Try to get HMM and Sequence lengths from the alignment information
330 4987 100       6666 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
331             # For Hmmsearch, if seq coverage ends in ']' it means that the alignment
332             # runs until the end. In that case add the END coordinate to @hitinfo
333             # to use it as Hit Length
334 4864 100 66     7254 if ( $seq_cov =~ m/\]$/
335 48         114 and scalar @{ $hitinfo[$hindex] } == 5
336             ) {
337 48         41 push @{ $hitinfo[$hindex] }, $seq_end ;
  48         93  
338             }
339             # For Hmmsearch, if hmm coverage ends in ']', it means that the alignment
340             # runs until the end. In that case use the END coordinate as Query Length
341 4864 100 66     14110 if ( $hmm_cov =~ m/\]$/
342             and not exists $self->{_values}->{'RESULT-query_length'}
343             ) {
344 5         33 $self->element(
345             { 'Name' => 'HMMER_query-len',
346             'Data' => $hmm_end
347             }
348             );
349             }
350             }
351             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
352             # For Hmmpfam, if hmm coverage ends in ']' it means that the alignment
353             # runs until the end. In that case add the END coordinate to @hitinfo
354             # to use it as Hit Length
355 123 100 100     240 if ( $hmm_cov =~ m/\]$/
356 73         156 and scalar @{ $hitinfo[$hindex] } == 5
357             ) {
358 69         64 push @{ $hitinfo[$hindex] }, $hmm_end ;
  69         113  
359             }
360             # For Hmmpfam, if seq coverage ends in ']', it means that the alignment
361             # runs until the end. In that case use the END coordinate as Query Length
362 123 100 66     213 if ( $seq_cov =~ m/\]$/
363             and not exists $self->{_values}->{'RESULT-query_length'}
364             ) {
365 2         11 $self->element(
366             { 'Name' => 'HMMER_query-len',
367             'Data' => $seq_end
368             }
369             );
370             }
371             }
372              
373 4987         10699 my @vals = ($seq_start, $seq_end,
374             $hmm_start, $hmm_end,
375             $score, $evalue);
376 4987         19369 push @hspinfo, [ $name, @vals ];
377             }
378             }
379             }
380             elsif ($buffer =~ /^Alignments of top/o) {
381 12         22 my ( $prelength, $count, $width );
382 12         15 $count = 0;
383 12         14 my %domaincounter;
384 12         18 my $second_tier = 0;
385 12         14 my $csline = '';
386              
387 12         27 while ( defined( $buffer = $self->_readline ) ) {
388 1204 50       1870 next if ( $buffer =~ /^Align/o );
389              
390 1204 100 100     3701 if ( $buffer =~ m/^Histogram/o
      66        
391             || $buffer =~ m!^//!o
392             || $buffer =~ m/^Query(?:\s+(?:sequence|HMM))?(?:\s+\d+)?:/o
393             ) {
394 12 50       18 if ( $self->in_element('hsp') ) {
395 12         31 $self->end_element( { 'Name' => 'Hsp' } );
396             }
397 12 50       37 if ( $self->within_element('hit') ) {
398 12         32 $self->end_element( { 'Name' => 'Hit' } );
399             }
400 12         50 $self->_pushback($buffer);
401 12         55 last;
402             }
403              
404 1192         1242 chomp $buffer;
405 1192 100       2724 if (
406             my ( $name, $domainct, $domaintotal,
407             $from, $to )
408             = ( $buffer =~
409             m/^\s*(.+):
410             \s+ domain \s+ (\d+) \s+ of \s+ (\d+) ,
411             \s+ from \s+ (\d+) \s+ to \s+ (\d+)/x
412             )
413             ) {
414 127         284 $domaincounter{$name}++;
415 127 100       274 if ( $self->within_element('hit') ) {
416 115 50       158 if ( $self->within_element('hsp') ) {
417 115         272 $self->end_element( { 'Name' => 'Hsp' } );
418             }
419 115         428 $self->end_element( { 'Name' => 'Hit' } );
420             }
421              
422 127         197 my $info = [ @{ $hitinfo[ $hitinfo{$name} ] } ];
  127         548  
423 127 50 33     445 if ( !defined $info
424             || $info->[0] ne $name
425             ) {
426 0         0 $self->warn(
427             "Somehow the Model table order does not match the order in the domains (got "
428             . $info->[0]
429             . ", expected $name). We're back loading this from the alignment information instead"
430             );
431 0         0 $info = [
432             $name, '',
433             $buffer =~ /score \s+ ([^,\s]+), \s+E\s+=\s+ (\S+)/ox,
434             $domaintotal
435             ];
436 0         0 push @hitinfo, $info;
437 0         0 $hitinfo{$name} = $#hitinfo;
438             }
439              
440 127         379 $self->start_element( { 'Name' => 'Hit' } );
441             $self->element(
442             {
443             'Name' => 'Hit_id',
444 127         192 'Data' => shift @{$info}
  127         385  
445             }
446             );
447             $self->element(
448             {
449             'Name' => 'Hit_desc',
450 127         227 'Data' => shift @{$info}
  127         334  
451             }
452             );
453             $self->element(
454             {
455             'Name' => 'Hit_score',
456 127         213 'Data' => shift @{$info}
  127         321  
457             }
458             );
459             $self->element(
460             {
461             'Name' => 'Hit_signif',
462 127         217 'Data' => shift @{$info}
  127         305  
463             }
464             );
465 127         190 my $dom_total = shift @{$info};
  127         199  
466 127 100       155 if (my $hit_end = shift @{$info}) {
  127         240  
467 75         165 $self->element(
468             {
469             'Name' => 'Hit_len',
470             'Data' => $hit_end
471             }
472             );
473             }
474              
475 127         385 $self->start_element( { 'Name' => 'Hsp' } );
476 127         221 my $HSPinfo = shift @hspinfo;
477 127         204 my $id = shift @$HSPinfo;
478              
479 127 50       209 if ( $id ne $name ) {
480 0         0 $self->throw(
481             "Somehow the domain list details do not match "
482             . "the table (got $id, expected $name)"
483             );
484             }
485              
486 127 100       335 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
487 4         11 $self->element(
488             {
489             'Name' => 'Hsp_hit-from',
490             'Data' => shift @$HSPinfo
491             }
492             );
493 4         13 $self->element(
494             {
495             'Name' => 'Hsp_hit-to',
496             'Data' => shift @$HSPinfo
497             }
498             );
499 4         12 $self->element(
500             {
501             'Name' => 'Hsp_query-from',
502             'Data' => shift @$HSPinfo
503             }
504             );
505 4         12 $self->element(
506             {
507             'Name' => 'Hsp_query-to',
508             'Data' => shift @$HSPinfo
509             }
510             );
511             }
512             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
513 123         327 $self->element(
514             {
515             'Name' => 'Hsp_query-from',
516             'Data' => shift @$HSPinfo
517             }
518             );
519 123         391 $self->element(
520             {
521             'Name' => 'Hsp_query-to',
522             'Data' => shift @$HSPinfo
523             }
524             );
525 123         415 $self->element(
526             {
527             'Name' => 'Hsp_hit-from',
528             'Data' => shift @$HSPinfo
529             }
530             );
531 123         339 $self->element(
532             {
533             'Name' => 'Hsp_hit-to',
534             'Data' => shift @$HSPinfo
535             }
536             );
537             }
538             $self->element(
539             {
540 127         392 'Name' => 'Hsp_score',
541             'Data' => shift @$HSPinfo
542             }
543             );
544 127         366 $self->element(
545             {
546             'Name' => 'Hsp_evalue',
547             'Data' => shift @$HSPinfo
548             }
549             );
550              
551 127 100       376 if ( $domaincounter{$name} == $domaintotal ) {
552 121         623 $hitinfo[ $hitinfo{$name} ] = undef;
553             }
554             }
555             else {
556              
557             # Might want to change this so that it
558             # accumulates all the of the alignment lines into
559             # three array slots and then tests for the
560             # end of the line
561 1065 100 66     7734 if ($buffer =~ m/^\s+(?:CS|RF)\s+/o && $count == 0) {
    100 100        
    100 100        
    100 100        
    100 100        
    100          
    50          
562             # Buffer the CS line now and process it later at
563             # midline point, where $prelength and width will be known
564 104         148 $csline = $buffer;
565 104         213 next;
566             }
567             elsif ($buffer =~ /^(\s+ \*->) (\S+)/ox) {
568             # start of domain
569 127         299 $prelength = CORE::length($1);
570 127         141 $width = 0;
571              
572             # deal with fact that start and stop is on same line
573 127         218 my $data = $2;
574 127 100       484 if ($data =~ s/<-?\*?\s*$//)
575             {
576 102         141 $width = CORE::length($data);
577             }
578              
579 127 100       352 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
580 4         14 $self->element(
581             {
582             'Name' => 'Hsp_qseq',
583             'Data' => $data
584             }
585             );
586             }
587             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
588 123         392 $self->element(
589             {
590             'Name' => 'Hsp_hseq',
591             'Data' => $data
592             }
593             );
594             }
595 127         237 $count = 0;
596 127         145 $second_tier = 0;
597             }
598             elsif ($buffer =~ /^(\s+) (\S+) <-?\*? \s*$/ox) {
599             # end of domain
600 25 100       55 $prelength -= 3 unless ( $second_tier++ );
601 25 100       65 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
602 4         14 $self->element(
603             {
604             'Name' => 'Hsp_qseq',
605             'Data' => $2
606             }
607             );
608             }
609             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
610 21         92 $self->element(
611             {
612             'Name' => 'Hsp_hseq',
613             'Data' => $2
614             }
615             );
616             }
617 25         61 $width = CORE::length($2);
618 25         26 $count = 0;
619             }
620             elsif ( ( $count != 1 && $buffer =~ /^\s+$/o )
621             || CORE::length($buffer) == 0
622             || $buffer =~ /^\s+\-?\*\s*$/
623             || $buffer =~ /^\s+\S+\s+\-\s+\-\s*$/ )
624             {
625 247         524 next;
626             }
627             elsif ( $count == 0 ) {
628 86 100       164 $prelength -= 3 unless ( $second_tier++ );
629 86 50       137 unless ( defined $prelength ) {
630              
631             # $self->warn("prelength not set");
632 0         0 next;
633             }
634 86 100       216 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
635 20         64 $self->element(
636             {
637             'Name' => 'Hsp_qseq',
638             'Data' => substr( $buffer, $prelength )
639             }
640             );
641             }
642             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
643 66         252 $self->element(
644             {
645             'Name' => 'Hsp_hseq',
646             'Data' => substr( $buffer, $prelength )
647             }
648             );
649             }
650             }
651             elsif ( $count == 1 ) {
652 238 50       392 if ( !defined $prelength ) {
653 0         0 $self->warn("prelength not set");
654             }
655 238 100       281 if ($width) {
656 127         519 $self->element(
657             {
658             'Name' => 'Hsp_midline',
659             'Data' => substr( $buffer, $prelength, $width )
660             }
661             );
662 127 100       301 if ($csline ne '') {
663 54         209 $self->element(
664             {
665             'Name' => 'Hsp_csline',
666             'Data' => substr( $csline, $prelength, $width )
667              
668             }
669             );
670 54         105 $csline = '';
671             }
672             }
673             else {
674 111         401 $self->element(
675             {
676             'Name' => 'Hsp_midline',
677             'Data' => substr( $buffer, $prelength )
678             }
679             );
680 111 100       224 if ($csline ne '') {
681 50         146 $self->element(
682             {
683             'Name' => 'Hsp_csline',
684             'Data' => substr( $csline, $prelength )
685             }
686             );
687 50         80 $csline = '';
688             }
689             }
690             }
691             elsif ( $count == 2 ) {
692 238 50       732 if ( $buffer =~ /^\s+(\S+)\s+(\d+|\-)\s+(\S*)\s+(\d+|\-)/o ) {
693 238 100       545 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
694 28         93 $self->element(
695             {
696             'Name' => 'Hsp_hseq',
697             'Data' => $3
698             }
699             );
700             }
701             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
702 210         760 $self->element(
703             {
704             'Name' => 'Hsp_qseq',
705             'Data' => $3
706             }
707             );
708             }
709             }
710             else {
711 0         0 $self->warn("unrecognized line ($count): $buffer\n");
712             }
713             }
714 714 100       2045 $count = 0 if $count++ >= 2;
715             }
716             }
717             }
718             elsif ( $buffer =~ /^Histogram/o || $buffer =~ m!^//!o ) {
719 16         23 my %domaincounter;
720              
721 16         47 while ( my $HSPinfo = shift @hspinfo ) {
722 4860         7050 my $id = shift @$HSPinfo;
723 4860         9133 $domaincounter{$id}++;
724              
725 4860         4966 my $info = [ @{ $hitinfo[ $hitinfo{$id} ] } ];
  4860         18421  
726 4860 50       7920 next unless defined $info;
727              
728 4860         12259 $self->start_element( { 'Name' => 'Hit' } );
729             $self->element(
730             {
731             'Name' => 'Hit_id',
732 4860         7517 'Data' => shift @{$info}
  4860         12844  
733             }
734             );
735             $self->element(
736             {
737             'Name' => 'Hit_desc',
738 4860         7935 'Data' => shift @{$info}
  4860         11998  
739             }
740             );
741             $self->element(
742             {
743             'Name' => 'Hit_score',
744 4860         7913 'Data' => shift @{$info}
  4860         11309  
745             }
746             );
747             $self->element(
748             {
749             'Name' => 'Hit_signif',
750 4860         7627 'Data' => shift @{$info}
  4860         10833  
751             }
752             );
753 4860         6592 my $domaintotal = shift @{$info};
  4860         5987  
754 4860 100       5094 if (my $hit_end = shift @{$info}) {
  4860         7892  
755 56         208 $self->element(
756             {
757             'Name' => 'Hit_len',
758             'Data' => $hit_end
759             }
760             );
761             }
762              
763             # Histogram is exclusive of Hmmsearch, not found in Hmmpfam,
764             # so just use Hmmsearch start/end order (first hit, then query)
765 4860         11633 $self->start_element( { 'Name' => 'Hsp' } );
766 4860         13902 $self->element(
767             {
768             'Name' => 'Hsp_hit-from',
769             'Data' => shift @$HSPinfo
770             }
771             );
772 4860         15266 $self->element(
773             {
774             'Name' => 'Hsp_hit-to',
775             'Data' => shift @$HSPinfo
776             }
777             );
778 4860         14221 $self->element(
779             {
780             'Name' => 'Hsp_query-from',
781             'Data' => shift @$HSPinfo
782             }
783             );
784 4860         13349 $self->element(
785             {
786             'Name' => 'Hsp_query-to',
787             'Data' => shift @$HSPinfo
788             }
789             );
790 4860         13719 $self->element(
791             {
792             'Name' => 'Hsp_score',
793             'Data' => shift @$HSPinfo
794             }
795             );
796 4860         13864 $self->element(
797             {
798             'Name' => 'Hsp_evalue',
799             'Data' => shift @$HSPinfo
800             }
801             );
802 4860         11888 $self->end_element( { 'Name' => 'Hsp' } );
803 4860         15397 $self->end_element( { 'Name' => 'Hit' } );
804              
805 4860 100       17275 if ( $domaincounter{$id} == $domaintotal ) {
806 3004         12162 $hitinfo[ $hitinfo{$id} ] = undef;
807             }
808             }
809 16         235 @hitinfo = ();
810 16         677 %hitinfo = ();
811 16         931 last;
812             }
813             # uncomment to see missed lines with verbose on
814             #else {
815             # $self->debug($buffer);
816             #}
817             }
818 444         794 $last = $buffer;
819             }
820 19 100       82 $self->end_element( { 'Name' => 'HMMER_Output' } ) unless !$seentop;
821 19         58 return $self->end_document();
822             }
823              
824             =head2 start_element
825              
826             Title : start_element
827             Usage : $eventgenerator->start_element
828             Function: Handles a start element event
829             Returns : none
830             Args : hashref with at least 2 keys 'Data' and 'Name'
831              
832              
833             =cut
834              
835             sub start_element {
836 90966     90966 1 90192 my ( $self, $data ) = @_;
837              
838             # we currently don't care about attributes
839 90966         87380 my $nm = $data->{'Name'};
840 90966         85095 my $type = $MODEMAP{$nm};
841 90966 100       105460 if ($type) {
842 9990 50       17881 if ( $self->_eventHandler->will_handle($type) ) {
843 9990         23293 my $func = sprintf( "start_%s", lc $type );
844 9990         15712 $self->_eventHandler->$func( $data->{'Attributes'} );
845             }
846 9990         14895 unshift @{ $self->{'_elements'} }, $type;
  9990         19148  
847             }
848 90966 100 100     143244 if ( defined $type
849             && $type eq 'result' )
850             {
851 16         33 $self->{'_values'} = {};
852 16         27 $self->{'_result'} = undef;
853             }
854             }
855              
856             =head2 end_element
857              
858             Title : start_element
859             Usage : $eventgenerator->end_element
860             Function: Handles an end element event
861             Returns : none
862             Args : hashref with at least 2 keys 'Data' and 'Name'
863              
864              
865             =cut
866              
867             sub end_element {
868 90966     90966 1 98986 my ( $self, $data ) = @_;
869 90966         94228 my $nm = $data->{'Name'};
870 90966         84433 my $type = $MODEMAP{$nm};
871 90966         71598 my $rc;
872              
873 90966 100       108399 if ( $nm eq 'HMMER_program' ) {
874 16 100       70 if ( $self->{'_last_data'} =~ /(HMM\S+)/i ) {
875 15         44 $self->{'_reporttype'} = uc $1;
876             }
877             }
878              
879             # Hsp are sort of weird, in that they end when another
880             # object begins so have to detect this in end_element for now
881 90966 100       103644 if ( $nm eq 'Hsp' ) {
882 4987         6939 foreach my $line (qw(Hsp_csline Hsp_qseq Hsp_midline Hsp_hseq)) {
883 19948         21736 my $data = $self->{'_last_hspdata'}->{$line};
884 19948 100 100     26314 if ($data && $line eq 'Hsp_hseq') {
885             # replace hmm '.' gap symbol by '-'
886 127         284 $data =~ s/\./-/g;
887             }
888             $self->element(
889             {
890 19948         40461 'Name' => $line,
891             'Data' => $data
892             }
893             );
894             # Since HMMER doesn't print some data explicitly,
895             # calculate it from the homology line (midline)
896 19948 100       36047 if ($line eq 'Hsp_midline') {
897 4987 100       6627 if ($data) {
898 127         133 my $length = length $data;
899 127         180 my $identical = ($data =~ tr/a-zA-Z//);
900 127         174 my $positive = ($data =~ tr/+//) + $identical;
901 127         383 $self->element(
902             {
903             'Name' => 'Hsp_align-len',
904             'Data' => $length
905             }
906             );
907 127         412 $self->element(
908             { 'Name' => 'Hsp_identity',
909             'Data' => $identical
910             }
911             );
912 127         343 $self->element(
913             { 'Name' => 'Hsp_positive',
914             'Data' => $positive
915             }
916             );
917             }
918             else {
919 4860         11196 $self->element(
920             { 'Name' => 'Hsp_identity',
921             'Data' => 0
922             }
923             );
924 4860         11976 $self->element(
925             { 'Name' => 'Hsp_positive',
926             'Data' => 0
927             }
928             );
929             }
930             }
931             }
932 4987         7335 $self->{'_last_hspdata'} = {};
933             }
934 90966 100       122135 if ($type) {
    50          
935 9990 50       20503 if ( $self->_eventHandler->will_handle($type) ) {
936 9990         24973 my $func = sprintf( "end_%s", lc $type );
937             $rc = $self->_eventHandler->$func( $self->{'_reporttype'},
938 9990         13757 $self->{'_values'} );
939             }
940 9990         11841 my $lastelem = shift @{ $self->{'_elements'} };
  9990         16933  
941              
942             # Flush corresponding values from the {_values} buffer
943 9990         12428 my $name = uc $type;
944 9990         10155 foreach my $key (keys %{ $self->{_values} }) {
  9990         31113  
945 194851 100       527504 delete $self->{_values}->{$key} if ($key =~ m/^$name-/);
946             }
947             }
948             elsif ( $MAPPING{$nm} ) {
949 80976 50       98129 if ( ref( $MAPPING{$nm} ) =~ /hash/i ) {
950 0         0 my $key = ( keys %{ $MAPPING{$nm} } )[0];
  0         0  
951             $self->{'_values'}->{$key}->{ $MAPPING{$nm}->{$key} } =
952 0         0 $self->{'_last_data'};
953             }
954             else {
955 80976         139195 $self->{'_values'}->{ $MAPPING{$nm} } = $self->{'_last_data'};
956             }
957             }
958             else {
959 0         0 $self->debug("unknown nm $nm, ignoring\n");
960             }
961 90966         99068 $self->{'_last_data'} = ''; # remove read data if we are at
962             # end of an element
963 90966 100 100     130575 $self->{'_result'} = $rc if ( defined $type && $type eq 'result' );
964 90966         113693 return $rc;
965             }
966              
967             =head2 element
968              
969             Title : element
970             Usage : $eventhandler->element({'Name' => $name, 'Data' => $str});
971             Function: Convience method that calls start_element, characters, end_element
972             Returns : none
973             Args : Hash ref with the keys 'Name' and 'Data'
974              
975              
976             =cut
977              
978             sub element {
979 80976     80976 1 89849 my ( $self, $data ) = @_;
980 80976         110212 $self->start_element($data);
981 80976         112821 $self->characters($data);
982 80976         112828 $self->end_element($data);
983             }
984              
985             =head2 characters
986              
987             Title : characters
988             Usage : $eventgenerator->characters($str)
989             Function: Send a character events
990             Returns : none
991             Args : string
992              
993              
994             =cut
995              
996             sub characters {
997 80976     80976 1 82990 my ( $self, $data ) = @_;
998              
999 80976 100 100     86829 if ( $self->in_element('hsp')
      100        
1000             && $data->{'Name'} =~ /Hsp\_(?:qseq|hseq|csline|midline)/o
1001             && defined $data->{'Data'} )
1002             {
1003 1253         3123 $self->{'_last_hspdata'}->{ $data->{'Name'} } .= $data->{'Data'};
1004             }
1005 80976 100 100     228011 return unless ( defined $data->{'Data'} && $data->{'Data'} !~ /^\s+$/o );
1006              
1007 61460         81871 $self->{'_last_data'} = $data->{'Data'};
1008             }
1009              
1010             =head2 within_element
1011              
1012             Title : within_element
1013             Usage : if( $eventgenerator->within_element($element) ) {}
1014             Function: Test if we are within a particular element
1015             This is different than 'in' because within can be tested
1016             for a whole block.
1017             Returns : boolean
1018             Args : string element name
1019              
1020              
1021             =cut
1022              
1023             sub within_element {
1024 254     254 1 370 my ( $self, $name ) = @_;
1025             return 0
1026             if ( !defined $name
1027             || !defined $self->{'_elements'}
1028 254 50 33     685 || scalar @{ $self->{'_elements'} } == 0 );
  254   33     573  
1029 254         283 foreach my $element ( @{ $self->{'_elements'} } ) {
  254         389  
1030 369 100       814 return 1 if ( $element eq $name );
1031             }
1032 12         28 return 0;
1033             }
1034              
1035             =head2 in_element
1036              
1037             Title : in_element
1038             Usage : if( $eventgenerator->in_element($element) ) {}
1039             Function: Test if we are in a particular element
1040             This is different than 'within' because 'in' only
1041             tests its immediete parent.
1042             Returns : boolean
1043             Args : string element name
1044              
1045              
1046             =cut
1047              
1048             sub in_element {
1049 80988     80988 1 84872 my ( $self, $name ) = @_;
1050 80988 50       109621 return 0 if !defined $self->{'_elements'}->[0];
1051 80988         295651 return ( $self->{'_elements'}->[0] eq $name );
1052             }
1053              
1054             =head2 start_document
1055              
1056             Title : start_document
1057             Usage : $eventgenerator->start_document
1058             Function: Handle a start document event
1059             Returns : none
1060             Args : none
1061              
1062              
1063             =cut
1064              
1065             sub start_document {
1066 19     19 1 30 my ($self) = @_;
1067 19         45 $self->{'_lasttype'} = '';
1068 19         44 $self->{'_values'} = {};
1069 19         37 $self->{'_result'} = undef;
1070 19         40 $self->{'_elements'} = [];
1071             }
1072              
1073             =head2 end_document
1074              
1075             Title : end_document
1076             Usage : $eventgenerator->end_document
1077             Function: Handles an end document event
1078             Returns : Bio::Search::Result::ResultI object
1079             Args : none
1080              
1081              
1082             =cut
1083              
1084             sub end_document {
1085 19     19 1 36 my ($self) = @_;
1086 19         147 return $self->{'_result'};
1087             }
1088              
1089             =head2 result_count
1090              
1091             Title : result_count
1092             Usage : my $count = $searchio->result_count
1093             Function: Returns the number of results we have processed
1094             Returns : integer
1095             Args : none
1096              
1097              
1098             =cut
1099              
1100             sub result_count {
1101 0     0 1   my $self = shift;
1102 0           return $self->{'_result_count'};
1103             }
1104              
1105             1;