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         1  
  1         31  
85              
86 1     1   3 use Bio::Factory::ObjectFactory;
  1         1  
  1         18  
87              
88 1         42 use vars qw(%MAPPING %MODEMAP
89 1     1   2 );
  1         2  
90              
91 1     1   3 use base qw(Bio::SearchIO::hmmer);
  1         1  
  1         123  
92              
93             BEGIN {
94              
95             # mapping of HMMER items to Bioperl hash keys
96 1     1   3 %MODEMAP = (
97             'HMMER_Output' => 'result',
98             'Hit' => 'hit',
99             'Hsp' => 'hsp'
100             );
101              
102 1         2967 %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 732 my ($self) = @_;
154 19         23 my $seentop = 0;
155 19         15 my $reporttype;
156 19         20 my ( $buffer, $last, @hitinfo, @hspinfo, %hspinfo, %hitinfo );
157 19         62 local $/ = "\n";
158              
159 19         46 my $verbose = $self->verbose; # cache for speed?
160 19         52 $self->start_document();
161              
162 19         52 while ( defined( $buffer = $self->_readline ) ) {
163 461         389 my $lineorig = $buffer;
164              
165 461         377 chomp $buffer;
166 461 100 100     2638 if ($buffer =~ /^HMMER\s+(\S+)\s+\((.+)\)/o) {
    100 33        
    100          
    100          
    100          
    100          
    100          
167 16         68 my ( $prog, $version ) = split( /\s+/, $buffer );
168 16 50       40 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         22 $self->{'_hmmidline'} = $buffer;
174 16         56 $self->start_element( { 'Name' => 'HMMER_Output' } );
175 16         35 $self->{'_result_count'}++;
176 16         15 $seentop = 1;
177 16 50       30 if ( defined $last ) {
178 16         46 ($reporttype) = split( /\s+/, $last );
179 16 100       39 $reporttype = uc($reporttype) if defined $reporttype;
180 16         68 $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         31 $self->{'_hmmfileline'} = $lineorig;
196 16         43 $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         26 $self->{'_hmmseqline'} = $lineorig;
205 16 100       47 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         51 'Name' => 'HMMER_seqfile',
216             'Data' => $buffer
217             }
218             );
219             }
220             elsif ($buffer =~ s/^Query(?:\s+(?:sequence|HMM))?(?:\s+\d+)?:\s+//o) {
221 17 100       34 if ( !$seentop ) {
222              
223             # we're in a multi-query report
224 1         6 $self->_pushback($lineorig);
225 1         4 $self->_pushback( $self->{'_hmmseqline'} );
226 1         3 $self->_pushback( $self->{'_hmmfileline'} );
227 1         3 $self->_pushback( $self->{'_hmmidline'} );
228 1         3 next;
229             }
230 16         37 $buffer =~ s/\s+$//;
231 16         49 $self->element(
232             {
233             'Name' => 'HMMER_query-def',
234             'Data' => $buffer
235             }
236             );
237             }
238             elsif ($buffer =~ s/^Accession:\s+//o) {
239 8         10 $buffer =~ s/\s+$//;
240 8         26 $self->element(
241             {
242             'Name' => 'HMMER_query-acc',
243             'Data' => $buffer
244             }
245             );
246             }
247             elsif ($buffer =~ s/^Description:\s+//o) {
248 8         20 $buffer =~ s/\s+$//;
249 8         22 $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     1387 if ($buffer =~ /^Scores for (?:complete sequences|sequence family)/o) {
    100          
    100          
    100          
262 16         39 while ( defined( $buffer = $self->_readline ) ) {
263 3173 100       5290 last if ($buffer =~ /^\s+$/);
264 3157 100 100     11801 next if ( $buffer =~ /^Model\s+Description/o
      100        
265             || $buffer =~ /^Sequence\s+Description/o
266             || $buffer =~ /^\-\-\-/o );
267              
268 3125         2201 chomp $buffer;
269 3125         12661 my @line = split( /\s+/, $buffer );
270 3125         3336 my ( $name, $domaintotal, $evalue, $score ) =
271             ( shift @line, pop @line, pop @line, pop @line );
272 3125         5254 my $desc = join( ' ', @line );
273 3125         5295 push @hitinfo, [ $name, $desc, $score, $evalue, $domaintotal ];
274 3125         8563 $hitinfo{$name} = $#hitinfo;
275             }
276             }
277             elsif ($buffer =~ /^Parsed for domains:/o) {
278 16         23 @hspinfo = ();
279              
280 16         32 while ( defined( $buffer = $self->_readline ) ) {
281 5035 100       8652 last if ($buffer =~ /^\s+$/);
282 5019 50       5445 if ($buffer =~ m!^//!) {
283 0         0 $self->_pushback($buffer);
284 0         0 last;
285             }
286 5019 100 100     12212 next if ( $buffer =~ /^(?:Model|Sequence)\s+Domain/ || $buffer =~ /^\-\-\-/ );
287              
288 4987         3758 chomp $buffer;
289 4987 50       29113 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         5431 my $hindex = $hitinfo{$name};
307 4987 50       6087 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         3840 my $info = $hitinfo[$hindex];
315 4987 50       5175 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       5679 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     7051 if ( $seq_cov =~ m/\]$/
335 48         116 and scalar @{ $hitinfo[$hindex] } == 5
336             ) {
337 48         37 push @{ $hitinfo[$hindex] }, $seq_end ;
  48         82  
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     16273 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     244 if ( $hmm_cov =~ m/\]$/
356 73         166 and scalar @{ $hitinfo[$hindex] } == 5
357             ) {
358 69         45 push @{ $hitinfo[$hindex] }, $hmm_end ;
  69         108  
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     224 if ( $seq_cov =~ m/\]$/
363             and not exists $self->{_values}->{'RESULT-query_length'}
364             ) {
365 2         8 $self->element(
366             { 'Name' => 'HMMER_query-len',
367             'Data' => $seq_end
368             }
369             );
370             }
371             }
372              
373 4987         7976 my @vals = ($seq_start, $seq_end,
374             $hmm_start, $hmm_end,
375             $score, $evalue);
376 4987         18711 push @hspinfo, [ $name, @vals ];
377             }
378             }
379             }
380             elsif ($buffer =~ /^Alignments of top/o) {
381 12         12 my ( $prelength, $count, $width );
382 12         13 $count = 0;
383 12         11 my %domaincounter;
384 12         4 my $second_tier = 0;
385 12         11 my $csline = '';
386              
387 12         23 while ( defined( $buffer = $self->_readline ) ) {
388 1204 50       1593 next if ( $buffer =~ /^Align/o );
389              
390 1204 100 100     4340 if ( $buffer =~ m/^Histogram/o
      66        
391             || $buffer =~ m!^//!o
392             || $buffer =~ m/^Query(?:\s+(?:sequence|HMM))?(?:\s+\d+)?:/o
393             ) {
394 12 50       20 if ( $self->in_element('hsp') ) {
395 12         29 $self->end_element( { 'Name' => 'Hsp' } );
396             }
397 12 50       34 if ( $self->within_element('hit') ) {
398 12         33 $self->end_element( { 'Name' => 'Hit' } );
399             }
400 12         43 $self->_pushback($buffer);
401 12         49 last;
402             }
403              
404 1192         925 chomp $buffer;
405 1192 100       2194 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         200 $domaincounter{$name}++;
415 127 100       188 if ( $self->within_element('hit') ) {
416 115 50       146 if ( $self->within_element('hsp') ) {
417 115         226 $self->end_element( { 'Name' => 'Hsp' } );
418             }
419 115         331 $self->end_element( { 'Name' => 'Hit' } );
420             }
421              
422 127         185 my $info = [ @{ $hitinfo[ $hitinfo{$name} ] } ];
  127         445  
423 127 50 33     434 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         298 $self->start_element( { 'Name' => 'Hit' } );
441             $self->element(
442             {
443             'Name' => 'Hit_id',
444 127         147 'Data' => shift @{$info}
  127         333  
445             }
446             );
447             $self->element(
448             {
449             'Name' => 'Hit_desc',
450 127         161 'Data' => shift @{$info}
  127         252  
451             }
452             );
453             $self->element(
454             {
455             'Name' => 'Hit_score',
456 127         163 'Data' => shift @{$info}
  127         241  
457             }
458             );
459             $self->element(
460             {
461             'Name' => 'Hit_signif',
462 127         171 'Data' => shift @{$info}
  127         240  
463             }
464             );
465 127         139 my $dom_total = shift @{$info};
  127         136  
466 127 100       76 if (my $hit_end = shift @{$info}) {
  127         191  
467 75         144 $self->element(
468             {
469             'Name' => 'Hit_len',
470             'Data' => $hit_end
471             }
472             );
473             }
474              
475 127         232 $self->start_element( { 'Name' => 'Hsp' } );
476 127         192 my $HSPinfo = shift @hspinfo;
477 127         157 my $id = shift @$HSPinfo;
478              
479 127 50       187 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       329 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
487 4         9 $self->element(
488             {
489             'Name' => 'Hsp_hit-from',
490             'Data' => shift @$HSPinfo
491             }
492             );
493 4         10 $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         10 $self->element(
506             {
507             'Name' => 'Hsp_query-to',
508             'Data' => shift @$HSPinfo
509             }
510             );
511             }
512             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
513 123         286 $self->element(
514             {
515             'Name' => 'Hsp_query-from',
516             'Data' => shift @$HSPinfo
517             }
518             );
519 123         324 $self->element(
520             {
521             'Name' => 'Hsp_query-to',
522             'Data' => shift @$HSPinfo
523             }
524             );
525 123         281 $self->element(
526             {
527             'Name' => 'Hsp_hit-from',
528             'Data' => shift @$HSPinfo
529             }
530             );
531 123         304 $self->element(
532             {
533             'Name' => 'Hsp_hit-to',
534             'Data' => shift @$HSPinfo
535             }
536             );
537             }
538             $self->element(
539             {
540 127         304 'Name' => 'Hsp_score',
541             'Data' => shift @$HSPinfo
542             }
543             );
544 127         324 $self->element(
545             {
546             'Name' => 'Hsp_evalue',
547             'Data' => shift @$HSPinfo
548             }
549             );
550              
551 127 100       367 if ( $domaincounter{$name} == $domaintotal ) {
552 121         504 $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     8831 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         96 $csline = $buffer;
565 104         185 next;
566             }
567             elsif ($buffer =~ /^(\s+ \*->) (\S+)/ox) {
568             # start of domain
569 127         226 $prelength = CORE::length($1);
570 127         129 $width = 0;
571              
572             # deal with fact that start and stop is on same line
573 127         148 my $data = $2;
574 127 100       401 if ($data =~ s/<-?\*?\s*$//)
575             {
576 102         85 $width = CORE::length($data);
577             }
578              
579 127 100       295 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
580 4         13 $self->element(
581             {
582             'Name' => 'Hsp_qseq',
583             'Data' => $data
584             }
585             );
586             }
587             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
588 123         298 $self->element(
589             {
590             'Name' => 'Hsp_hseq',
591             'Data' => $data
592             }
593             );
594             }
595 127         159 $count = 0;
596 127         106 $second_tier = 0;
597             }
598             elsif ($buffer =~ /^(\s+) (\S+) <-?\*? \s*$/ox) {
599             # end of domain
600 25 100       50 $prelength -= 3 unless ( $second_tier++ );
601 25 100       60 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         83 $self->element(
611             {
612             'Name' => 'Hsp_hseq',
613             'Data' => $2
614             }
615             );
616             }
617 25         50 $width = CORE::length($2);
618 25         23 $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         481 next;
626             }
627             elsif ( $count == 0 ) {
628 86 100       182 $prelength -= 3 unless ( $second_tier++ );
629 86 50       132 unless ( defined $prelength ) {
630              
631             # $self->warn("prelength not set");
632 0         0 next;
633             }
634 86 100       187 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
635 20         59 $self->element(
636             {
637             'Name' => 'Hsp_qseq',
638             'Data' => substr( $buffer, $prelength )
639             }
640             );
641             }
642             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
643 66         211 $self->element(
644             {
645             'Name' => 'Hsp_hseq',
646             'Data' => substr( $buffer, $prelength )
647             }
648             );
649             }
650             }
651             elsif ( $count == 1 ) {
652 238 50       311 if ( !defined $prelength ) {
653 0         0 $self->warn("prelength not set");
654             }
655 238 100       252 if ($width) {
656 127         418 $self->element(
657             {
658             'Name' => 'Hsp_midline',
659             'Data' => substr( $buffer, $prelength, $width )
660             }
661             );
662 127 100       258 if ($csline ne '') {
663 54         128 $self->element(
664             {
665             'Name' => 'Hsp_csline',
666             'Data' => substr( $csline, $prelength, $width )
667              
668             }
669             );
670 54         73 $csline = '';
671             }
672             }
673             else {
674 111         344 $self->element(
675             {
676             'Name' => 'Hsp_midline',
677             'Data' => substr( $buffer, $prelength )
678             }
679             );
680 111 100       197 if ($csline ne '') {
681 50         135 $self->element(
682             {
683             'Name' => 'Hsp_csline',
684             'Data' => substr( $csline, $prelength )
685             }
686             );
687 50         68 $csline = '';
688             }
689             }
690             }
691             elsif ( $count == 2 ) {
692 238 50       624 if ( $buffer =~ /^\s+(\S+)\s+(\d+|\-)\s+(\S*)\s+(\d+|\-)/o ) {
693 238 100       452 if ($self->{'_reporttype'} eq 'HMMSEARCH') {
    50          
694 28         105 $self->element(
695             {
696             'Name' => 'Hsp_hseq',
697             'Data' => $3
698             }
699             );
700             }
701             elsif ($self->{'_reporttype'} eq 'HMMPFAM') {
702 210         670 $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       1902 $count = 0 if $count++ >= 2;
715             }
716             }
717             }
718             elsif ( $buffer =~ /^Histogram/o || $buffer =~ m!^//!o ) {
719 16         23 my %domaincounter;
720              
721 16         51 while ( my $HSPinfo = shift @hspinfo ) {
722 4860         5301 my $id = shift @$HSPinfo;
723 4860         8185 $domaincounter{$id}++;
724              
725 4860         3433 my $info = [ @{ $hitinfo[ $hitinfo{$id} ] } ];
  4860         16010  
726 4860 50       7028 next unless defined $info;
727              
728 4860         11186 $self->start_element( { 'Name' => 'Hit' } );
729             $self->element(
730             {
731             'Name' => 'Hit_id',
732 4860         5540 'Data' => shift @{$info}
  4860         10724  
733             }
734             );
735             $self->element(
736             {
737             'Name' => 'Hit_desc',
738 4860         6026 'Data' => shift @{$info}
  4860         9838  
739             }
740             );
741             $self->element(
742             {
743             'Name' => 'Hit_score',
744 4860         5992 'Data' => shift @{$info}
  4860         9341  
745             }
746             );
747             $self->element(
748             {
749             'Name' => 'Hit_signif',
750 4860         6230 'Data' => shift @{$info}
  4860         9172  
751             }
752             );
753 4860         5508 my $domaintotal = shift @{$info};
  4860         4693  
754 4860 100       3246 if (my $hit_end = shift @{$info}) {
  4860         6684  
755 56         139 $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         8384 $self->start_element( { 'Name' => 'Hsp' } );
766 4860         11484 $self->element(
767             {
768             'Name' => 'Hsp_hit-from',
769             'Data' => shift @$HSPinfo
770             }
771             );
772 4860         10962 $self->element(
773             {
774             'Name' => 'Hsp_hit-to',
775             'Data' => shift @$HSPinfo
776             }
777             );
778 4860         11371 $self->element(
779             {
780             'Name' => 'Hsp_query-from',
781             'Data' => shift @$HSPinfo
782             }
783             );
784 4860         12153 $self->element(
785             {
786             'Name' => 'Hsp_query-to',
787             'Data' => shift @$HSPinfo
788             }
789             );
790 4860         11359 $self->element(
791             {
792             'Name' => 'Hsp_score',
793             'Data' => shift @$HSPinfo
794             }
795             );
796 4860         11372 $self->element(
797             {
798             'Name' => 'Hsp_evalue',
799             'Data' => shift @$HSPinfo
800             }
801             );
802 4860         9789 $self->end_element( { 'Name' => 'Hsp' } );
803 4860         13234 $self->end_element( { 'Name' => 'Hit' } );
804              
805 4860 100       18146 if ( $domaincounter{$id} == $domaintotal ) {
806 3004         10565 $hitinfo[ $hitinfo{$id} ] = undef;
807             }
808             }
809 16         216 @hitinfo = ();
810 16         576 %hitinfo = ();
811 16         881 last;
812             }
813             # uncomment to see missed lines with verbose on
814             #else {
815             # $self->debug($buffer);
816             #}
817             }
818 444         737 $last = $buffer;
819             }
820 19 100       90 $self->end_element( { 'Name' => 'HMMER_Output' } ) unless !$seentop;
821 19         54 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 56731 my ( $self, $data ) = @_;
837              
838             # we currently don't care about attributes
839 90966         63042 my $nm = $data->{'Name'};
840 90966         68060 my $type = $MODEMAP{$nm};
841 90966 100       96873 if ($type) {
842 9990 50       17001 if ( $self->_eventHandler->will_handle($type) ) {
843 9990         18755 my $func = sprintf( "start_%s", lc $type );
844 9990         11341 $self->_eventHandler->$func( $data->{'Attributes'} );
845             }
846 9990         10564 unshift @{ $self->{'_elements'} }, $type;
  9990         13451  
847             }
848 90966 100 100     154289 if ( defined $type
849             && $type eq 'result' )
850             {
851 16         23 $self->{'_values'} = {};
852 16         26 $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 64039 my ( $self, $data ) = @_;
869 90966         64766 my $nm = $data->{'Name'};
870 90966         68714 my $type = $MODEMAP{$nm};
871 90966         49234 my $rc;
872              
873 90966 100       109603 if ( $nm eq 'HMMER_program' ) {
874 16 100       62 if ( $self->{'_last_data'} =~ /(HMM\S+)/i ) {
875 15         35 $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       93385 if ( $nm eq 'Hsp' ) {
882 4987         5669 foreach my $line (qw(Hsp_csline Hsp_qseq Hsp_midline Hsp_hseq)) {
883 19948         15886 my $data = $self->{'_last_hspdata'}->{$line};
884 19948 100 100     26182 if ($data && $line eq 'Hsp_hseq') {
885             # replace hmm '.' gap symbol by '-'
886 127         182 $data =~ s/\./-/g;
887             }
888             $self->element(
889             {
890 19948         34522 '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       34755 if ($line eq 'Hsp_midline') {
897 4987 100       4806 if ($data) {
898 127         104 my $length = length $data;
899 127         128 my $identical = ($data =~ tr/a-zA-Z//);
900 127         132 my $positive = ($data =~ tr/+//) + $identical;
901 127         266 $self->element(
902             {
903             'Name' => 'Hsp_align-len',
904             'Data' => $length
905             }
906             );
907 127         292 $self->element(
908             { 'Name' => 'Hsp_identity',
909             'Data' => $identical
910             }
911             );
912 127         259 $self->element(
913             { 'Name' => 'Hsp_positive',
914             'Data' => $positive
915             }
916             );
917             }
918             else {
919 4860         8695 $self->element(
920             { 'Name' => 'Hsp_identity',
921             'Data' => 0
922             }
923             );
924 4860         10256 $self->element(
925             { 'Name' => 'Hsp_positive',
926             'Data' => 0
927             }
928             );
929             }
930             }
931             }
932 4987         5166 $self->{'_last_hspdata'} = {};
933             }
934 90966 100       121261 if ($type) {
    50          
935 9990 50       17972 if ( $self->_eventHandler->will_handle($type) ) {
936 9990         19754 my $func = sprintf( "end_%s", lc $type );
937             $rc = $self->_eventHandler->$func( $self->{'_reporttype'},
938 9990         12403 $self->{'_values'} );
939             }
940 9990         8185 my $lastelem = shift @{ $self->{'_elements'} };
  9990         12871  
941              
942             # Flush corresponding values from the {_values} buffer
943 9990         9957 my $name = uc $type;
944 9990         6162 foreach my $key (keys %{ $self->{_values} }) {
  9990         28418  
945 194851 100       468452 delete $self->{_values}->{$key} if ($key =~ m/^$name-/);
946             }
947             }
948             elsif ( $MAPPING{$nm} ) {
949 80976 50       81632 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         126682 $self->{'_values'}->{ $MAPPING{$nm} } = $self->{'_last_data'};
956             }
957             }
958             else {
959 0         0 $self->debug("unknown nm $nm, ignoring\n");
960             }
961 90966         80685 $self->{'_last_data'} = ''; # remove read data if we are at
962             # end of an element
963 90966 100 100     138667 $self->{'_result'} = $rc if ( defined $type && $type eq 'result' );
964 90966         86083 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 59218 my ( $self, $data ) = @_;
980 80976         73813 $self->start_element($data);
981 80976         75327 $self->characters($data);
982 80976         85835 $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 54415 my ( $self, $data ) = @_;
998              
999 80976 100 100     77635 if ( $self->in_element('hsp')
      100        
1000             && $data->{'Name'} =~ /Hsp\_(?:qseq|hseq|csline|midline)/o
1001             && defined $data->{'Data'} )
1002             {
1003 1253         2288 $self->{'_last_hspdata'}->{ $data->{'Name'} } .= $data->{'Data'};
1004             }
1005 80976 100 100     230762 return unless ( defined $data->{'Data'} && $data->{'Data'} !~ /^\s+$/o );
1006              
1007 61460         59765 $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 193 my ( $self, $name ) = @_;
1025             return 0
1026             if ( !defined $name
1027             || !defined $self->{'_elements'}
1028 254 50 33     706 || scalar @{ $self->{'_elements'} } == 0 );
  254   33     566  
1029 254         176 foreach my $element ( @{ $self->{'_elements'} } ) {
  254         313  
1030 369 100       750 return 1 if ( $element eq $name );
1031             }
1032 12         27 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 55000 my ( $self, $name ) = @_;
1050 80988 50       99154 return 0 if !defined $self->{'_elements'}->[0];
1051 80988         304834 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 15 my ($self) = @_;
1067 19         28 $self->{'_lasttype'} = '';
1068 19         33 $self->{'_values'} = {};
1069 19         29 $self->{'_result'} = undef;
1070 19         29 $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 25 my ($self) = @_;
1086 19         124 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;