?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Bio::DB::Universal: move into its own distribution
[bioperl-live.git] / Bio / SearchIO / fasta.pm
blob0a7435ad7bc1a3b1749e700bfbcef1d9b60ba11f
1 #
2 # BioPerl module for Bio::SearchIO::fasta
3 #
4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
5 #
6 # Cared for by Jason Stajich <jason-at-bioperl.org>
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::fasta - A SearchIO parser for FASTA results
17
18 =head1 SYNOPSIS
19
20 # Do not use this object directly, use it through the SearchIO system
21 use Bio::SearchIO;
22 my $searchio = Bio::SearchIO->new(-format => 'fasta',
23 -file => 'report.FASTA');
24 while( my $result = $searchio->next_result ) {
25 # ... do what you would normally doi with Bio::SearchIO.
26 }
27
28 =head1 DESCRIPTION
29
30 This object contains the event based parsing code for FASTA format
31 reports. It creates L<Bio::Search::HSP::FastaHSP> objects instead of
32 L<Bio::Search::HSP::GenericHSP> for the HSP objects.
33
34 This module will parse -m 9 -d 0 output as well as default m 1 output
35 from FASTA as well as SSEARCH.
36
37 Also see the SearchIO HOWTO:
38 L<http://bioperl.org/howtos/SearchIO_HOWTO.html>.
39
40 =head1 FEEDBACK
41
42 =head2 Mailing Lists
43
44 User feedback is an integral part of the evolution of this and other
45 Bioperl modules. Send your comments and suggestions preferably to
46 the Bioperl mailing list. Your participation is much appreciated.
47
48 bioperl-l@bioperl.org - General discussion
49 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
50
51 =head2 Support
52
53 Please direct usage questions or support issues to the mailing list:
54
55 I<bioperl-l@bioperl.org>
56
57 rather than to the module maintainer directly. Many experienced and
58 reponsive experts will be able look at the problem and quickly
59 address it. Please include a thorough description of the problem
60 with code and data examples if at all possible.
61
62 =head2 Reporting Bugs
63
64 Report bugs to the Bioperl bug tracking system to help us keep track
65 of the bugs and their resolution. Bug reports can be submitted via the
66 web:
67
68 https://github.com/bioperl/bioperl-live/issues
69
70 =head1 AUTHOR - Jason Stajich, Aaron Mackey, William Pearson
71
72 Email jason-at-bioperl.org
73
74 =head1 APPENDIX
75
76 The rest of the documentation details each of the object methods.
77 Internal methods are usually preceded with a _
78
79 =cut
80
81 # Let the code begin...
82
83 package Bio::SearchIO::fasta;
84 use vars qw(%MODEMAP %MAPPING $IDLENGTH);
85 use strict;
86
87 # Object preamble - inherits from Bio::Root::RootI
88
89 use Bio::Factory::ObjectFactory;
90
91 BEGIN {
92
93 # Set IDLENGTH to a new value if you have
94 # compile FASTA with a different ID length
95 # (actually newest FASTA allows the setting of this
96 # via -C parameter, default is 6)
97 $IDLENGTH = 6;
98
99 # mapping of NCBI Blast terms to Bioperl hash keys
100 %MODEMAP = (
101 'FastaOutput' => 'result',
102 'Hit' => 'hit',
103 'Hsp' => 'hsp'
104 );
105
106 # This should really be done more intelligently, like with
107 # XSLT
108
109 %MAPPING = (
110 'Hsp_bit-score' => 'HSP-bits',
111 'Hsp_score' => 'HSP-score',
112 'Hsp_sw-score' => 'HSP-swscore',
113 'Hsp_evalue' => 'HSP-evalue',
114 'Hsp_evalue2' => 'HSP-evalue2',
115 'Hsp_query-from' => 'HSP-query_start',
116 'Hsp_query-to' => 'HSP-query_end',
117 'Hsp_hit-from' => 'HSP-hit_start',
118 'Hsp_hit-to' => 'HSP-hit_end',
119 'Hsp_positive' => 'HSP-conserved',
120 'Hsp_identity' => 'HSP-identical',
121 'Hsp_gaps' => 'HSP-hsp_gaps',
122 'Hsp_hitgaps' => 'HSP-hit_gaps',
123 'Hsp_querygaps' => 'HSP-query_gaps',
124 'Hsp_qseq' => 'HSP-query_seq',
125 'Hsp_hseq' => 'HSP-hit_seq',
126 'Hsp_midline' => 'HSP-homology_seq',
127 'Hsp_align-len' => 'HSP-hsp_length',
128 'Hsp_query-frame' => 'HSP-query_frame',
129 'Hsp_hit-frame' => 'HSP-hit_frame',
130
131 'Hit_id' => 'HIT-name',
132 'Hit_len' => 'HIT-length',
133 'Hit_accession' => 'HIT-accession',
134 'Hit_def' => 'HIT-description',
135 'Hit_signif' => 'HIT-significance',
136 'Hit_score' => 'HIT-score',
137
138 'FastaOutput_program' => 'RESULT-algorithm_name',
139 'FastaOutput_version' => 'RESULT-algorithm_version',
140 'FastaOutput_query-def' => 'RESULT-query_name',
141 'FastaOutput_querydesc' => 'RESULT-query_description',
142 'FastaOutput_query-len' => 'RESULT-query_length',
143 'FastaOutput_db' => 'RESULT-database_name',
144 'FastaOutput_db-len' => 'RESULT-database_entries',
145 'FastaOutput_db-let' => 'RESULT-database_letters',
146
147 'Parameters_matrix' => { 'RESULT-parameters' => 'matrix' },
148 'Parameters_expect' => { 'RESULT-parameters' => 'expect' },
149 'Parameters_include' => { 'RESULT-parameters' => 'include' },
150 'Parameters_sc-match' => { 'RESULT-parameters' => 'match' },
151 'Parameters_sc-mismatch' => { 'RESULT-parameters' => 'mismatch' },
152 'Parameters_gap-open' => { 'RESULT-parameters' => 'gapopen' },
153 'Parameters_gap-ext' => { 'RESULT-parameters' => 'gapext' },
154 'Parameters_word-size' => { 'RESULT-parameters' => 'wordsize' },
155 'Parameters_ktup' => { 'RESULT-parameters' => 'ktup' },
156 'Parameters_filter' => { 'RESULT-parameters' => 'filter' },
157 'Statistics_db-num' => { 'RESULT-statistics' => 'dbentries' },
158 'Statistics_db-len' => { 'RESULT-statistics' => 'dbletters' },
159 'Statistics_hsp-len' => { 'RESULT-statistics' => 'hsplength' },
160 'Statistics_eff-space' => { 'RESULT-statistics' => 'effectivespace' },
161 'Statistics_kappa' => { 'RESULT-statistics' => 'kappa' },
162 'Statistics_lambda' => { 'RESULT-statistics' => 'lambda' },
163 'Statistics_entropy' => { 'RESULT-statistics' => 'entropy' },
164 );
165 }
166
167 use base qw(Bio::SearchIO);
168
169 =head2 new
170
171 Title : new
172 Usage : my $obj = Bio::SearchIO::fasta->new();
173 Function: Builds a new Bio::SearchIO::fasta object
174 Returns : Bio::SearchIO::fasta
175 Args : -idlength - set ID length to something other
176 than the default (6), this is only
177 necessary if you have compiled FASTA
178 with a new default id length to display
179 in the HSP alignment blocks
180
181 =cut
182
183 sub _initialize {
184 my ( $self, @args ) = @_;
185 $self->SUPER::_initialize(@args);
186 return unless @args;
187 my ($idlength) = $self->_rearrange( [qw(IDLENGTH)], @args );
188 $self->idlength( $idlength || $IDLENGTH );
189 $self->_eventHandler->register_factory(
190 'hsp',
191 Bio::Factory::ObjectFactory->new(
192 -type => 'Bio::Search::HSP::FastaHSP',
193 -interface => 'Bio::Search::HSP::HSPI'
194 )
195 );
196 return 1;
197 }
198
199 =head2 next_result
200
201 Title : next_result
202 Usage : my $hit = $searchio->next_result;
203 Function: Returns the next Result from a search
204 Returns : Bio::Search::Result::ResultI object
205 Args : none
206
207 =cut
208
209 sub next_result {
210 my ($self) = @_;
211 local $/ = "\n";
212 local $_;
213
214 my $data = '';
215 my $seentop = 0;
216 my $current_hsp;
217 my $m9HSP = 0;
218 $self->start_document();
219 my @hit_signifs;
220 while ( defined( $_ = $self->_readline ) ) {
221 next if ( !$self->in_element('hsp')
222 && /^\s+$/ ); # skip empty lines
223 if (
224 m/(\S+)\s+searches\s+a\s+(protein\s+or\s+DNA\s+)?sequence/oxi
225 || /(\S+)\s+compares\s+a/
226 || /(\S+)\s+performs\s+a/
227 || /(\S+)\s+produces\s/
228 || /(\S+)\s+finds\s+/ # for lalign, but does not work because no "The best scores are:"
229 || ( m/^\#\s+/ # has a command log line
230 && ( $_ = $self->_readline )
231 && /(\S+)\s+searches\s+a\s+(protein\s+or\s+DNA\s+)?sequence/oxi
232 || /(\S+)\s+compares\s+a/
233 || /(\S+)\s+performs\s+a/
234 || /(\S+)\s+produces\s/
235 || /(\S+)\s+finds\s+/ # for lalign, but does not work because no "The best scores are:"
236 )
237 )
238 {
239 if ($seentop) {
240 $self->_pushback($_);
241 $self->end_element( { 'Name' => 'FastaOutput' } );
242 return $self->end_document();
243 }
244 $self->{'_reporttype'} = $1;
245 $self->start_element( { 'Name' => 'FastaOutput' } );
246 $self->{'_result_count'}++;
247 $seentop = 1;
248 #$self->debug( "reporttype is " . $self->{'_reporttype'} . "\n" );
249 $self->element(
250 {
251 'Name' => 'FastaOutput_program',
252 'Data' => $self->{'_reporttype'}
253 }
254 );
255 my $version;
256 # version 35 version string on same line
257 if (/version/) {
258 ($version) = (/version\s+(\S+)/);
259 }
260 # earlier versions, it's on the next line
261 else {
262 $_ = $self->_readline();
263 ($version) = (/version\s+(\S+)/);
264 }
265 $version = '' unless defined $version;
266 $self->{'_version'} = $version;
267 $self->element(
268 {
269 'Name' => 'FastaOutput_version',
270 'Data' => $version
271 }
272 );
273
274 my ( $last, $leadin, $type, $querylen, $querytype, $querydef );
275
276 while ( defined( $_ = $self->_readline() ) ) {
277 if (
278 /^ (
279 (?:\s+>) | # fa33 lead-in
280 (?:\s*\d+\s*>>>) # fa34 mlib lead-in
281 )
282 (.*)
283 /x
284 )
285 {
286 ( $leadin, $querydef ) = ( $1, $2 );
287 if ( $leadin =~ m/>>>/ ) {
288 if ( $querydef =~
289 /^(.*?)\s+(?:\-\s+)?(\d+)\s+(aa|nt).*$/o )
290 {
291 ( $querydef, $querylen, $querytype ) =
292 ( $1, $2, $3 );
293 last;
294 }
295 }
296 else {
297 if ( $last =~ /(\S+)[:,]\s*(\d+)\s+(aa|nt)/ ) {
298 ( $querylen, $querytype ) = ( $2, $3 );
299 $querydef ||= $1;
300 last;
301 }
302 }
303 }
304 elsif (m/^\s*vs\s+\S+/o) {
305 if ( $last =~ /(\S+)[,:]\s+(\d+)\s+(aa|nt)/o ) {
306 ( $querydef, $querylen, $querytype ) = ( $1, $2, $3 );
307 last;
308 }
309 }
310 $last = $_;
311 }
312 if ( $self->{'_reporttype'}
313 && $self->{'_reporttype'} eq 'FASTA' )
314 {
315 if ( $querytype eq 'nt' ) {
316 $self->{'_reporttype'} = 'FASTN';
317 }
318 elsif ( $querytype eq 'aa' ) {
319 $self->{'_reporttype'} = 'FASTP';
320 }
321 }
322 my ( $name, $descr ) = $querydef =~ m/^(\S+)\s*(.*?)\s*$/o;
323 $self->element(
324 {
325 'Name' => 'FastaOutput_query-def',
326 'Data' => $name
327 }
328 );
329 $self->element(
330 {
331 'Name' => 'FastaOutput_querydesc',
332 'Data' => $descr
333 }
334 );
335 if ($querylen) {
336 $self->element(
337 {
338 'Name' => 'FastaOutput_query-len',
339 'Data' => $querylen
340 }
341 );
342 }
343 else {
344 $self->warn("unable to find and set query length");
345 }
346 if (
347 $last =~ /^\s*vs\s+(\S+)/
348 || ( $last =~ /^searching\s+(\S+)\s+library/ )
349 || ( $last =~ /^Library:\s+(\S+)\s+/ )
350 || (
351 defined $_
352 && ( /^\s*vs\s+(\S+)/
353 || /^Library:\s+(\S+)\s+/ )
354 )
355 || ( defined( $_ = $self->_readline() )
356 && ( /^\s*vs\s+(\S+)/ || /^Library:\s+(\S+)/ ) )
357 )
358 {
359 $self->element(
360 {
361 'Name' => 'FastaOutput_db',
362 'Data' => $1
363 }
364 );
365 }
366 elsif (m/^\s+opt(?:\s+E\(\))?$/o) {
367
368 # histogram ... read over it more rapidly than the larger outer loop:
369 while ( defined( $_ = $self->_readline ) ) {
370 last if m/^>\d+/;
371 }
372 }
373 }
374 elsif (/(\d+)\s+residues\s+in\s+(\d+)\s+(?:library\s+)?sequences/) {
375 $self->element(
376 {
377 'Name' => 'FastaOutput_db-let',
378 'Data' => $1
379 }
380 );
381 $self->element(
382 {
383 'Name' => 'FastaOutput_db-len',
384 'Data' => $2
385 }
386 );
387 $self->element(
388 {
389 'Name' => 'Statistics_db-len',
390 'Data' => $1
391 }
392 );
393 $self->element(
394 {
395 'Name' => 'Statistics_db-num',
396 'Data' => $2
397 }
398 );
399 }
400 elsif (/Lambda=\s*(\S+)/) {
401 $self->element(
402 {
403 'Name' => 'Statistics_lambda',
404 'Data' => $1
405 }
406 );
407 }
408 elsif (/K=\s*(\S+)/) {
409 $self->element(
410 {
411 'Name' => 'Statistics_kappa',
412 'Data' => $1
413 }
414 );
415 }
416 elsif (/^\s*(Smith-Waterman)/) {
417
418 $self->{'_reporttype'} = $1;
419
420 m/\[\s*(\S+)\s+matrix \([^\)]+\)(xS)?\],/;
421
422 $self->element(
423 {
424 'Name' => 'Parameters_matrix',
425 'Data' => $1
426 }
427 );
428 $self->element(
429 {
430 'Name' => 'Parameters_filter',
431 'Data' => defined $2 ? 1 : 0,
432 }
433 );
434 if (/\s*gap\-penalty:\s*(\-?\d+)\/(\-?\d+)/) {
435 $self->element(
436 {
437 'Name' => 'Parameters_gap-open',
438 'Data' => $1,
439 }
440 );
441
442 $self->element(
443 {
444 'Name' => 'Parameters_gap-ext',
445 'Data' => $2,
446 }
447 );
448 }
449 elsif (/\s*open\/ext:\s*(\-?\d+)\/(\-?\d+)/) {
450 $self->element(
451 {
452 'Name' => 'Parameters_gap-open',
453 'Data' => $1,
454 }
455 );
456
457 $self->element(
458 {
459 'Name' => 'Parameters_gap-ext',
460 'Data' => $2,
461 }
462 );
463 }
464
465 $self->element(
466 {
467 'Name' => 'FastaOutput_program',
468 'Data' => $self->{'_reporttype'}
469 }
470 );
471 }
472 elsif (/The best( related| unrelated)? scores are:/) {
473 my $rel = $1;
474 my @labels = split;
475 @labels = map {
476 if ( $_ =~ m/^E\((\d+)\)$/o )
477 {
478 $self->element(
479 { 'Name' => 'Statistics_eff-space', 'Data' => $1 } );
480 "evalue";
481 }
482 else {
483 # canonicalize changed column headers
484 if ($_ eq "gapl") {
485 $_ = "lgaps";
486 } elsif ($_ eq "gapq") {
487 $_ = "qgaps";
488 } elsif ($_ eq "E2()") {
489 $_ = "evalue2";
490 }
491
492 $_;
493 }
494 } @labels[ $rel ? 5 : 4 .. $#labels ];
495
496 while ( defined( $_ = $self->_readline() )
497 && !/^\s+$/ )
498 {
499 my @line = split;
500
501 if ( $line[-1] =~ m/\=/o && $labels[-1] ne 'aln_code' ) {
502 # unlabelled alignment hit;
503 push @labels, "aln_code";
504 }
505
506 if ($line[0] eq "+-") {
507 $m9HSP = 1;
508 # parse HSP, add to last parsed Hit
509 my %hspData;
510
511 my @line = split;
512
513 @hspData{@labels} = splice( @line, @line - @labels );
514 $hspData{lframe} = $hit_signifs[-1]->{lframe};
515
516 push @{$hit_signifs[-1]->{HSPs}}, \%hspData;
517
518 next;
519 }
520 elsif ($line[0] eq '>>><<<') {
521 last;
522 }
523
524 my (%data, %hspData);
525 @data{@labels} = @hspData{@labels} = splice( @line, @line - @labels );
526 if ( $line[-1] =~ m/\[([1-6rf])\]/o ) {
527 my $fr = $1;
528 $hspData{lframe} = $data{lframe} = (
529 $fr =~ /\d/o
530 ? ( $fr <= 3 ? "+$fr" : "-@{[$fr-3]}" )
531 : ( $fr eq 'f' ? '+1' : '-1' )
532 );
533 pop @line;
534 }
535 else {
536 $hspData{lframe} = $data{lframe} = '0';
537 }
538
539 if ( $line[-1] =~ m/^\(?(\d+)\)$/ ) {
540 $data{hit_len} = $1;
541 pop @line;
542 if ( $line[-1] =~ m/^\($/ ) {
543 pop @line;
544 }
545 }
546 else {
547 $data{hit_len} = 0;
548 }
549
550 # rebuild the first part of the line, preserving spaces:
551 ($_) = m/^(\S+(?:\s+\S+){$#line})/;
552
553 my ( $id, $desc ) = split( /\s+/, $_, 2 );
554 my @pieces = split( /\|/, $id );
555 my $acc = pop @pieces;
556 $acc =~ s/\.\d+$//;
557
558 @data{qw(id desc acc)} = ( $id, $desc, $acc );
559 $data{HSPs} = [ \%hspData ];
560
561 push @hit_signifs, \%data;
562 }
563 }
564 elsif (
565 /^\s*([T]?FAST[XYAF]).+,\s*(\S+)\s*matrix[^\]]+?(xS)?\]\s*ktup:\s*(\d+)/
566 )
567 {
568
569 $self->element(
570 {
571 'Name' => 'Parameters_matrix',
572 'Data' => $2
573 }
574 );
575 $self->element(
576 {
577 'Name' => 'Parameters_filter',
578 'Data' => defined $3 ? 1 : 0,
579 }
580 );
581 $self->element(
582 {
583 'Name' => 'Parameters_ktup',
584 'Data' => $4
585 }
586 );
587 $self->{'_reporttype'} = $1
588 if ( $self->{'_reporttype'} !~ /FAST[PN]/i );
589
590 #
591 # get gap-pen line for FASTA33, which is not on the matrix line
592 #
593 # FASTA (3.36 June 2000) function [optimized, BL50 matrix (15:-5)] ktup: 2
594 # join: 36, opt: 24, gap-pen: -12/ -2, width: 16
595 #
596 $_ = $self->_readline();
597 if (/(?:gap\-pen|open\/ext):\s+([\-\+]?\d+)\s*\/\s*([\-\+]?\d+)/) {
598 $self->element(
599 {
600 'Name' => 'Parameters_gap-open',
601 'Data' => $1
602 }
603 );
604 $self->element(
605 {
606 'Name' => 'Parameters_gap-ext',
607 'Data' => $2
608 }
609 );
610 }
611
612 $self->element(
613 {
614 'Name' => 'FastaOutput_program',
615 'Data' => $self->{'_reporttype'}
616 }
617 );
618
619 }
620 elsif (/^Algorithm:\s+(\S+)\s+.*\s*\(([^)]+)\)\s+(\S+)/) {
621 $self->{'_reporttype'} = $1
622 if ( $self->{'_reporttype'} !~ /FAST[PN]/i );
623 }
624 elsif ( /^Parameters:/ ) { # FASTA 35.04/FASTA 36
625 m/Parameters:\s+(\S+)\s+matrix\s+\([^\)]+\)(xS)?,?\s/;
626 $self->element(
627 {
628 'Name' => 'Parameters_matrix',
629 'Data' => $1
630 }
631 );
632 $self->element(
633 {
634 'Name' => 'Parameters_filter',
635 'Data' => defined $2 ? $2 : 0,
636 }
637 );
638 if (/ktup:\s(\d+)/) {
639 $self->element(
640 {
641 'Name' => 'Parameters_ktup',
642 'Data' => $1
643 }
644 );
645 if (/ktup:\s\d+$/) {
646 $_ = $self->_readline();
647 }
648
649 }
650 if (/(?:gap\-pen|open\/ext):\s+([\-\+]?\d+)\s*\/\s*([\-\+]?\d+)/) {
651 $self->element(
652 {
653 'Name' => 'Parameters_gap-open',
654 'Data' => $1
655 }
656 );
657 $self->element(
658 {
659 'Name' => 'Parameters_gap-ext',
660 'Data' => $2
661 }
662 );
663 }
664 $self->element(
665 {
666 'Name' => 'FastaOutput_program',
667 'Data' => $self->{'_reporttype'}
668 }
669 );
670 }
671 elsif (
672 /^\s+ktup:\s*(\d+),/
673 )
674 {
675 $self->element(
676 {
677 'Name' => 'Parameters_ktup',
678 'Data' => $1
679 }
680 );
681 }
682 elsif (/^(>--)$/ || /^>>(?!>)(.+?)\s+(?:\((\d+)\s*(aa|nt)\))?$/) {
683 if ( $self->in_element('hsp') ) {
684 $self->end_element( { 'Name' => 'Hsp' } );
685 }
686
687 my $firstHSP = 0;
688 if ($1 ne ">--") {
689 $firstHSP = 1;
690
691 my ($hit_id, $len, $alphabet) = ($1, $2, $3);
692 if (!$len || !$alphabet) {
693 WRAPPED:
694 while (defined($_ = $self->_readline)) {
695 if (/(.*?)\s+\((\d+)\s*(aa|nt)\)/) {
696 ($len, $alphabet) = ($2, $3);
697 $hit_id .= $1 ? " ".$1 : '';
698 last WRAPPED;
699 }
700 if (/^>>(?!>)/) { # too far, throw
701 $self->throw("Couldn't find length, bailing");
702 }
703 }
704 }
705 if ( $self->in_element('hit') ) {
706 $self->end_element( { 'Name' => 'Hit' } );
707 shift @hit_signifs if @hit_signifs;
708 }
709
710 $self->start_element( { 'Name' => 'Hit' } );
711 $self->element(
712 {
713 'Name' => 'Hit_len',
714 'Data' => $len
715 }
716 );
717 my ( $id, $desc ) = split( /\s+/, $hit_id, 2 );
718 $self->element(
719 {
720 'Name' => 'Hit_id',
721 'Data' => $id
722 }
723 );
724
725 #$self->debug("Hit ID is $id\n");
726 my @pieces = split( /\|/, $id );
727 my $acc = pop @pieces;
728 $acc =~ s/\.\d+$//;
729 $self->element(
730 {
731 'Name' => 'Hit_accession',
732 'Data' => $acc
733 }
734 );
735 $self->element(
736 {
737 'Name' => 'Hit_def',
738 'Data' => $desc
739 }
740 );
741 }
742 else {
743 # push @{$hit_signifs[0]->{HSPs}}, $current_hsp;
744 }
745
746
747 $_ = $self->_readline();
748 my ( $score, $bits, $e, $e2 ) = /Z-score: \s* (\S+) \s*
749 (?: bits: \s* (\S+) \s+ )?
750 (?: E|expect ) \s* \((?:\d+)?\) :? \s*(\S+)
751 (?: \s* E2 \s* \(\) :? \s*(\S+) )?
752 /ox;
753 $bits = $score unless defined $bits;
754
755 my ($v);
756
757 if ($firstHSP && !$m9HSP) {
758 $v = shift @{$hit_signifs[0]->{HSPs}}
759 if (@hit_signifs && @{$hit_signifs[0]->{HSPs}});
760 $current_hsp = $v;
761 }
762 else {
763 $v = $current_hsp;
764 }
765
766 if ( defined $v ) {
767 @{$v}{qw(evalue evalue2 bits z-sc)} = ( $e, $e2, $bits, $score );
768 }
769
770 if ($firstHSP) {
771 $self->element(
772 {
773 'Name' => 'Hit_signif',
774 'Data' => $v ? $v->{evalue} : $e
775 }
776 );
777 $self->element(
778 {
779 'Name' => 'Hit_score',
780 'Data' => $v ? $v->{bits} : $bits
781 }
782 );
783 }
784
785 $self->start_element( { 'Name' => 'Hsp' } );
786
787 $self->element(
788 {
789 'Name' => 'Hsp_score',
790 'Data' => $v ? $v->{'z-sc'} : $score
791 }
792 );
793 $self->element(
794 {
795 'Name' => 'Hsp_evalue',
796 'Data' => $v ? $v->{evalue} : $e
797 }
798 );
799 $self->element(
800 {
801 'Name' => 'Hsp_evalue2',
802 'Data' => $v && exists($v->{evalue2}) ? $v->{evalue2} : $e2
803 }
804 ) if (($v && exists($v->{evalue2})) || defined $e2);
805
806 $self->element(
807 {
808 'Name' => 'Hsp_bit-score',
809 'Data' => $v ? $v->{bits} : $bits
810 }
811 );
812 $_ = $self->_readline();
813
814 if (s/global\/.* score:\s*(\d+)\;?//) {
815 $self->element(
816 {
817 'Name' => 'Hsp_sw-score',
818 'Data' => $1
819 }
820 );
821 }
822 if (s/Smith-Waterman score:\s*(\d+)\;?//) {
823 $self->element(
824 {
825 'Name' => 'Hsp_sw-score',
826 'Data' => $1
827 }
828 );
829 }
830 if (
831 / (\d*\.?\d+)\% \s* identity
832 (?:\s* \(\s*(\S+)\% \s* (?:ungapped|similar) \) )?
833 \s* in \s* (\d+) \s+ (?:aa|nt) \s+ overlap \s*
834 \( (\d+) \- (\d+) : (\d+) \- (\d+) \)
835 /x
836 )
837 {
838 my ( $identper, $gapper, $len, $querystart, $queryend,
839 $hitstart, $hitend )
840 = ( $1, $2, $3, $4, $5, $6, $7 );
841 my $ident = sprintf( "%.0f", ( $identper / 100 ) * $len );
842 my $positive = sprintf( "%.0f", ( $gapper / 100 ) * $len );
843
844 $self->element(
845 {
846 'Name' => 'Hsp_identity',
847 'Data' => $ident
848 }
849 );
850 $self->element(
851 {
852 'Name' => 'Hsp_positive',
853 'Data' => $positive
854 }
855 );
856 $self->element(
857 {
858 'Name' => 'Hsp_align-len',
859 'Data' => $len
860 }
861 );
862
863 $self->element(
864 {
865 'Name' => 'Hsp_query-from',
866 'Data' => $querystart
867 }
868 );
869 $self->element(
870 {
871 'Name' => 'Hsp_query-to',
872 'Data' => $queryend
873 }
874 );
875 $self->element(
876 {
877 'Name' => 'Hsp_hit-from',
878 'Data' => $hitstart
879 }
880 );
881 $self->element(
882 {
883 'Name' => 'Hsp_hit-to',
884 'Data' => $hitend
885 }
886 );
887
888 }
889
890 if ($v) {
891 $self->element(
892 { 'Name' => 'Hsp_querygaps', 'Data' => $v->{qgaps} } )
893 if exists $v->{qgaps};
894 $self->element(
895 { 'Name' => 'Hsp_hitgaps', 'Data' => $v->{lgaps} } )
896 if exists $v->{lgaps};
897
898 if ( $self->{'_reporttype'} =~ m/^FAST[NXY]$/o ) {
899 if ( 8 == scalar grep { exists $v->{$_} }
900 qw(an0 ax0 pn0 px0 an1 ax1 pn1 px1) )
901 {
902 if ( $v->{ax0} < $v->{an0} ) {
903 $self->element(
904 {
905 'Name' => 'Hsp_query-frame',
906 'Data' =>
907 "-@{[(($v->{px0} - $v->{ax0}) % 3) + 1]}"
908 }
909 );
910 }
911 else {
912 $self->element(
913 {
914 'Name' => 'Hsp_query-frame',
915 'Data' =>
916 "+@{[(($v->{an0} - $v->{pn0}) % 3) + 1]}"
917 }
918 );
919 }
920 if ( $v->{ax1} < $v->{an1} ) {
921 $self->element(
922 {
923 'Name' => 'Hsp_hit-frame',
924 'Data' =>
925 "-@{[(($v->{px1} - $v->{ax1}) % 3) + 1]}"
926 }
927 );
928 }
929 else {
930 $self->element(
931 {
932 'Name' => 'Hsp_hit-frame',
933 'Data' =>
934 "+@{[(($v->{an1} - $v->{pn1}) % 3) + 1]}"
935 }
936 );
937 }
938 }
939 else {
940 $self->element(
941 {
942 'Name' => 'Hsp_query-frame',
943 'Data' => $v->{lframe}
944 }
945 );
946 $self->element(
947 { 'Name' => 'Hsp_hit-frame', 'Data' => 0 } );
948 }
949 }
950 else {
951 $self->element(
952 { 'Name' => 'Hsp_query-frame', 'Data' => 0 } );
953 $self->element(
954 { 'Name' => 'Hsp_hit-frame', 'Data' => $v->{lframe} } );
955 }
956
957 }
958 else {
959 $self->warn("unable to parse FASTA score line: $_");
960 }
961 }
962 elsif (/\d+\s*residues\s*in\s*\d+\s*query\s*sequences/) {
963 if ( $self->in_element('hsp') ) {
964 $self->end_element( { 'Name' => 'Hsp' } );
965 }
966 if ( $self->in_element('hit') ) {
967 $self->end_element( { 'Name' => 'Hit' } );
968 shift @hit_signifs if @hit_signifs;
969
970 }
971
972 # $_ = $self->_readline();
973 # my ( $liblen,$libsize) = /(\d+)\s+residues\s*in(\d+)\s*library/;
974 # fast forward to the end of the file as there is
975 # nothing else left to do with this file and want to be sure and
976 # reset it
977 while ( defined( $_ = $self->_readline() ) ) {
978 last if (/^Function used was/);
979 if (
980 /(\S+)\s+searches\s+a\s+(protein\s+or\s+DNA\s+)?
981 sequence/oxi || /(\S+)\s+compares\s+a/oi
982 )
983 {
984 $self->_pushback($_);
985 }
986 }
987
988 $self->_processHits(@hit_signifs) if @hit_signifs;
989
990 $self->end_element( { 'Name' => 'FastaOutput' } );
991 return $self->end_document();
992 }
993 elsif (/^\s*\d+\s*>>>/) {
994 if ( $self->within_element('FastaOutput') ) {
995 if ( $self->in_element('hsp') ) {
996 $self->end_element( { 'Name' => 'Hsp' } );
997 }
998 if ( $self->in_element('hit') ) {
999 $self->end_element( { 'Name' => 'Hit' } );
1000 shift @hit_signifs if @hit_signifs;
1001
1002 }
1003
1004 $self->_processHits(@hit_signifs) if (@hit_signifs);
1005
1006 $self->end_element( { 'Name' => 'FastaOutput' } );
1007 $self->_pushback($_);
1008 return $self->end_document();
1009 }
1010 else {
1011 $self->start_element( { 'Name' => 'FastaOutput' } );
1012 $self->{'_result_count'}++;
1013 $seentop = 1;
1014 $self->element(
1015 {
1016 'Name' => 'FastaOutput_program',
1017 'Data' => $self->{'_reporttype'}
1018 }
1019 );
1020 $self->element(
1021 {
1022 'Name' => 'FastaOutput_version',
1023 'Data' => $self->{'_version'}
1024 }
1025 );
1026
1027 my ( $type, $querylen, $querytype, $querydef );
1028
1029 if (/^\s*\d+\s*>>>(.*)/) {
1030 $querydef = $1;
1031 if ( $querydef =~ /^(.*?)\s+(?:\-\s+)?(\d+)\s+(aa|nt).*$/o )
1032 {
1033 ( $querydef, $querylen, $querytype ) = ( $1, $2, $3 );
1034 }
1035 }
1036
1037 if ( $self->{'_reporttype'}
1038 && $self->{'_reporttype'} eq 'FASTA' )
1039 {
1040 if ( $querytype eq 'nt' ) {
1041 $self->{'_reporttype'} = 'FASTN';
1042 }
1043 elsif ( $querytype eq 'aa' ) {
1044 $self->{'_reporttype'} = 'FASTP';
1045 }
1046 }
1047 my ( $name, $descr ) =
1048 ( $querydef =~ m/^(\S+)(?:\s+(.*))?\s*$/o );
1049 $self->element(
1050 {
1051 'Name' => 'FastaOutput_query-def',
1052 'Data' => $name
1053 }
1054 );
1055 $self->element(
1056 {
1057 'Name' => 'FastaOutput_querydesc',
1058 'Data' => $descr
1059 }
1060 );
1061 if ($querylen) {
1062 $self->element(
1063 {
1064 'Name' => 'FastaOutput_query-len',
1065 'Data' => $querylen
1066 }
1067 );
1068 }
1069 else {
1070 $self->warn("unable to find and set query length");
1071 }
1072 if ( defined( $_ = $self->_readline() )
1073 && ( /^\s*vs\s+(\S+)/ || /^Library:\s+(\S+)/ ) )
1074 {
1075 $self->element(
1076 {
1077 'Name' => 'FastaOutput_db',
1078 'Data' => $1
1079 }
1080 );
1081 }
1082
1083 }
1084 }
1085 elsif ( $self->in_element('hsp') ) {
1086 my @data = ( [], [], [] );
1087 my $count = 0;
1088 my $len = $self->idlength + 1;
1089 my ($seq1_id);
1090 while ( defined($_) ) {
1091 chomp;
1092 #$self->debug("$count $_\n");
1093 if (/residues in \d+\s+query\s+sequences/o) {
1094 $self->_pushback($_);
1095 last;
1096 }
1097 elsif (/^>>>(\*\*\*|\/\/\/|<<<)/o) {
1098 $self->end_element( { Name => "Hsp" } );
1099 last;
1100 }
1101 elsif (/^>>/o) {
1102 $self->_pushback($_);
1103 last;
1104 }
1105 elsif (/^\s*\d+\s*>>>/o) {
1106 $self->_pushback($_);
1107 last;
1108 }
1109 if ( $count == 0 ) {
1110 if (/^(\S+)\s+/) {
1111 $self->_pushback($_);
1112 $count = 2;
1113 }
1114 elsif ( /^\s+\d+/ || /^\s+$/ ) {
1115
1116 # do nothing, this is really a 0 line
1117 }
1118 elsif ( length($_) == 0 ) {
1119 $count = -1;
1120 }
1121 else {
1122 $self->_pushback($_);
1123 $count = 0;
1124 }
1125 }
1126 elsif ( $count == 1 || $count == 3 ) {
1127 if (/^(\S+)\s+/) {
1128 $len = CORE::length($1) if $len < CORE::length($1);
1129 s/\s+$//; # trim trailing spaces,we don't want them
1130 push @{ $data[ $count - 1 ] }, substr( $_, $len );
1131 }
1132 elsif (/^\s+(\d+)/) {
1133 $count = -1;
1134 $self->_pushback($_);
1135 }
1136 elsif ( /^\s+$/ || length($_) == 0 ) {
1137 $count = 5;
1138
1139 # going to skip these
1140 }
1141 elsif ( /\s+\S+fasta3\d\s+/) {
1142 # this is something that looks like a path but contains
1143 # the fasta3x executable string, such as:
1144 # /usr/local/fasta3/bin/fasta35 -n -U -Q -H -A -E 2.0 -C 19 -m 0 -m 9i test.fa ../other_mirs.fa -O test.fasta35
1145 last;
1146 }
1147 else {
1148 $self->throw(
1149 "Unrecognized alignment line ($count) '$_'");
1150 }
1151 }
1152 elsif ( $count == 2 ) {
1153 if (/^\s+\d+\s+/) {
1154 $self->warn("$_\n") if $self->verbose > 0;
1155
1156 # we are on a Subject part of the alignment
1157 # but we THOUGHT we were on the Query
1158 # move that last line to the proper place
1159 push @{ $data[2] }, pop @{ $data[0] };
1160 $count = 4;
1161 }
1162 else {
1163
1164 # toss the first IDLENGTH characters of the line
1165 if ( length($_) >= $len ) {
1166 push @{ $data[ $count - 1 ] }, substr( $_, $len );
1167 }
1168 elsif (/^$/) { # fix issue #255
1169 $count = 5;
1170 }
1171 }
1172 }
1173 last if ( $count++ >= 5 );
1174 $_ = $self->_readline();
1175 }
1176 if ( @{ $data[0] } || @{ $data[2] } ) {
1177 $self->characters(
1178 {
1179 'Name' => 'Hsp_qseq',
1180 'Data' => join( '', @{ $data[0] } )
1181 }
1182 );
1183 $self->characters(
1184 {
1185 'Name' => 'Hsp_midline',
1186 'Data' => join( '', @{ $data[1] } )
1187 }
1188 );
1189 $self->characters(
1190 {
1191 'Name' => 'Hsp_hseq',
1192 'Data' => join( '', @{ $data[2] } )
1193 }
1194 );
1195 }
1196 }
1197 else {
1198 if ( !$seentop ) {
1199 $self->debug($_);
1200 #$self->warn("unrecognized FASTA Family report file!");
1201 #return;
1202 }
1203 }
1204 }
1205 if ( $self->in_element('result') ) {
1206 if ( $self->in_element('hsp') ) {
1207 $self->end_element( { 'Name' => 'Hsp' } );
1208 }
1209 if ( $self->in_element('hit') ) {
1210 $self->end_element( { 'Name' => 'Hit' } );
1211 shift @hit_signifs if @hit_signifs;
1212 }
1213 $self->end_element( { 'Name' => 'FastaOutput' } );
1214 }
1215 return $self->end_document();
1216 }
1217
1218 =head2 start_element
1219
1220 Title : start_element
1221 Usage : $eventgenerator->start_element
1222 Function: Handles a start element event
1223 Returns : none
1224 Args : hashref with at least 2 keys 'Data' and 'Name'
1225
1226
1227 =cut
1228
1229 sub start_element {
1230 my ( $self, $data ) = @_;
1231
1232 # we currently don't care about attributes
1233 my $nm = $data->{'Name'};
1234 if ( my $type = $MODEMAP{$nm} ) {
1235 $self->_mode($type);
1236 if ( my $handler = $self->_will_handle($type) ) {
1237 my $func = sprintf( "start_%s", lc $type );
1238 $handler->$func( $data->{'Attributes'} );
1239 }
1240 unshift @{ $self->{'_elements'} }, $type;
1241 }
1242 if ( $nm eq 'FastaOutput' ) {
1243 $self->{'_values'} = {};
1244 $self->{'_result'} = undef;
1245 $self->{'_mode'} = '';
1246 }
1247
1248 }
1249
1250 =head2 end_element
1251
1252 Title : start_element
1253 Usage : $eventgenerator->end_element
1254 Function: Handles an end element event
1255 Returns : none
1256 Args : hashref with at least 2 keys 'Data' and 'Name'
1257
1258
1259 =cut
1260
1261 sub end_element {
1262 my ( $self, $data ) = @_;
1263 my $nm = $data->{'Name'};
1264 my $rc;
1265
1266 # Hsp are sort of weird, in that they end when another
1267 # object begins so have to detect this in end_element for now
1268 if ( $nm eq 'Hsp' ) {
1269 foreach (qw(Hsp_qseq Hsp_midline Hsp_hseq)) {
1270 $self->element(
1271 {
1272 'Name' => $_,
1273 'Data' => $self->{'_last_hspdata'}->{$_}
1274 }
1275 );
1276 }
1277 $self->{'_last_hspdata'} = {};
1278 }
1279
1280 if ( my $type = $MODEMAP{$nm} ) {
1281 if ( my $handler = $self->_will_handle($type) ) {
1282 my $func = sprintf( "end_%s", lc $type );
1283 $rc = $handler->$func( $self->{'_reporttype'}, $self->{'_values'} );
1284 }
1285 shift @{ $self->{'_elements'} };
1286
1287 }
1288 elsif ( $MAPPING{$nm} ) {
1289 if ( ref( $MAPPING{$nm} ) =~ /hash/i ) {
1290 my $key = ( keys %{ $MAPPING{$nm} } )[0];
1291 $self->{'_values'}->{$key}->{ $MAPPING{$nm}->{$key} } =
1292 $self->{'_last_data'};
1293 }
1294 else {
1295 $self->{'_values'}->{ $MAPPING{$nm} } = $self->{'_last_data'};
1296 }
1297 }
1298 else {
1299 $self->warn("unknown nm $nm, ignoring\n");
1300 }
1301 $self->{'_last_data'} = ''; # remove read data if we are at
1302 # end of an element
1303 $self->{'_result'} = $rc if ( $nm eq 'FastaOutput' );
1304 return $rc;
1305
1306 }
1307
1308 =head2 element
1309
1310 Title : element
1311 Usage : $eventhandler->element({'Name' => $name, 'Data' => $str});
1312 Function: Convience method that calls start_element, characters, end_element
1313 Returns : none
1314 Args : Hash ref with the keys 'Name' and 'Data'
1315
1316
1317 =cut
1318
1319 sub element {
1320 my ( $self, $data ) = @_;
1321 $self->start_element($data);
1322 $self->characters($data);
1323 $self->end_element($data);
1324 }
1325
1326 =head2 characters
1327
1328 Title : characters
1329 Usage : $eventgenerator->characters($str)
1330 Function: Send a character events
1331 Returns : none
1332 Args : string
1333
1334
1335 =cut
1336
1337 sub characters {
1338 my ( $self, $data ) = @_;
1339
1340 return unless ( defined $data->{'Data'} );
1341 if ( $data->{'Data'} =~ /^\s+$/ ) {
1342 return unless $data->{'Name'} =~ /Hsp\_(midline|qseq|hseq)/;
1343 }
1344
1345 if ( $self->in_element('hsp')
1346 && $data->{'Name'} =~ /Hsp\_(qseq|hseq|midline)/ )
1347 {
1348
1349 $self->{'_last_hspdata'}->{ $data->{'Name'} } .= $data->{'Data'};
1350 }
1351
1352 $self->{'_last_data'} = $data->{'Data'};
1353 }
1354
1355 =head2 _mode
1356
1357 Title : _mode
1358 Usage : $obj->_mode($newval)
1359 Function:
1360 Example :
1361 Returns : value of _mode
1362 Args : newvalue (optional)
1363
1364
1365 =cut
1366
1367 sub _mode {
1368 my ( $self, $value ) = @_;
1369 if ( defined $value ) {
1370 $self->{'_mode'} = $value;
1371 }
1372 return $self->{'_mode'};
1373 }
1374
1375 =head2 within_element
1376
1377 Title : within_element
1378 Usage : if( $eventgenerator->within_element($element) ) {}
1379 Function: Test if we are within a particular element
1380 This is different than 'in' because within can be tested
1381 for a whole block.
1382 Returns : boolean
1383 Args : string element name
1384
1385
1386 =cut
1387
1388 sub within_element {
1389 my ( $self, $name ) = @_;
1390 return 0
1391 if (!defined $name && !defined $self->{'_elements'}
1392 || scalar @{ $self->{'_elements'} } == 0 );
1393 foreach ( @{ $self->{'_elements'} } ) {
1394 if ( $_ eq $name || $_ eq $MODEMAP{$name} ) {
1395 return 1;
1396 }
1397 }
1398 return 0;
1399 }
1400
1401 =head2 in_element
1402
1403 Title : in_element
1404 Usage : if( $eventgenerator->in_element($element) ) {}
1405 Function: Test if we are in a particular element
1406 This is different than 'in' because within can be tested
1407 for a whole block.
1408 Returns : boolean
1409 Args : string element name
1410
1411
1412 =cut
1413
1414 sub in_element {
1415 my ( $self, $name ) = @_;
1416 return 0 if !defined $self->{'_elements'}->[0];
1417 return (
1418 $self->{'_elements'}->[0] eq $name
1419 || ( exists $MODEMAP{$name}
1420 && $self->{'_elements'}->[0] eq $MODEMAP{$name} )
1421 );
1422 }
1423
1424 =head2 start_document
1425
1426 Title : start_document
1427 Usage : $eventgenerator->start_document
1428 Function: Handles a start document event
1429 Returns : none
1430 Args : none
1431
1432
1433 =cut
1434
1435 sub start_document {
1436 my ($self) = @_;
1437 $self->{'_lasttype'} = '';
1438 $self->{'_values'} = {};
1439 $self->{'_result'} = undef;
1440 $self->{'_mode'} = '';
1441 $self->{'_elements'} = [];
1442 }
1443
1444 =head2 end_document
1445
1446 Title : end_document
1447 Usage : $eventgenerator->end_document
1448 Function: Handles an end document event
1449 Returns : Bio::Search::Result::ResultI object
1450 Args : none
1451
1452
1453 =cut
1454
1455 sub end_document {
1456 my ( $self, @args ) = @_;
1457 return $self->{'_result'};
1458 }
1459
1460 =head2 idlength
1461
1462 Title : idlength
1463 Usage : $obj->idlength($newval)
1464 Function: Internal storage of the length of the ID desc
1465 in the HSP alignment blocks. Defaults to
1466 $IDLENGTH class variable value
1467 Returns : value of idlength
1468 Args : newvalue (optional)
1469
1470
1471 =cut
1472
1473 sub idlength {
1474 my ( $self, $value ) = @_;
1475 if ( defined $value ) {
1476 $self->{'_idlength'} = $value;
1477 }
1478 return $self->{'_idlength'} || $IDLENGTH;
1479 }
1480
1481 =head2 result_count
1482
1483 Title : result_count
1484 Usage : my $count = $searchio->result_count
1485 Function: Returns the number of results we have processed
1486 Returns : integer
1487 Args : none
1488
1489 =cut
1490
1491 sub result_count {
1492 my $self = shift;
1493 return $self->{'_result_count'};
1494 }
1495
1496 sub attach_EventHandler {
1497 my ( $self, $handler ) = @_;
1498
1499 $self->SUPER::attach_EventHandler($handler);
1500
1501 # Optimization: caching the EventHandler since it is used a lot
1502 # during the parse.
1503
1504 $self->{'_handler_cache'} = $handler;
1505 return;
1506 }
1507
1508 =head2 _will_handle
1509
1510 Title : _will_handle
1511 Usage : Private method. For internal use only.
1512 if( $self->_will_handle($type) ) { ... }
1513 Function: Provides an optimized way to check whether or not an element of a
1514 given type is to be handled.
1515 Returns : Reference to EventHandler object if the element type is to be handled.
1516 undef if the element type is not to be handled.
1517 Args : string containing type of element.
1518
1519 Optimizations:
1520
1521 =over 2
1522
1523 =item 1
1524
1525 Using the cached pointer to the EventHandler to minimize repeated
1526 lookups.
1527
1528 =item 2
1529
1530 Caching the will_handle status for each type that is encountered so
1531 that it only need be checked by calling
1532 handler-E<gt>will_handle($type) once.
1533
1534 =back
1535
1536 This does not lead to a major savings by itself (only 5-10%). In
1537 combination with other optimizations, or for large parse jobs, the
1538 savings good be significant.
1539
1540 To test against the unoptimized version, remove the parentheses from
1541 around the third term in the ternary " ? : " operator and add two
1542 calls to $self-E<gt>_eventHandler().
1543
1544 =cut
1545
1546 sub _will_handle {
1547 my ( $self, $type ) = @_;
1548 my $handler = $self->{'_handler_cache'};
1549 my $will_handle =
1550 defined( $self->{'_will_handle_cache'}->{$type} )
1551 ? $self->{'_will_handle_cache'}->{$type}
1552 : ( $self->{'_will_handle_cache'}->{$type} =
1553 $handler->will_handle($type) );
1554
1555 return $will_handle ? $handler : undef;
1556 }
1557
1558 =head2 _processHits
1559
1560 Title : _processHits
1561 Usage : Private method. For internal use only.
1562 Function: Process/report any hits/hsps we saw in the top table, not in alignments.
1563 Returns : nothing.
1564 Args : array of hits to process.
1565
1566 =cut
1567
1568 sub _processHits {
1569
1570 my ($self, @hit_signifs) = @_;
1571
1572 # process remaining best hits
1573 for my $hit (@hit_signifs) {
1574
1575 # Hsp_score Hsp_evalue Hsp_bit-score
1576 # Hsp_sw-score Hsp_gaps Hsp_identity Hsp_positive
1577 # Hsp_align-len Hsp_query-from Hsp_query-to
1578 # Hsp_hit-from Hsp_hit-to Hsp_qseq Hsp_midline
1579
1580 $self->start_element( { 'Name' => 'Hit' } );
1581 $self->element(
1582 {
1583 'Name' => 'Hit_len',
1584 'Data' => $hit->{hit_len}
1585 }
1586 ) if exists $hit->{hit_len};
1587 $self->element(
1588 {
1589 'Name' => 'Hit_id',
1590 'Data' => $hit->{id}
1591 }
1592 ) if exists $hit->{id};
1593 $self->element(
1594 {
1595 'Name' => 'Hit_accession',
1596 'Data' => $hit->{acc}
1597 }
1598 ) if exists $hit->{acc};
1599 $self->element(
1600 {
1601 'Name' => 'Hit_def',
1602 'Data' => $hit->{desc}
1603 }
1604 ) if exists $hit->{desc};
1605 $self->element(
1606 {
1607 'Name' => 'Hit_signif',
1608 'Data' => $hit->{evalue}
1609 }
1610 ) if exists $hit->{evalue};
1611 $self->element(
1612 {
1613 'Name' => 'Hit_score',
1614 'Data' => $hit->{bits}
1615 }
1616 ) if exists $hit->{bits};
1617
1618 for my $hsp (@{$hit->{HSPs}}) {
1619
1620 $self->start_element( { 'Name' => 'Hsp' } );
1621 $self->element({'Name' => 'Hsp_score', 'Data' => $hsp->{'z-sc'}})
1622 if exists $hsp->{'z-sc'};
1623 $self->element({'Name' => 'Hsp_evalue', 'Data' => $hsp->{evalue} } )
1624 if exists $hsp->{evalue};
1625 $self->element({'Name' => 'Hsp_evalue2', 'Data' => $hsp->{evalue2} } )
1626 if exists $hsp->{evalue2};
1627
1628 $self->element({'Name' => 'Hsp_bit-score', 'Data' => $hsp->{bits} } )
1629 if exists $hsp->{bits};
1630 $self->element({'Name' => 'Hsp_sw-score', 'Data' => $hsp->{'n-w'} } )
1631 if exists $hsp->{'n-w'};
1632 $self->element({'Name' => 'Hsp_sw-score', 'Data' => $hsp->{sw} } )
1633 if exists $hsp->{sw};
1634 $self->element({'Name' => 'Hsp_gaps', 'Data' => $hsp->{'%_gid'} } )
1635 if exists $hsp->{'%_gid'};
1636 $self->element({
1637 'Name' => 'Hsp_identity',
1638 'Data' =>
1639 sprintf( "%.0f", $hsp->{'%_id'} * $hsp->{alen} )
1640 }) if ( exists $hsp->{'%_id'} && exists $hsp->{alen} );
1641
1642 if ( exists $hsp->{'%_gid'} ) {
1643 $self->element(
1644 {
1645 'Name' => 'Hsp_positive',
1646 'Data' =>
1647 sprintf( "%.0f", $hsp->{'%_gid'} * $hsp->{alen} )
1648 }
1649 ) if exists $hsp->{'%_gid'} && exists $hsp->{alen};
1650 } else {
1651 $self->element(
1652 {
1653 'Name' => 'Hsp_positive',
1654 'Data' =>
1655 sprintf( "%.0f", $hsp->{'%_id'} * $hsp->{alen} )
1656 }
1657 ) if ( exists $hsp->{'%_id'} && exists $hsp->{alen} );
1658 }
1659
1660 $self->element(
1661 {
1662 'Name' => 'Hsp_align-len', 'Data' => $hsp->{alen} } )
1663 if exists $hsp->{alen};
1664 $self->element(
1665 {
1666 'Name' => 'Hsp_query-from', 'Data' => $hsp->{an0} } )
1667 if exists $hsp->{an0};
1668 $self->element(
1669 {
1670 'Name' => 'Hsp_query-to', 'Data' => $hsp->{ax0} } )
1671 if exists $hsp->{ax0};
1672 $self->element(
1673 {
1674 'Name' => 'Hsp_hit-from', 'Data' => $hsp->{an1} } )
1675 if exists $hsp->{an1};
1676 $self->element(
1677 {
1678 'Name' => 'Hsp_hit-to', 'Data' => $hsp->{ax1} } )
1679 if exists $hsp->{ax1};
1680
1681 $self->element(
1682 {
1683 'Name' => 'Hsp_querygaps', 'Data' => $hsp->{qgaps} } )
1684 if exists $hsp->{qgaps};
1685 $self->element(
1686 {
1687 'Name' => 'Hsp_hitgaps', 'Data' => $hsp->{lgaps} } )
1688 if exists $hsp->{lgaps};
1689
1690 if ( $self->{'_reporttype'} =~ m/^FAST[NXY]$/o ) {
1691 if ( 8 == scalar grep { exists $hsp->{$_} }
1692 qw(an0 ax0 pn0 px0 an1 ax1 pn1 px1) ) {
1693 if ( $hsp->{ax0} < $hsp->{an0} ) {
1694 $self->element(
1695 {
1696 'Name' => 'Hsp_query-frame',
1697 'Data' =>
1698 "-@{[(($hsp->{px0} - $hsp->{ax0}) % 3) + 1]}"
1699 }
1700 );
1701 } else {
1702 $self->element(
1703 {
1704 'Name' => 'Hsp_query-frame',
1705 'Data' =>
1706 "+@{[(($hsp->{an0} - $hsp->{pn0}) % 3) + 1]}"
1707 }
1708 );
1709 }
1710 if ( $hsp->{ax1} < $hsp->{an1} ) {
1711 $self->element(
1712 {
1713 'Name' => 'Hsp_hit-frame',
1714 'Data' =>
1715 "-@{[(($hsp->{px1} - $hsp->{ax1}) % 3) + 1]}"
1716 }
1717 );
1718 } else {
1719 $self->element(
1720 {
1721 'Name' => 'Hsp_hit-frame',
1722 'Data' =>
1723 "+@{[(($hsp->{an1} - $hsp->{pn1}) % 3) + 1]}"
1724 }
1725 );
1726 }
1727 } else {
1728 $self->element(
1729 {
1730 'Name' => 'Hsp_query-frame',
1731 'Data' => $hsp->{lframe}
1732 }
1733 );
1734 $self->element(
1735 {
1736 'Name' => 'Hsp_hit-frame', 'Data' => 0 } );
1737 }
1738 } else {
1739 $self->element(
1740 {
1741 'Name' => 'Hsp_query-frame', 'Data' => 0 } );
1742 $self->element(
1743 {
1744 'Name' => 'Hsp_hit-frame',
1745 'Data' => $hsp->{lframe}
1746 }
1747 );
1748 }
1749
1750 $self->end_element( { 'Name' => 'Hsp' } );
1751 }
1752
1753 $self->end_element( { 'Name' => 'Hit' } );
1754
1755 }
1756 }
1757
1758
1759
1760 1;
1761
BioPerl core