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   4 use strict;
  1         2  
  1         30  
85              
86 1     1   3 use Bio::Factory::ObjectFactory;
  1         1  
  1         18  
87              
88 1         51 use vars qw(%MAPPING %MODEMAP
89 1     1   3 );
  1         1  
90              
91 1     1   3 use base qw(Bio::SearchIO::hmmer);
  1         1  
  1         139  
92              
93             BEGIN {
94              
95             # mapping of HMMER items to Bioperl hash keys
96 1     1   5 %MODEMAP = (
97             'HMMER_Output' => 'result',
98             'Hit' => 'hit',
99             'Hsp' => 'hsp'
100             );
101              
102 1         3020 %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 775 my ($self) = @_;
154 19         28 my $seentop = 0;
155 19         20 my $reporttype;
156 19         25 my ( $buffer, $last, @hitinfo, @hspinfo, %hspinfo, %hitinfo );
157 19         75 local $/ = "\n";
158              
159 19         57 my $verbose = $self->verbose; # cache for speed?
160 19         63 $self->start_document();
161              
162 19         58 while ( defined( $buffer = $self->_readline ) ) {
163 461         375 my $lineorig = $buffer;
164              
165 461         407 chomp $buffer;
166 461 100 100     2811 if ($buffer =~ /^HMMER\s+(\S+)\s+\((.+)\)/o) {
    100 33        
    100          
    100          
    100          
    100          
    100          
167 16         81 my ( $prog, $version ) = split( /\s+/, $buffer );
168 16 50       36 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         33 $self->{'_hmmidline'} = $buffer;
174 16         77 $self->start_element( { 'Name' => 'HMMER_Output' } );
175 16         36 $self->{'_result_count'}++;
176 16         19 $seentop = 1;
177 16 50       70 if ( defined $last ) {
178 16         69 ($reporttype) = split( /\s+/, $last );
179 16 100       52 $reporttype = uc($reporttype) if defined $reporttype;
180 16         81 $self->element(
181             {
182             'Name' => 'HMMER_program',
183             'Data' => $reporttype
184             }
185             );
186             }
187             $self->element(
188             {
189 16         62 'Name' => 'HMMER_version',
190             'Data' => $version
191             }
192             );
193             }
194             elsif ($buffer =~ s/^HMM file:\s+//o) {
195 16         43 $self->{'_hmmfileline'} = $lineorig;
196 16         73 $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         33 $self->{'_hmmseqline'} = $lineorig;
205 16 100       53 if ( $1 eq 'database' ) {
206 5         15 $self->element(
207             {
208             'Name' => 'HMMER_db',
209             'Data' => $buffer
210             }
211             );
212             }
213             $self->element(
214             {
215 16         61 'Name' => 'HMMER_seqfile',
216             'Data' => $buffer
217             }
218             );
219             }
220             elsif ($buffer =~ s/^Query(?:\s+(?:sequence|HMM))?(?:\s+\d+)?:\s+//o) {
221 17 100       37 if ( !$seentop ) {
222              
223             # we're in a multi-query report
224 1         6 $self->_pushback($lineorig);
225 1         8 $self->_pushback( $self->{'_hmmseqline'} );
226 1         3 $self->_pushback( $self->{'_hmmfileline'} );
227 1         4 $self->_pushback( $self->{'_hmmidline'} );
228 1         4 next;
229             }
230 16         39 $buffer =~ s/\s+$//;
231 16         62 $self->element(
232             {
233             'Name' => 'HMMER_query-def',
234             'Data' => $buffer
235             }
236             );
237             }
238             elsif ($buffer =~ s/^Accession:\s+//o) {
239 8         14 $buffer =~ s/\s+$//;
240 8         27 $self->element(
241             {
242             'Name' => 'HMMER_query-acc',
243             'Data' => $buffer
244             }
245             );
246             }
247             elsif ($buffer =~ s/^Description:\s+//o) {
248 8         24 $buffer =~ s/\s+$//;
249 8         32 $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     1472 if ($buffer =~ /^Scores for (?:complete sequences|sequence family)/o) {
    100          
    100          
    100          
262 16         37 while ( defined( $buffer = $self->_readline ) ) {
263 3173 100       5700 last if ($buffer =~ /^\s+$/);
264 3157 100 100     11873 next if ( $buffer =~ /^Model\s+Description/o
      100        
265             || $buffer =~ /^Sequence\s+Description/o
266             || $buffer =~ /^\-\-\-/o );
267              
268 3125         2235 chomp $buffer;
269 3125         12950 my @line = split( /\s+/, $buffer );
270 3125         3422 my ( $name, $domaintotal, $evalue, $score ) =
271             ( shift @line, pop @line, pop @line, pop @line );
272 3125         4911 my $desc = join( ' ', @line );
273 3125         5326 push @hitinfo, [ $name, $desc, $score, $evalue, $domaintotal ];
274 3125         8582 $hitinfo{$name} = $#hitinfo;
275             }
276             }
277             elsif ($buffer =~ /^Parsed for domains:/o) {
278 16         23 @hspinfo = ();
279              
280 16         45 while ( defined( $buffer = $self->_readline ) ) {
281 5035 100       8310 last if ($buffer =~ /^\s+$/);
282 5019 50       5452 if ($buffer =~ m!^//!) {
283 0         0 $self->_pushback($buffer);
284 0         0 last;
285             }
286 5019 100 100     12978 next if ( $buffer =~ /^(?:Model|Sequence)\s+Domain/ || $buffer =~ /^\-\-\-/ );
287              
288 4987         3635 chomp $buffer;
289 4987 50       28989 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         5883 my $hindex = $hitinfo{$name};
307 4987 50       5563 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         4086 my $info = $hitinfo[$hindex];
315 4987 50       5162 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       5675 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     6933 if ( $seq_cov =~ m/\]$/
335 48         125 and scalar @{ $hitinfo[$hindex] } == 5
336             ) {
337 48         34 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     16002 if ( $hmm_cov =~ m/\]$/
342             and not exists $self->{_values}->{'RESULT-query_length'}
343             ) {
344 5         34 $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     230 if ( $hmm_cov =~ m/\]$/
356 73         217 and scalar @{ $hitinfo[$hindex] } == 5
357             ) {
358 69         48 push @{ $hitinfo[$hindex] }, $hmm_end ;
  69         120  
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     221 if ( $seq_cov =~ m/\]$/
363             and not exists $self->{_values}->{'RESULT-query_length'}
364             ) {
365 2         30 $self->element(
366             { 'Name' => 'HMMER_query-len',
367             'Data' => $seq_end
368             }
369             );
370             }
371             }
372              
373 4987         8141 my @vals = ($seq_start, $seq_end,
374             $hmm_start, $hmm_end,
375             $score, $evalue);
376 4987         18676 push @hspinfo, [ $name, @vals ];
377             }
378             }
379             }
380             elsif ($buffer =~ /^Alignments of top/o) {
381 12         17 my ( $prelength, $count, $width );
382 12         13 $count = 0;
383 12         22 my %domaincounter;
384 12         11 my $second_tier = 0;
385 12         12 my $csline = '';
386              
387 12         30 while ( defined( $buffer = $self->_readline ) ) {
388 1204 50       1789 next if ( $buffer =~ /^Align/o );
389              
390 1204 100 100     4999 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         43 $self->end_element( { 'Name' => 'Hsp' } );
396             }
397 12 50       38 if ( $self->within_element('hit') ) {
398 12         36 $self->end_element( { 'Name' => 'Hit' } );
399             }
400 12         51 $self->_pushback($buffer);
401 12         49 last;
402             }
403              
404 1192         1189 chomp $buffer;
405 1192 100       2630 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         256 $domaincounter{$name}++;
415 127 100       262 if ( $self->within_element('hit') ) {
416 115 50       157 if ( $self->within_element('hsp') ) {
417 115         348 $self->end_element( { 'Name' => 'Hsp' } );
418             }
419 115         431 $self->end_element( { 'Name' => 'Hit' } );
420             }
421              
422 127         196 my $info = [ @{ $hitinfo[ $hitinfo{$name} ] } ];
  127         603  
423 127 50 33     520 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         383 $self->start_element( { 'Name' => 'Hit' } );
441             $self->element(
442             {
443             'Name' => 'Hit_id',
444 127         201 'Data' => shift @{$info}
  127         351  
445             }
446             );
447             $self->element(
448             {
449             'Name' => 'Hit_desc',
450 127         198 'Data' => shift @{$info}
  127         293  
451             }
452             );
453             $self->element(
454             {
455             'Name' => 'Hit_score',
456 127         172 'Data' => shift @{$info}
  127         283  
457             }
458             );
459             $self->element(
460             {
461             'Name' => 'Hit_signif',
462 127         248 'Data' => shift @{$info}
  127         284  
463             }
464             );
465 127         157 my $dom_total = shift @{$info};
  127         190  
466 127 100       159 if (my $hit_end = shift @{$info}) {
  127         275  
467 75         238 $self->element(
468             {
469             'Name' => 'Hit_len',
470             'Data' => $hit_end
471             }
472             );
473             }
474              
475 127         327 $self->start_element( { 'Name' => 'Hsp' } );
476 127         196 my $HSPinfo = shift @hspinfo;
477 127         180 my $id = shift @$HSPinfo;
478              
479 127 50       242 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       400 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
487 4         10 $self->element(
488             {
489             'Name' => 'Hsp_hit-from',
490             'Data' => shift @$HSPinfo
491             }
492             );
493 4         11 $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         11 $self->element(
506             {
507             'Name' => 'Hsp_query-to',
508             'Data' => shift @$HSPinfo
509             }
510             );
511             }
512             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
513 123         379 $self->element(
514             {
515             'Name' => 'Hsp_query-from',
516             'Data' => shift @$HSPinfo
517             }
518             );
519 123         419 $self->element(
520             {
521             'Name' => 'Hsp_query-to',
522             'Data' => shift @$HSPinfo
523             }
524             );
525 123         358 $self->element(
526             {
527             'Name' => 'Hsp_hit-from',
528             'Data' => shift @$HSPinfo
529             }
530             );
531 123         395 $self->element(
532             {
533             'Name' => 'Hsp_hit-to',
534             'Data' => shift @$HSPinfo
535             }
536             );
537             }
538             $self->element(
539             {
540 127         421 'Name' => 'Hsp_score',
541             'Data' => shift @$HSPinfo
542             }
543             );
544 127         382 $self->element(
545             {
546             'Name' => 'Hsp_evalue',
547             'Data' => shift @$HSPinfo
548             }
549             );
550              
551 127 100       459 if ( $domaincounter{$name} == $domaintotal ) {
552 121         691 $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     10112 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         155 $csline = $buffer;
565 104         295 next;
566             }
567             elsif ($buffer =~ /^(\s+ \*->) (\S+)/ox) {
568             # start of domain
569 127         323 $prelength = CORE::length($1);
570 127         134 $width = 0;
571              
572             # deal with fact that start and stop is on same line
573 127         187 my $data = $2;
574 127 100       572 if ($data =~ s/<-?\*?\s*$//)
575             {
576 102         127 $width = CORE::length($data);
577             }
578              
579 127 100       408 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
580 4         10 $self->element(
581             {
582             'Name' => 'Hsp_qseq',
583             'Data' => $data
584             }
585             );
586             }
587             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
588 123         385 $self->element(
589             {
590             'Name' => 'Hsp_hseq',
591             'Data' => $data
592             }
593             );
594             }
595 127         181 $count = 0;
596 127         135 $second_tier = 0;
597             }
598             elsif ($buffer =~ /^(\s+) (\S+) <-?\*? \s*$/ox) {
599             # end of domain
600 25 100       57 $prelength -= 3 unless ( $second_tier++ );
601 25 100       80 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
602 4         13 $self->element(
603             {
604             'Name' => 'Hsp_qseq',
605             'Data' => $2
606             }
607             );
608             }
609             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
610 21         83 $self->element(
611             {
612             'Name' => 'Hsp_hseq',
613             'Data' => $2
614             }
615             );
616             }
617 25         54 $width = CORE::length($2);
618 25         46 $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         557 next;
626             }
627             elsif ( $count == 0 ) {
628 86 100       182 $prelength -= 3 unless ( $second_tier++ );
629 86 50       149 unless ( defined $prelength ) {
630              
631             # $self->warn("prelength not set");
632 0         0 next;
633             }
634 86 100       210 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
635 20         55 $self->element(
636             {
637             'Name' => 'Hsp_qseq',
638             'Data' => substr( $buffer, $prelength )
639             }
640             );
641             }
642             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
643 66         269 $self->element(
644             {
645             'Name' => 'Hsp_hseq',
646             'Data' => substr( $buffer, $prelength )
647             }
648             );
649             }
650             }
651             elsif ( $count == 1 ) {
652 238 50       465 if ( !defined $prelength ) {
653 0         0 $self->warn("prelength not set");
654             }
655 238 100       285 if ($width) {
656 127         505 $self->element(
657             {
658             'Name' => 'Hsp_midline',
659             'Data' => substr( $buffer, $prelength, $width )
660             }
661             );
662 127 100       310 if ($csline ne '') {
663 54         210 $self->element(
664             {
665             'Name' => 'Hsp_csline',
666             'Data' => substr( $csline, $prelength, $width )
667              
668             }
669             );
670 54         90 $csline = '';
671             }
672             }
673             else {
674 111         423 $self->element(
675             {
676             'Name' => 'Hsp_midline',
677             'Data' => substr( $buffer, $prelength )
678             }
679             );
680 111 100       229 if ($csline ne '') {
681 50         160 $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       699 if ( $buffer =~ /^\s+(\S+)\s+(\d+|\-)\s+(\S*)\s+(\d+|\-)/o ) {
693 238 100       578 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
694 28         84 $self->element(
695             {
696             'Name' => 'Hsp_hseq',
697             'Data' => $3
698             }
699             );
700             }
701             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
702 210         862 $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       2343 $count = 0 if $count++ >= 2;
715             }
716             }
717             }
718             elsif ( $buffer =~ /^Histogram/o || $buffer =~ m!^//!o ) {
719 16         25 my %domaincounter;
720              
721 16         56 while ( my $HSPinfo = shift @hspinfo ) {
722 4860         5735 my $id = shift @$HSPinfo;
723 4860         8024 $domaincounter{$id}++;
724              
725 4860         3248 my $info = [ @{ $hitinfo[ $hitinfo{$id} ] } ];
  4860         17538  
726 4860 50       7403 next unless defined $info;
727              
728 4860         11002 $self->start_element( { 'Name' => 'Hit' } );
729             $self->element(
730             {
731             'Name' => 'Hit_id',
732 4860         5374 'Data' => shift @{$info}
  4860         10023  
733             }
734             );
735             $self->element(
736             {
737             'Name' => 'Hit_desc',
738 4860         6137 'Data' => shift @{$info}
  4860         8855  
739             }
740             );
741             $self->element(
742             {
743             'Name' => 'Hit_score',
744 4860         5651 'Data' => shift @{$info}
  4860         8881  
745             }
746             );
747             $self->element(
748             {
749             'Name' => 'Hit_signif',
750 4860         6830 'Data' => shift @{$info}
  4860         8820  
751             }
752             );
753 4860         5727 my $domaintotal = shift @{$info};
  4860         5018  
754 4860 100       3313 if (my $hit_end = shift @{$info}) {
  4860         6674  
755 56         141 $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         15473 $self->start_element( { 'Name' => 'Hsp' } );
766 4860         11251 $self->element(
767             {
768             'Name' => 'Hsp_hit-from',
769             'Data' => shift @$HSPinfo
770             }
771             );
772 4860         12019 $self->element(
773             {
774             'Name' => 'Hsp_hit-to',
775             'Data' => shift @$HSPinfo
776             }
777             );
778 4860         11465 $self->element(
779             {
780             'Name' => 'Hsp_query-from',
781             'Data' => shift @$HSPinfo
782             }
783             );
784 4860         11691 $self->element(
785             {
786             'Name' => 'Hsp_query-to',
787             'Data' => shift @$HSPinfo
788             }
789             );
790 4860         11501 $self->element(
791             {
792             'Name' => 'Hsp_score',
793             'Data' => shift @$HSPinfo
794             }
795             );
796 4860         11322 $self->element(
797             {
798             'Name' => 'Hsp_evalue',
799             'Data' => shift @$HSPinfo
800             }
801             );
802 4860         9770 $self->end_element( { 'Name' => 'Hsp' } );
803 4860         12360 $self->end_element( { 'Name' => 'Hit' } );
804              
805 4860 100       18263 if ( $domaincounter{$id} == $domaintotal ) {
806 3004         10934 $hitinfo[ $hitinfo{$id} ] = undef;
807             }
808             }
809 16         232 @hitinfo = ();
810 16         547 %hitinfo = ();
811 16         902 last;
812             }
813             # uncomment to see missed lines with verbose on
814             #else {
815             # $self->debug($buffer);
816             #}
817             }
818 444         769 $last = $buffer;
819             }
820 19 100       84 $self->end_element( { 'Name' => 'HMMER_Output' } ) unless !$seentop;
821 19         69 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 61223 my ( $self, $data ) = @_;
837              
838             # we currently don't care about attributes
839 90966         62408 my $nm = $data->{'Name'};
840 90966         66897 my $type = $MODEMAP{$nm};
841 90966 100       98412 if ($type) {
842 9990 50       16596 if ( $self->_eventHandler->will_handle($type) ) {
843 9990         19726 my $func = sprintf( "start_%s", lc $type );
844 9990         13556 $self->_eventHandler->$func( $data->{'Attributes'} );
845             }
846 9990         11212 unshift @{ $self->{'_elements'} }, $type;
  9990         15175  
847             }
848 90966 100 100     155914 if ( defined $type
849             && $type eq 'result' )
850             {
851 16         27 $self->{'_values'} = {};
852 16         69 $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 64661 my ( $self, $data ) = @_;
869 90966         69454 my $nm = $data->{'Name'};
870 90966         66480 my $type = $MODEMAP{$nm};
871 90966         52018 my $rc;
872              
873 90966 100       104381 if ( $nm eq 'HMMER_program' ) {
874 16 100       82 if ( $self->{'_last_data'} =~ /(HMM\S+)/i ) {
875 15         47 $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       95637 if ( $nm eq 'Hsp' ) {
882 4987         6261 foreach my $line (qw(Hsp_csline Hsp_qseq Hsp_midline Hsp_hseq)) {
883 19948         16021 my $data = $self->{'_last_hspdata'}->{$line};
884 19948 100 100     28375 if ($data && $line eq 'Hsp_hseq') {
885             # replace hmm '.' gap symbol by '-'
886 127         261 $data =~ s/\./-/g;
887             }
888             $self->element(
889             {
890 19948         33873 '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       35330 if ($line eq 'Hsp_midline') {
897 4987 100       5469 if ($data) {
898 127         148 my $length = length $data;
899 127         182 my $identical = ($data =~ tr/a-zA-Z//);
900 127         146 my $positive = ($data =~ tr/+//) + $identical;
901 127         298 $self->element(
902             {
903             'Name' => 'Hsp_align-len',
904             'Data' => $length
905             }
906             );
907 127         347 $self->element(
908             { 'Name' => 'Hsp_identity',
909             'Data' => $identical
910             }
911             );
912 127         337 $self->element(
913             { 'Name' => 'Hsp_positive',
914             'Data' => $positive
915             }
916             );
917             }
918             else {
919 4860         8691 $self->element(
920             { 'Name' => 'Hsp_identity',
921             'Data' => 0
922             }
923             );
924 4860         10033 $self->element(
925             { 'Name' => 'Hsp_positive',
926             'Data' => 0
927             }
928             );
929             }
930             }
931             }
932 4987         5880 $self->{'_last_hspdata'} = {};
933             }
934 90966 100       123790 if ($type) {
    50          
935 9990 50       17983 if ( $self->_eventHandler->will_handle($type) ) {
936 9990         20561 my $func = sprintf( "end_%s", lc $type );
937             $rc = $self->_eventHandler->$func( $self->{'_reporttype'},
938 9990         13577 $self->{'_values'} );
939             }
940 9990         8964 my $lastelem = shift @{ $self->{'_elements'} };
  9990         12249  
941              
942             # Flush corresponding values from the {_values} buffer
943 9990         9808 my $name = uc $type;
944 9990         6290 foreach my $key (keys %{ $self->{_values} }) {
  9990         27337  
945 194851 100       475599 delete $self->{_values}->{$key} if ($key =~ m/^$name-/);
946             }
947             }
948             elsif ( $MAPPING{$nm} ) {
949 80976 50       83731 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         132873 $self->{'_values'}->{ $MAPPING{$nm} } = $self->{'_last_data'};
956             }
957             }
958             else {
959 0         0 $self->debug("unknown nm $nm, ignoring\n");
960             }
961 90966         75598 $self->{'_last_data'} = ''; # remove read data if we are at
962             # end of an element
963 90966 100 100     137819 $self->{'_result'} = $rc if ( defined $type && $type eq 'result' );
964 90966         87576 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 58947 my ( $self, $data ) = @_;
980 80976         76406 $self->start_element($data);
981 80976         77027 $self->characters($data);
982 80976         85077 $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 54147 my ( $self, $data ) = @_;
998              
999 80976 100 100     72647 if ( $self->in_element('hsp')
      100        
1000             && $data->{'Name'} =~ /Hsp\_(?:qseq|hseq|csline|midline)/o
1001             && defined $data->{'Data'} )
1002             {
1003 1253         2865 $self->{'_last_hspdata'}->{ $data->{'Name'} } .= $data->{'Data'};
1004             }
1005 80976 100 100     239706 return unless ( defined $data->{'Data'} && $data->{'Data'} !~ /^\s+$/o );
1006              
1007 61460         69708 $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 314 my ( $self, $name ) = @_;
1025             return 0
1026             if ( !defined $name
1027             || !defined $self->{'_elements'}
1028 254 50 33     901 || scalar @{ $self->{'_elements'} } == 0 );
  254   33     642  
1029 254         228 foreach my $element ( @{ $self->{'_elements'} } ) {
  254         395  
1030 369 100       873 return 1 if ( $element eq $name );
1031             }
1032 12         34 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 57407 my ( $self, $name ) = @_;
1050 80988 50       103983 return 0 if !defined $self->{'_elements'}->[0];
1051 80988         321167 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 33 my ($self) = @_;
1067 19         41 $self->{'_lasttype'} = '';
1068 19         35 $self->{'_values'} = {};
1069 19         30 $self->{'_result'} = undef;
1070 19         35 $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 29 my ($self) = @_;
1086 19         135 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;