File Coverage

Bio/SearchIO/blast.pm

Criterion	Covered	Total	%
statement	571	623	91.6
branch	339	390	86.9
condition	122	154	79.2
subroutine	22	24	91.6
pod	12	13	92.3
total	1066	1204	88.5

line	stmt	bran	cond	sub	pod	time	code
1							#
2							# BioPerl module for Bio::SearchIO::blast
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							# 20030409 - sac
15							# PSI-BLAST full parsing support. Rollout of new
16							# model which will remove Steve's old psiblast driver
17							# 20030424 - jason
18							# Megablast parsing fix as reported by Neil Saunders
19							# 20030427 - jason
20							# Support bl2seq parsing
21							# 20031124 - jason
22							# Parse more blast statistics, lambda, entropy, etc
23							# from WU-BLAST in frame-specific manner
24							# 20060216 - cjf - fixed blast parsing for BLAST v2.2.13 output
25							# 20071104 - dmessina - added support for WUBLAST -echofilter
26							# 20071121 - cjf - fixed several bugs (bugs 2391, 2399, 2409)
27
28							=head1 NAME
29
30							Bio::SearchIO::blast - Event generator for event based parsing of
31							blast reports
32
33							=head1 SYNOPSIS
34
35							# Do not use this object directly - it is used as part of the
36							# Bio::SearchIO system.
37
38							use Bio::SearchIO;
39							my $searchio = Bio::SearchIO->new(-format => 'blast',
40							-file => 't/data/ecolitst.bls');
41							while( my $result = $searchio->next_result ) {
42							while( my $hit = $result->next_hit ) {
43							while( my $hsp = $hit->next_hsp ) {
44							# ...
45							}
46							}
47							}
48
49							=head1 DESCRIPTION
50
51							This object encapsulated the necessary methods for generating events
52							suitable for building Bio::Search objects from a BLAST report file.
53							Read the L for more information about how to use this.
54
55							This driver can parse:
56
57							=over 4
58
59							=item *
60
61							NCBI produced plain text BLAST reports from blastall, this also
62							includes PSIBLAST, PSITBLASTN, RPSBLAST, and bl2seq reports. NCBI XML
63							BLAST output is parsed with the blastxml SearchIO driver
64
65							=item *
66
67							WU-BLAST all reports
68
69							=item *
70
71							Jim Kent's BLAST-like output from his programs (BLASTZ, BLAT)
72
73							=item *
74
75							BLAST-like output from Paracel BTK output
76
77							=back
78
79							=head2 bl2seq parsing
80
81							Since I cannot differentiate between BLASTX and TBLASTN since bl2seq
82							doesn't report the algorithm used - I assume it is BLASTX by default -
83							you can supply the program type with -report_type in the SearchIO
84							constructor i.e.
85
86							my $parser = Bio::SearchIO->new(-format => 'blast',
87							-file => 'bl2seq.tblastn.report',
88							-report_type => 'tblastn');
89
90							This only really affects where the frame and strand information are
91							put - they will always be on the $hsp-Equery instead of on the
92							$hsp-Ehit part of the feature pair for blastx and tblastn bl2seq
93							produced reports. Hope that's clear...
94
95							=head1 FEEDBACK
96
97							=head2 Mailing Lists
98
99							User feedback is an integral part of the evolution of this and other
100							Bioperl modules. Send your comments and suggestions preferably to
101							the Bioperl mailing list. Your participation is much appreciated.
102
103							bioperl-l@bioperl.org - General discussion
104							http://bioperl.org/wiki/Mailing_lists - About the mailing lists
105
106							=head2 Support
107
108							Please direct usage questions or support issues to the mailing list:
109
110							I
111
112							rather than to the module maintainer directly. Many experienced and
113							reponsive experts will be able look at the problem and quickly
114							address it. Please include a thorough description of the problem
115							with code and data examples if at all possible.
116
117							=head2 Reporting Bugs
118
119							Report bugs to the Bioperl bug tracking system to help us keep track
120							of the bugs and their resolution. Bug reports can be submitted via the
121							web:
122
123							https://github.com/bioperl/bioperl-live/issues
124
125							=head1 AUTHOR - Jason Stajich
126
127							Email Jason Stajich jason-at-bioperl.org
128
129							=head1 CONTRIBUTORS
130
131							Steve Chervitz sac-at-bioperl.org
132
133							=head1 APPENDIX
134
135							The rest of the documentation details each of the object methods.
136							Internal methods are usually preceded with a _
137
138							=cut
139
140							# Let the code begin...'
141
142							package Bio::SearchIO::blast;
143
144	12			12		4302	use Bio::SearchIO::IteratedSearchResultEventBuilder;
	12					30
	12					344
145	12			12		72	use strict;
	12					18
	12					313
146	12					979	use vars qw(%MAPPING %MODEMAP
147							$DEFAULT_BLAST_WRITER_CLASS
148							$MAX_HSP_OVERLAP
149							$DEFAULT_SIGNIF
150							$DEFAULT_SCORE
151							$DEFAULTREPORTTYPE
152	12			12		58	);
	12					22
153
154
155	12			12		62	use base qw(Bio::SearchIO);
	12					35
	12					769
156	12			12		71	use Data::Dumper;
	12					19
	12					6094
157
158							BEGIN {
159
160							# mapping of NCBI Blast terms to Bioperl hash keys
161	12			12		78	%MODEMAP = (
162							'BlastOutput' => 'result',
163							'Iteration' => 'iteration',
164							'Hit' => 'hit',
165							'Hsp' => 'hsp'
166							);
167
168							# This should really be done more intelligently, like with
169							# XSLT
170
171	12					609	%MAPPING = (
172							'Hsp_bit-score' => 'HSP-bits',
173							'Hsp_score' => 'HSP-score',
174							'Hsp_evalue' => 'HSP-evalue',
175							'Hsp_n', => 'HSP-n',
176							'Hsp_pvalue' => 'HSP-pvalue',
177							'Hsp_query-from' => 'HSP-query_start',
178							'Hsp_query-to' => 'HSP-query_end',
179							'Hsp_hit-from' => 'HSP-hit_start',
180							'Hsp_hit-to' => 'HSP-hit_end',
181							'Hsp_positive' => 'HSP-conserved',
182							'Hsp_identity' => 'HSP-identical',
183							'Hsp_gaps' => 'HSP-hsp_gaps',
184							'Hsp_hitgaps' => 'HSP-hit_gaps',
185							'Hsp_querygaps' => 'HSP-query_gaps',
186							'Hsp_qseq' => 'HSP-query_seq',
187							'Hsp_hseq' => 'HSP-hit_seq',
188							'Hsp_midline' => 'HSP-homology_seq',
189							'Hsp_align-len' => 'HSP-hsp_length',
190							'Hsp_query-frame' => 'HSP-query_frame',
191							'Hsp_hit-frame' => 'HSP-hit_frame',
192							'Hsp_links' => 'HSP-links',
193							'Hsp_group' => 'HSP-hsp_group',
194							'Hsp_features' => 'HSP-hit_features',
195
196							'Hit_id' => 'HIT-name',
197							'Hit_len' => 'HIT-length',
198							'Hit_accession' => 'HIT-accession',
199							'Hit_def' => 'HIT-description',
200							'Hit_signif' => 'HIT-significance',
201							# For NCBI blast, the description line contains bits.
202							# For WU-blast, the description line contains score.
203							'Hit_score' => 'HIT-score',
204							'Hit_bits' => 'HIT-bits',
205
206							'Iteration_iter-num' => 'ITERATION-number',
207							'Iteration_converged' => 'ITERATION-converged',
208
209							'BlastOutput_program' => 'RESULT-algorithm_name',
210							'BlastOutput_version' => 'RESULT-algorithm_version',
211							'BlastOutput_algorithm-reference' => 'RESULT-algorithm_reference',
212							'BlastOutput_rid' => 'RESULT-rid',
213							'BlastOutput_query-def' => 'RESULT-query_name',
214							'BlastOutput_query-len' => 'RESULT-query_length',
215							'BlastOutput_query-acc' => 'RESULT-query_accession',
216							'BlastOutput_query-gi' => 'RESULT-query_gi',
217							'BlastOutput_querydesc' => 'RESULT-query_description',
218							'BlastOutput_db' => 'RESULT-database_name',
219							'BlastOutput_db-len' => 'RESULT-database_entries',
220							'BlastOutput_db-let' => 'RESULT-database_letters',
221							'BlastOutput_inclusion-threshold' => 'RESULT-inclusion_threshold',
222
223							'Parameters_matrix' => { 'RESULT-parameters' => 'matrix' },
224							'Parameters_expect' => { 'RESULT-parameters' => 'expect' },
225							'Parameters_include' => { 'RESULT-parameters' => 'include' },
226							'Parameters_sc-match' => { 'RESULT-parameters' => 'match' },
227							'Parameters_sc-mismatch' => { 'RESULT-parameters' => 'mismatch' },
228							'Parameters_gap-open' => { 'RESULT-parameters' => 'gapopen' },
229							'Parameters_gap-extend' => { 'RESULT-parameters' => 'gapext' },
230							'Parameters_filter' => { 'RESULT-parameters' => 'filter' },
231							'Parameters_allowgaps' => { 'RESULT-parameters' => 'allowgaps' },
232							'Parameters_full_dbpath' => { 'RESULT-parameters' => 'full_dbpath' },
233							'Statistics_db-len' => { 'RESULT-statistics' => 'dbentries' },
234							'Statistics_db-let' => { 'RESULT-statistics' => 'dbletters' },
235							'Statistics_hsp-len' =>
236							{ 'RESULT-statistics' => 'effective_hsplength' },
237							'Statistics_query-len' => { 'RESULT-statistics' => 'querylength' },
238							'Statistics_eff-space' => { 'RESULT-statistics' => 'effectivespace' },
239							'Statistics_eff-spaceused' =>
240							{ 'RESULT-statistics' => 'effectivespaceused' },
241							'Statistics_eff-dblen' =>
242							{ 'RESULT-statistics' => 'effectivedblength' },
243							'Statistics_kappa' => { 'RESULT-statistics' => 'kappa' },
244							'Statistics_lambda' => { 'RESULT-statistics' => 'lambda' },
245							'Statistics_entropy' => { 'RESULT-statistics' => 'entropy' },
246							'Statistics_gapped_kappa' => { 'RESULT-statistics' => 'kappa_gapped' },
247							'Statistics_gapped_lambda' =>
248							{ 'RESULT-statistics' => 'lambda_gapped' },
249							'Statistics_gapped_entropy' =>
250							{ 'RESULT-statistics' => 'entropy_gapped' },
251
252							'Statistics_framewindow' =>
253							{ 'RESULT-statistics' => 'frameshiftwindow' },
254							'Statistics_decay' => { 'RESULT-statistics' => 'decayconst' },
255
256							'Statistics_hit_to_db' => { 'RESULT-statistics' => 'Hits_to_DB' },
257							'Statistics_num_suc_extensions' =>
258							{ 'RESULT-statistics' => 'num_successful_extensions' },
259							'Statistics_length_adjustment' => { 'RESULT-statistics' => 'length_adjustment' },
260							'Statistics_number_of_hsps_better_than_expect_value_cutoff_without_gapping' =>
261							{ 'RESULT-statistics' => 'number_of_hsps_better_than_expect_value_cutoff_without_gapping' },
262							'Statistics_number_of_hsps_gapped' => { 'RESULT-statistics' => 'number_of_hsps_gapped' },
263							'Statistics_number_of_hsps_successfully_gapped' => { 'RESULT-statistics' => 'number_of_hsps_successfully_gapped' },
264
265							# WU-BLAST stats
266							'Statistics_DFA_states' => { 'RESULT-statistics' => 'num_dfa_states' },
267							'Statistics_DFA_size' => { 'RESULT-statistics' => 'dfa_size' },
268							'Statistics_noprocessors' =>
269							{ 'RESULT-statistics' => 'no_of_processors' },
270							'Statistics_neighbortime' =>
271							{ 'RESULT-statistics' => 'neighborhood_generate_time' },
272							'Statistics_starttime' => { 'RESULT-statistics' => 'start_time' },
273							'Statistics_endtime' => { 'RESULT-statistics' => 'end_time' },
274							);
275
276							# add WU-BLAST Frame-Based Statistics
277	12					52	for my $frame ( 0 .. 3 ) {
278	48					63	for my $strand ( '+', '-' ) {
279	96					124	for my $ind (
280							qw(length efflength E S W T X X_gapped E2
281							E2_gapped S2)
282							)
283							{
284	1056					2916	$MAPPING{"Statistics_frame$strand$frame\_$ind"} =
285							{ 'RESULT-statistics' => "Frame$strand$frame\_$ind" };
286							}
287	96					127	for my $val (qw(lambda kappa entropy )) {
288	288					319	for my $type (qw(used computed gapped)) {
289	864					1233	my $key = "Statistics_frame$strand$frame\_$val\_$type";
290	864					1787	my $val =
291							{ 'RESULT-statistics' =>
292							"Frame$strand$frame\_$val\_$type" };
293	864					1619	$MAPPING{$key} = $val;
294							}
295							}
296							}
297							}
298
299							# add Statistics
300	12					29	for my $stats (
301							qw(T A X1 X2 X3 S1 S2 X1_bits X2_bits X3_bits
302							S1_bits S2_bits num_extensions
303							num_successful_extensions
304							seqs_better_than_cutoff
305							posted_date
306							search_cputime total_cputime
307							search_actualtime total_actualtime
308							no_of_processors ctxfactor)
309							)
310							{
311	264					359	my $key = "Statistics_$stats";
312	264					375	my $val = { 'RESULT-statistics' => $stats };
313	264					542	$MAPPING{$key} = $val;
314							}
315
316							# add WU-BLAST Parameters
317	12					22	for my $param (
318							qw(span span1 span2 links warnings notes hspsepsmax
319							hspsepqmax topcomboN topcomboE postsw cpus wordmask
320							filter sort_by_pvalue sort_by_count sort_by_highscore
321							sort_by_totalscore sort_by_subjectlength noseqs gi qtype
322							qres V B Z Y M N)
323							)
324							{
325	348					416	my $key = "Parameters_$param";
326	348					477	my $val = { 'RESULT-parameters' => $param };
327	348					567	$MAPPING{$key} = $val;
328							}
329
330	12					26	$DEFAULT_BLAST_WRITER_CLASS = 'Bio::SearchIO::Writer::HitTableWriter';
331	12					19	$MAX_HSP_OVERLAP = 2; # Used when tiling multiple HSPs.
332	12					42863	$DEFAULTREPORTTYPE = 'BLASTP'; # for bl2seq
333							}
334
335							=head2 new
336
337							Title : new
338							Usage : my $obj = Bio::SearchIO::blast->new(%args);
339							Function: Builds a new Bio::SearchIO::blast object
340							Returns : Bio::SearchIO::blast
341							Args : Key-value pairs:
342							-fh/-file => filehandle/filename to BLAST file
343							-format => 'blast'
344							-report_type => 'blastx', 'tblastn', etc -- only for bl2seq
345							reports when you want to distinguish between
346							tblastn and blastx reports (this only controls
347							where the frame information is put - on the query
348							or subject object.
349							-inclusion_threshold => e-value threshold for inclusion in the
350							PSI-BLAST score matrix model (blastpgp)
351							-signif => float or scientific notation number to be used
352							as a P- or Expect value cutoff
353							-score => integer or scientific notation number to be used
354							as a blast score value cutoff
355							-bits => integer or scientific notation number to be used
356							as a bit score value cutoff
357							-hit_filter => reference to a function to be used for
358							filtering hits based on arbitrary criteria.
359							All hits of each BLAST report must satisfy
360							this criteria to be retained.
361							If a hit fails this test, it is ignored.
362							This function should take a
363							Bio::Search::Hit::BlastHit.pm object as its first
364							argument and return true
365							if the hit should be retained.
366							Sample filter function:
367							-hit_filter => sub { $hit = shift;
368							$hit->gaps == 0; },
369							(Note: -filt_func is synonymous with -hit_filter)
370							-overlap => integer. The amount of overlap to permit between
371							adjacent HSPs when tiling HSPs. A reasonable value is 2.
372							Default = $Bio::SearchIO::blast::MAX_HSP_OVERLAP.
373
374							The following criteria are not yet supported:
375							(these are probably best applied within this module rather than in the
376							event handler since they would permit the parser to take some shortcuts.)
377
378							-check_all_hits => boolean. Check all hits for significance against
379							significance criteria. Default = false.
380							If false, stops processing hits after the first
381							non-significant hit or the first hit that fails
382							the hit_filter call. This speeds parsing,
383							taking advantage of the fact that the hits
384							are processed in the order they appear in the report.
385							-min_query_len => integer to be used as a minimum for query sequence length.
386							Reports with query sequences below this length will
387							not be processed. Default = no minimum length.
388							-best => boolean. Only process the best hit of each report;
389							default = false.
390
391							=cut
392
393							sub _initialize {
394	92			92		280	my ( $self, @args ) = @_;
395	92					409	$self->SUPER::_initialize(@args);
396
397							# Blast reports require a specialized version of the SREB due to the
398							# possibility of iterations (PSI-BLAST). Forwarding all arguments to it. An
399							# issue here is that we want to set new default object factories if none are
400							# supplied.
401
402	92					662	my $handler = Bio::SearchIO::IteratedSearchResultEventBuilder->new(@args);
403	92					262	$self->attach_EventHandler($handler);
404
405							# 2006-04-26 move this to the attach_handler function in this module so we
406							# can really reset the handler
407							# Optimization: caching
408							# the EventHandler since it is used a lot during the parse.
409
410							# $self->{'_handler_cache'} = $handler;
411
412	92					326	my ($rpttype ) = $self->_rearrange(
413							[
414							qw(
415							REPORT_TYPE)
416							],
417							@args
418							);
419	92	100				359	defined $rpttype && ( $self->{'_reporttype'} = $rpttype );
420							}
421
422							sub attach_EventHandler {
423	184			184	1	345	my ($self,$handler) = @_;
424
425	184					656	$self->SUPER::attach_EventHandler($handler);
426
427							# Optimization: caching the EventHandler since it is used a lot
428							# during the parse.
429
430	184					399	$self->{'_handler_cache'} = $handler;
431	184					247	return;
432							}
433
434							=head2 next_result
435
436							Title : next_result
437							Usage : my $hit = $searchio->next_result;
438							Function: Returns the next Result from a search
439							Returns : Bio::Search::Result::ResultI object
440							Args : none
441
442							=cut
443
444							sub next_result {
445	102			102	1	2601	my ($self) = @_;
446	102					278	my $v = $self->verbose;
447	102					179	my $data = '';
448	102					171	my $flavor = '';
449	102					210	$self->{'_seentop'} = 0; # start next report at top
450
451	102					301	my ( $reporttype, $seenquery, $reportline, $reportversion );
452	102					0	my ( $seeniteration, $found_again );
453	102					392	my $incl_threshold = $self->inclusion_threshold;
454	102					148	my $bl2seq_fix;
455	102					378	$self->start_document(); # let the fun begin...
456	102					126	my (@hit_signifs);
457	102					144	my $gapped_stats = 0; # for switching between gapped/ungapped
458							# lambda, K, H
459	102					187	local $_ = "\n"; #consistency
460							PARSER:
461	102					367	while ( defined( $_ = $self->_readline ) ) {
462	20670	100				49471	next if (/^\s+$/); # skip empty lines
463	15331	100				26126	next if (/CPU time:/);
464	15329	50				20019	next if (/^>\s*$/);
465	15329	100				19590	next if (/[]+\s+No hits found\s+[]+/);
466	15325	100	66			164319	if (
		100	100
		100	66
		100	100
		100	66
		100	100
		100	100
		100	100
		100	100
		100	100
		100	100
		50	100
		50	100
		100	100
		100	100
		100	100
		100	66
		100	100
		100	66
		100	66
		100	66
		100
467							/^((?:\S+?)?BLAST[NPX]?)\s+(.+)$/i # NCBI BLAST, PSIBLAST
468							# RPSBLAST, MEGABLAST
469							\|\| /^(P?GENEWISE\|HFRAME\|SWN\|TSWN)\s+(.+)/i #Paracel BTK
470							)
471							{
472	85					384	($reporttype, $reportversion) = ($1, $2);
473							# need to keep track of whether this is WU-BLAST
474	85	50	33			369	if ($reportversion && $reportversion =~ m{WashU$}) {
475	0					0	$self->{'_wublast'}++;
476							}
477	85					613	$self->debug("blast.pm: Start of new report: $reporttype, $reportversion\n");
478	85	100				237	if ( $self->{'_seentop'} ) {
479							# This handles multi-result input streams
480	5					17	$self->_pushback($_);
481	5					11	last PARSER;
482							}
483	80					282	$self->_start_blastoutput;
484	80	100				196	if ($reporttype =~ /RPS-BLAST/) {
485	4					9	$reporttype .= '(BLASTP)'; # default RPS-BLAST type
486							}
487	80					127	$reportline = $_; # to fix the fact that RPS-BLAST output is wrong
488	80					406	$self->element(
489							{
490							'Name' => 'BlastOutput_program',
491							'Data' => $reporttype
492							}
493							);
494
495	80					322	$self->element(
496							{
497							'Name' => 'BlastOutput_version',
498							'Data' => $reportversion
499							}
500							);
501	80					296	$self->element(
502							{
503							'Name' => 'BlastOutput_inclusion-threshold',
504							'Data' => $incl_threshold
505							}
506							);
507							}
508							# parse the BLAST algorithm reference
509							elsif(/^Reference:\s+(.*)$/) {
510							# want to preserve newlines for the BLAST algorithm reference
511	73					262	my $algorithm_reference = "$1\n";
512	73					181	$_ = $self->_readline;
513							# while the current line, does not match an empty line, a RID:, a
514							# Database:, or a query definition line (Query=) we are still
515							# looking at the algorithm_reference, append it to what we parsed so
516							# far
517	73		100			651	while($_ !~ /^$/ && $_ !~ /^RID:/ && $_ !~ /^Database:/ && $_ !~ /^Query=/) {
			100
			66
518	146					399	$algorithm_reference .= "$_";
519	146					249	$_ = $self->_readline;
520							}
521							# if we exited the while loop, we saw an empty line, a RID:, or a
522							# Database:, so push it back
523	73					290	$self->_pushback($_);
524	73					285	$self->element(
525							{
526							'Name' => 'BlastOutput_algorithm-reference',
527							'Data' => $algorithm_reference
528							}
529							);
530							}
531							# parse BLAST RID (Request ID)
532							elsif(/^RID:\s+(.*)$/) {
533	9					28	my $rid = $1;
534	9					42	$self->element(
535							{
536							'Name' => 'BlastOutput_rid',
537							'Data' => $rid
538							}
539							);
540							}
541							# added Windows workaround for bug 1985
542							elsif (/^(Searching\|Results from round)/) {
543	67	100				331	next unless $1 =~ /Results from round/;
544	2					9	$self->debug("blast.pm: Possible psi blast iterations found...\n");
545
546	2	100				6	$self->in_element('hsp')
547							&& $self->end_element( { 'Name' => 'Hsp' } );
548	2	100				7	$self->in_element('hit')
549							&& $self->end_element( { 'Name' => 'Hit' } );
550	2	100				7	if ( defined $seeniteration ) {
551	1	50				2	$self->within_element('iteration')
552							&& $self->end_element( { 'Name' => 'Iteration' } );
553	1					5	$self->start_element( { 'Name' => 'Iteration' } );
554							}
555							else {
556	1					4	$self->start_element( { 'Name' => 'Iteration' } );
557							}
558	2					9	$seeniteration = 1;
559							}
560							elsif (/^Query=\s(.)$/) {
561	102					265	my $q = $1;
562	102					489	$self->debug("blast.pm: Query= found...$_\n");
563	102					180	my $size = 0;
564	102	100				273	if ( defined $seenquery ) {
565	8					31	$self->_pushback($_);
566	8	100				27	$self->_pushback($reportline) if $reportline;
567	8					16	last PARSER;
568							}
569	94	100				233	if ( !defined $reporttype ) {
570	14					52	$self->_start_blastoutput;
571	14	50				450	if ( defined $seeniteration ) {
572	0	0				0	$self->in_element('iteration')
573							&& $self->end_element( { 'Name' => 'Iteration' } );
574	0					0	$self->start_element( { 'Name' => 'Iteration' } );
575							}
576							else {
577	14					52	$self->start_element( { 'Name' => 'Iteration' } );
578							}
579	14					35	$seeniteration = 1;
580							}
581	94					159	$seenquery = $q;
582	94					235	$_ = $self->_readline;
583	94					250	while ( defined($_) ) {
584	142	50				321	if (/^Database:/) {
585	0					0	$self->_pushback($_);
586	0					0	last;
587							}
588							# below line fixes length issue with BLAST v2.2.13; still works
589							# with BLAST v2.2.12
590	142	100	100			818	if ( /$(\-?[\d,]+)\s+letters.*$/ \|\| /^Length=(\-?[\d,]+)/ ) {
591	94					217	$size = $1;
592	94					181	$size =~ s/,//g;
593	94					185	last;
594							}
595							else {
596							# bug 2391
597	48	100	100			446	$q .= ($q =~ /\w$/ && $_ =~ /^\w/) ? " $_" : $_;
598	48					382	$q =~ s/\s+/ /g; # this catches the newline as well
599	48					508	$q =~ s/^ \| $//g;
600							}
601
602	48					129	$_ = $self->_readline;
603							}
604	94					197	chomp($q);
605	94					411	my ( $nm, $desc ) = split( /\s+/, $q, 2 );
606	94	100				549	$self->element(
607							{
608							'Name' => 'BlastOutput_query-def',
609							'Data' => $nm
610							}
611							) if $nm;
612	94					431	$self->element(
613							{
614							'Name' => 'BlastOutput_query-len',
615							'Data' => $size
616							}
617							);
618	94	100				416	defined $desc && $desc =~ s/\s+$//;
619	94					404	$self->element(
620							{
621							'Name' => 'BlastOutput_querydesc',
622							'Data' => $desc
623							}
624							);
625	94					407	my ( $gi, $acc, $version ) = $self->_get_seq_identifiers($nm);
626	94	100	66			404	$version = defined($version) && length($version) ? ".$version" : "";
627	94	100				488	$self->element(
628							{
629							'Name' => 'BlastOutput_query-acc',
630							'Data' => "$acc$version"
631							}
632							) if $acc;
633
634							# these elements are dropped with some multiquery reports; add
635							# back here
636							$self->element(
637							{
638							'Name' => 'BlastOutput_db-len',
639							'Data' => $self->{'_blsdb_length'}
640							}
641	94	100				339	) if $self->{'_blsdb_length'};
642							$self->element(
643							{
644							'Name' => 'BlastOutput_db-let',
645							'Data' => $self->{'_blsdb_letters'}
646							}
647	94	100				301	) if $self->{'_blsdb_letters'};
648							$self->element(
649							{
650							'Name' => 'BlastOutput_db',
651							'Data' => $self->{'_blsdb'}
652							}
653	94	100				416	) if $self->{'_blsdb_letters'};
654							}
655							# added check for WU-BLAST -echofilter option (bug 2388)
656							elsif (/^>Unfiltered[+-]1$/) {
657							# skip all of the lines of unfiltered sequence
658	1					4	while($_ !~ /^Database:/) {
659	69					146	$self->debug("Bypassing features line: $_");
660	69					92	$_ = $self->_readline;
661							}
662	1					5	$self->_pushback($_);
663							}
664							elsif (/Sequences producing significant alignments:/) {
665	56					215	$self->debug("blast.pm: Processing NCBI-BLAST descriptions\n");
666	56					110	$flavor = 'ncbi';
667
668							# PSI-BLAST parsing needs to be fixed to specifically look
669							# for old vs new per iteration, as sorting based on duplication
670							# leads to bugs, see bug 1986
671
672							# The next line is not necessarily whitespace in psiblast reports.
673							# Also note that we must look for the end of this section by testing
674							# for a line with a leading >. Blank lines occur with this section
675							# for psiblast.
676	56	100				137	if ( !$self->in_element('iteration') ) {
677	50					205	$self->start_element( { 'Name' => 'Iteration' } );
678							}
679
680							# changed 8/28/2008 to exit hit table if blank line is found after an
681							# appropriate line
682	56					173	my $h_regex;
683							my $seen_block;
684							DESCLINE:
685	56					174	while ( defined( my $descline = $self->_readline() ) ) {
686	1354	100				3011	if ($descline =~ m{^\s*$}) {
687	105	100				355	last DESCLINE if $seen_block;
688	55					157	next DESCLINE;
689							}
690							# any text match is part of block...
691	1249					1193	$seen_block++;
692							# GCG multiline oddness...
693	1249	100				1737	if ($descline =~ /^(\S+)\s+Begin:\s\d+\s+End:\s+\d+/xms) {
694	3					6	my ($id, $nextline) = ($1, $self->_readline);
695	3					9	$nextline =~ s{^!}{};
696	3					7	$descline = "$id $nextline";
697							}
698							# NCBI style hit table (no N)
699	1249	100				204384	if ($descline =~ /(?
		50
700							(\d*\.?(?:[\+\-eE]+)?\d+) # number (float or scientific notation)
701							\s+ # space
702							(\d*\.?(?:[\+\-eE]+)?\d+) # number (float or scientific notation)
703							\s*$/xms) {
704
705	1243					2686	my ( $score, $evalue ) = ($1, $2);
706
707							# Some data clean-up so e-value will appear numeric to perl
708	1243					1663	$evalue =~ s/^e/1e/i;
709
710							# This to handle no-HSP case
711	1243					3212	my @line = split ' ',$descline;
712
713							# we want to throw away the score, evalue
714	1243					1316	pop @line, pop @line;
715
716							# and N if it is present (of course they are not
717							# really in that order, but it doesn't matter
718	1243	50				2245	if ($3) { pop @line }
	0					0
719
720							# add the last 2 entries s.t. we can reconstruct
721							# a minimal Hit object at the end of the day
722	1243					5443	push @hit_signifs, [ $evalue, $score, shift @line, join( ' ', @line ) ];
723							} elsif ($descline =~ /^CONVERGED/i) {
724	0					0	$self->element(
725							{
726							'Name' => 'Iteration_converged',
727							'Data' => 1
728							}
729							);
730							} else {
731	6					20	$self->_pushback($descline); # Catch leading > (end of section)
732	6					16	last DESCLINE;
733							}
734							}
735							}
736							elsif (/Sequences producing High-scoring Segment Pairs:/) {
737
738							# This block is for WU-BLAST, so we don't have to check for psi-blast stuff
739							# skip the next line
740	26					98	$self->debug("blast.pm: Processing WU-BLAST descriptions\n");
741	26					68	$_ = $self->_readline();
742	26					44	$flavor = 'wu';
743
744	26	100				57	if ( !$self->in_element('iteration') ) {
745	23					90	$self->start_element( { 'Name' => 'Iteration' } );
746							}
747
748	26		66			107	while ( defined( $_ = $self->_readline() )
749							&& !/^\s+$/ )
750							{
751	923					2308	my @line = split;
752	923					900	pop @line; # throw away first number which is for 'N'col
753
754							# add the last 2 entries to array s.t. we can reconstruct
755							# a minimal Hit object at the end of the day
756	923					3469	push @hit_signifs,
757							[ pop @line, pop @line, shift @line, join( ' ', @line ) ];
758							}
759
760							}
761							elsif (/^Database:\s(.+?)\s$/) {
762
763	77					472	$self->debug("blast.pm: Database: $1\n");
764	77					198	my $db = $1;
765	77					206	while ( defined( $_ = $self->_readline ) ) {
766	155	100				546	if (
767							/^\s+(\-?[\d\,]+\|\S+)\s+sequences\;
768							\s+(\-?[\d,]+\|\S+)\s+ # Deal with NCBI 2.2.8 OSX problems
769							total\s+letters/ox
770							)
771							{
772	75					288	my ( $s, $l ) = ( $1, $2 );
773	75					260	$s =~ s/,//g;
774	75					235	$l =~ s/,//g;
775	75					340	$self->element(
776							{
777							'Name' => 'BlastOutput_db-len',
778							'Data' => $s
779							}
780							);
781	75					324	$self->element(
782							{
783							'Name' => 'BlastOutput_db-let',
784							'Data' => $l
785							}
786							);
787							# cache for next round in cases with multiple queries
788	75					191	$self->{'_blsdb'} = $db;
789	75					193	$self->{'_blsdb_length'} = $s;
790	75					171	$self->{'_blsdb_letters'} = $l;
791	75					139	last;
792							}
793							else {
794	80					89	chomp;
795	80					182	$db .= $_;
796							}
797							}
798							$self->element(
799							{
800	77					296	'Name' => 'BlastOutput_db',
801							'Data' => $db
802							}
803							);
804							}
805
806							# move inside of a hit
807							elsif (/^(?:Subject=\|>)\s(\S+)\s(.*)?/) {
808	922					2003	chomp;
809
810	922					4001	$self->debug("blast.pm: Hit: $1\n");
811	922	100				1866	$self->in_element('hsp')
812							&& $self->end_element( { 'Name' => 'Hsp' } );
813	922	100				2192	$self->in_element('hit')
814							&& $self->end_element( { 'Name' => 'Hit' } );
815
816							# special case when bl2seq reports don't have a leading
817							# Query=
818	922	100				2395	if ( !$self->within_element('result') ) {
		100
819	2					12	$self->_start_blastoutput;
820	2					10	$self->start_element( { 'Name' => 'Iteration' } );
821							}
822							elsif ( !$self->within_element('iteration') ) {
823	2					9	$self->start_element( { 'Name' => 'Iteration' } );
824							}
825	922					2735	$self->start_element( { 'Name' => 'Hit' } );
826	922					3210	my $id = $1;
827	922					1517	my $restofline = $2;
828
829	922					4372	$self->debug("Starting a hit: $1 $2\n");
830	922					3470	$self->element(
831							{
832							'Name' => 'Hit_id',
833							'Data' => $id
834							}
835							);
836	922					3089	my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
837	922					2857	$self->element(
838							{
839							'Name' => 'Hit_accession',
840							'Data' => $acc
841							}
842							);
843							# add hit significance (from the hit table)
844							# this is where Bug 1986 went awry
845
846							# Changed for Bug2409; hit->significance and hit->score/bits derived
847							# from HSPs, not hit table unless necessary
848
849							HITTABLE:
850	922					2416	while (my $v = shift @hit_signifs) {
851	642					990	my $tableid = $v->[2];
852	642	100				7542	if ($tableid !~ m{\Q$id\E}) {
853	14					46	$self->debug("Hit table ID $tableid doesn't match current hit id $id, checking next hit table entry...\n");
854	14					39	next HITTABLE;
855							} else {
856	628					1441	last HITTABLE;
857							}
858							}
859	922					2896	while ( defined( $_ = $self->_readline() ) ) {
860	1461	100				3777	next if (/^\s+$/);
861	1458					1711	chomp;
862	1458	100				4003	if (/Length\s=\s([\d,]+)/) {
863	922					1747	my $l = $1;
864	922					1291	$l =~ s/\,//g;
865	922					2842	$self->element(
866							{
867							'Name' => 'Hit_len',
868							'Data' => $l
869							}
870							);
871	922					1787	last;
872							}
873							else {
874	536	100				1255	if ($restofline !~ /\s$/) { # bug #3235
875	519					875	s/^\s(?!\s)/\x01/; #new line to concatenate desc lines with
876							}
877	536	100	100			2414	$restofline .= ($restofline =~ /\w$/ && $_ =~ /^\w/) ? " $_" : $_;
878	536					4457	$restofline =~ s/\s+/ /g; # this catches the newline as well
879	536					6774	$restofline =~ s/^ \| $//g;
880							}
881							}
882	922					5309	$restofline =~ s/\s+/ /g;
883	922					2553	$self->element(
884							{
885							'Name' => 'Hit_def',
886							'Data' => $restofline
887							}
888							);
889							}
890							elsif (/\s+(Plus\|Minus) Strand HSPs:/i) {
891	19					52	next;
892							}
893							elsif (
894							( $self->in_element('hit') \|\| $self->in_element('hsp') )
895							&& # paracel genewise BTK
896							m/Score\s=\s(\S+)\sbits\s # Bit score
897							(?:$(\d+)$)?, # Raw score
898							\s+Log\-Length\sScore\s=\s(\d+) # Log-Length score
899							/ox
900							)
901							{
902	0	0				0	$self->in_element('hsp')
903							&& $self->end_element( { 'Name' => 'Hsp' } );
904	0					0	$self->start_element( { 'Name' => 'Hsp' } );
905
906	0					0	$self->debug( "Got paracel genewise HSP score=$1\n");
907
908							# Some data clean-up so e-value will appear numeric to perl
909	0					0	my ( $bits, $score, $evalue ) = ( $1, $2, $3 );
910	0					0	$evalue =~ s/^e/1e/i;
911	0					0	$self->element(
912							{
913							'Name' => 'Hsp_score',
914							'Data' => $score
915							}
916							);
917	0					0	$self->element(
918							{
919							'Name' => 'Hsp_bit-score',
920							'Data' => $bits
921							}
922							);
923	0					0	$self->element(
924							{
925							'Name' => 'Hsp_evalue',
926							'Data' => $evalue
927							}
928							);
929							}
930							elsif (
931							( $self->in_element('hit') \|\| $self->in_element('hsp') )
932							&& # paracel hframe BTK
933							m/Score\s=\s([^,\s]+), # Raw score
934							\sExpect\s=\s*([^,\s]+), # E-value
935							\sP(?:$\S+$)?\s=\s*([^,\s]+) # P-value
936							/ox
937							)
938							{
939	0	0				0	$self->in_element('hsp')
940							&& $self->end_element( { 'Name' => 'Hsp' } );
941	0					0	$self->start_element( { 'Name' => 'Hsp' } );
942
943	0					0	$self->debug( "Got paracel hframe HSP score=$1\n");
944
945							# Some data clean-up so e-value will appear numeric to perl
946	0					0	my ( $score, $evalue, $pvalue ) = ( $1, $2, $3 );
947	0	0				0	$evalue = "1$evalue" if $evalue =~ /^e/;
948	0	0				0	$pvalue = "1$pvalue" if $pvalue =~ /^e/;
949
950	0					0	$self->element(
951							{
952							'Name' => 'Hsp_score',
953							'Data' => $score
954							}
955							);
956	0					0	$self->element(
957							{
958							'Name' => 'Hsp_evalue',
959							'Data' => $evalue
960							}
961							);
962	0					0	$self->element(
963							{
964							'Name' => 'Hsp_pvalue',
965							'Data' => $pvalue
966							}
967							);
968							}
969							elsif (
970							( $self->in_element('hit') \|\| $self->in_element('hsp') )
971							&& # wublast
972							m/Score\s=\s(\S+)\s* # Bit score
973							$([\d\.]+)\s*bits$, # Raw score
974							\sExpect\s=\s*([^,\s]+), # E-value
975							\s(?:Sum)?\s # SUM
976							P(?:$\d+$)?\s=\s([^,\s]+) # P-value
977							(?:\s,\s+Group\s\=\s*(\d+))? # HSP Group
978							/ox
979							)
980							{ # wu-blast HSP parse
981	423	100				790	$self->in_element('hsp')
982							&& $self->end_element( { 'Name' => 'Hsp' } );
983	423					1440	$self->start_element( { 'Name' => 'Hsp' } );
984
985							# Some data clean-up so e-value will appear numeric to perl
986	423					2284	my ( $score, $bits, $evalue, $pvalue, $group ) =
987							( $1, $2, $3, $4, $5 );
988	423					764	$evalue =~ s/^e/1e/i;
989	423					533	$pvalue =~ s/^e/1e/i;
990
991	423					1319	$self->element(
992							{
993							'Name' => 'Hsp_score',
994							'Data' => $score
995							}
996							);
997	423					1282	$self->element(
998							{
999							'Name' => 'Hsp_bit-score',
1000							'Data' => $bits
1001							}
1002							);
1003	423					1206	$self->element(
1004							{
1005							'Name' => 'Hsp_evalue',
1006							'Data' => $evalue
1007							}
1008							);
1009	423					1171	$self->element(
1010							{
1011							'Name' => 'Hsp_pvalue',
1012							'Data' => $pvalue
1013							}
1014							);
1015
1016	423	100				1539	if ( defined $group ) {
1017	29					74	$self->element(
1018							{
1019							'Name' => 'Hsp_group',
1020							'Data' => $group
1021							}
1022							);
1023							}
1024
1025							}
1026							elsif (
1027							( $self->in_element('hit') \|\| $self->in_element('hsp') )
1028							&& # ncbi blast, works with 2.2.17
1029							m/^\sFeatures\s\w+\sthis\spart/xmso
1030							# If the line begins with "Features in/flanking this part of subject sequence:"
1031							)
1032							{
1033	70	100				139	$self->in_element('hsp')
1034							&& $self->end_element( { 'Name' => 'Hsp' } );
1035	70					125	my $featline;
1036	70					192	$_ = $self->_readline;
1037	70					239	while($_ !~ /^\s*$/) {
1038	119					167	chomp;
1039	119					188	$featline .= $_;
1040	119					187	$_ = $self->_readline;
1041							}
1042	70					195	$self->_pushback($_);
1043	70					728	$featline =~ s{(?:^\s+\|\s+^)}{}g;
1044	70					121	$featline =~ s{\n}{;}g;
1045	70					241	$self->start_element( { 'Name' => 'Hsp' } );
1046	70					343	$self->element(
1047							{
1048							'Name' => 'Hsp_features',
1049							'Data' => $featline
1050							}
1051							);
1052	70					254	$self->{'_seen_hsp_features'} = 1;
1053							}
1054							elsif (
1055							( $self->in_element('hit') \|\| $self->in_element('hsp') )
1056							&& # ncbi blast, works with 2.2.17
1057							m/Score\s=\s(\S+)\sbits\s # Bit score
1058							(?:$(\d+)$)?, # Missing for BLAT pseudo-BLAST fmt
1059							\sExpect(?:$(\d+\+?)$)?\s=\s*([^,\s]+) # E-value
1060							/ox
1061							)
1062							{ # parse NCBI blast HSP
1063	1200	100				2356	if( !$self->{'_seen_hsp_features'} ) {
1064	1130	100				1661	$self->in_element('hsp')
1065							&& $self->end_element( { 'Name' => 'Hsp' } );
1066	1130					3236	$self->start_element( { 'Name' => 'Hsp' } );
1067							}
1068
1069							# Some data clean-up so e-value will appear numeric to perl
1070	1200					5699	my ( $bits, $score, $n, $evalue ) = ( $1, $2, $3, $4 );
1071	1200					2478	$evalue =~ s/^e/1e/i;
1072	1200					4112	$self->element(
1073							{
1074							'Name' => 'Hsp_score',
1075							'Data' => $score
1076							}
1077							);
1078	1200					3648	$self->element(
1079							{
1080							'Name' => 'Hsp_bit-score',
1081							'Data' => $bits
1082							}
1083							);
1084	1200					3632	$self->element(
1085							{
1086							'Name' => 'Hsp_evalue',
1087							'Data' => $evalue
1088							}
1089							);
1090	1200	100				2376	$self->element(
1091							{
1092							'Name' => 'Hsp_n',
1093							'Data' => $n
1094							}
1095							) if defined $n;
1096	1200	50				1743	$score = '' unless defined $score; # deal with BLAT which
1097							# has no score only bits
1098	1200					4193	$self->debug("Got NCBI HSP score=$score, evalue $evalue\n");
1099							}
1100							elsif (
1101							$self->in_element('hsp')
1102							&& m/Identities\s=\s(\d+)\s\/\s(\d+)\s[\d\%]+\s
1103							(?:,\sPositives\s=\s(\d+)\/(\d+)\s[\d\%]+\s*)? # pos only valid for Protein alignments
1104							(?:\,\sGaps\s=\s*(\d+)\/(\d+))? # Gaps
1105							/oxi
1106							)
1107							{
1108	1623					6704	$self->element(
1109							{
1110							'Name' => 'Hsp_identity',
1111							'Data' => $1
1112							}
1113							);
1114	1623					5681	$self->element(
1115							{
1116							'Name' => 'Hsp_align-len',
1117							'Data' => $2
1118							}
1119							);
1120	1623	100				3624	if ( defined $3 ) {
1121	1251					2974	$self->element(
1122							{
1123							'Name' => 'Hsp_positive',
1124							'Data' => $3
1125							}
1126							);
1127							}
1128							else {
1129	372					940	$self->element(
1130							{
1131							'Name' => 'Hsp_positive',
1132							'Data' => $1
1133							}
1134							);
1135							}
1136	1623	100				3652	if ( defined $6 ) {
1137	787					1928	$self->element(
1138							{
1139							'Name' => 'Hsp_gaps',
1140							'Data' => $5
1141							}
1142							);
1143							}
1144
1145	1623					4576	$self->{'_Query'} = { 'begin' => 0, 'end' => 0 };
1146	1623					3329	$self->{'_Sbjct'} = { 'begin' => 0, 'end' => 0 };
1147
1148	1623	100				6214	if (/(Frame\s=\s.+)$/) {
1149
1150							# handle wu-blast Frame listing on same line
1151	44					117	$self->_pushback($1);
1152							}
1153							}
1154							elsif ( $self->in_element('hsp')
1155							&& /Strand\s=\s(Plus\|Minus)\s\/\s(Plus\|Minus)/i )
1156							{
1157
1158							# consume this event ( we infer strand from start/end)
1159	372	100				772	if (!defined($reporttype)) {
1160	2					5	$self->{'_reporttype'} = $reporttype = 'BLASTN';
1161	2					3	$bl2seq_fix = 1; # special case to resubmit the algorithm
1162							# reporttype
1163							}
1164	372					1176	next;
1165							}
1166							elsif ( $self->in_element('hsp')
1167							&& /Links\s=\s(\S+)/ox )
1168							{
1169	342					1132	$self->element(
1170							{
1171							'Name' => 'Hsp_links',
1172							'Data' => $1
1173							}
1174							);
1175							}
1176							elsif ( $self->in_element('hsp')
1177							&& /Frame\s=\s([\+\-][1-3])\s(\/\s([\+\-][1-3]))?/ )
1178							{
1179	308		50			897	my $frame1 = $1 \|\| 0;
1180	308		100			775	my $frame2 = $2 \|\| 0;
1181							# this is for bl2seq only
1182	308	100				547	if ( not defined $reporttype ) {
1183	4					8	$bl2seq_fix = 1;
1184	4	100	66			24	if ( $frame1 && $frame2 ) {
1185	1					2	$reporttype = 'TBLASTX'
1186							}
1187							else {
1188							# We can distinguish between BLASTX and TBLASTN from the report
1189							# (and assign $frame1 properly) by using the start/end from query.
1190							# If the report is BLASTX, the coordinates distance from query
1191							# will be 3 times the length of the alignment shown (coordinates in nt,
1192							# alignment in aa); if not then subject is the nucleotide sequence (TBLASTN).
1193							# Will have to fast-forward to query alignment line and then go back.
1194	3					30	my $fh = $self->_fh;
1195	3					13	my $file_pos = tell $fh;
1196
1197	3					7	my $a_position = '';
1198	3					7	my $ali_length = '';
1199	3					7	my $b_position = '';
1200	3					14	while (my $line = <$fh>) {
1201	6	100				35	if ($line =~ m/^(?:Query\|Sbjct):?\s+(\-?\d+)?\s*(\S+)\s+(\-?\d+)?/) {
1202	3					14	$a_position = $1;
1203	3					10	my $alignment = $2;
1204	3					8	$b_position = $3;
1205
1206	12			12		4448	use Bio::LocatableSeq;
	12					37
	12					58998
1207	3					9	my $gap_symbols = $Bio::LocatableSeq::GAP_SYMBOLS;
1208	3					61	$alignment =~ s/[$gap_symbols]//g;
1209	3					11	$ali_length = length($alignment);
1210	3					9	last;
1211							}
1212							}
1213	3	100				22	my $coord_length = ($a_position < $b_position) ? ($b_position - $a_position + 1)
1214							: ($a_position - $b_position + 1);
1215	3	100				14	($coord_length == ($ali_length * 3)) ? ($reporttype = 'BLASTX') : ($reporttype = 'TBLASTN');
1216
1217							# Rewind filehandle to its original position to continue parsing
1218	3					26	seek $fh, $file_pos, 0;
1219							}
1220	4					13	$self->{'_reporttype'} = $reporttype;
1221							}
1222
1223	308					387	my ( $queryframe, $hitframe );
1224	308	100	66			720	if ( $reporttype eq 'TBLASTX' ) {
		100	33
		50
1225	187					297	( $queryframe, $hitframe ) = ( $frame1, $frame2 );
1226	187					625	$hitframe =~ s/\/\s*//g;
1227							}
1228							elsif ( $reporttype eq 'TBLASTN' \|\| $reporttype eq 'PSITBLASTN') {
1229	95					179	( $hitframe, $queryframe ) = ( $frame1, 0 );
1230							}
1231							elsif ( $reporttype eq 'BLASTX' \|\| $reporttype eq 'RPS-BLAST(BLASTP)') {
1232	26					48	( $queryframe, $hitframe ) = ( $frame1, 0 );
1233							# though NCBI doesn't report it, this is a special BLASTX-like
1234							# RPS-BLAST; should be handled differently
1235	26	50				41	if ($reporttype eq 'RPS-BLAST(BLASTP)') {
1236	0					0	$self->element(
1237							{
1238							'Name' => 'BlastOutput_program',
1239							'Data' => 'RPS-BLAST(BLASTX)'
1240							}
1241							);
1242							}
1243							}
1244							$self->element(
1245							{
1246	308					1033	'Name' => 'Hsp_query-frame',
1247							'Data' => $queryframe
1248							}
1249							);
1250
1251	308					896	$self->element(
1252							{
1253							'Name' => 'Hsp_hit-frame',
1254							'Data' => $hitframe
1255							}
1256							);
1257							}
1258							elsif (/^Parameters:/
1259							\|\| /^\s+Database:\s+?/
1260							\|\| /^\s+Subset/
1261							\|\| /^\s*Lambda/
1262							\|\| /^\s*Histogram/
1263							\|\| ( $self->in_element('hsp') && /WARNING\|NOTE/ ) )
1264							{
1265
1266							# Note: Lambda check was necessary to parse
1267							# t/data/ecoli_domains.rpsblast AND to parse bl2seq
1268	83					310	$self->debug("blast.pm: found parameters section \n");
1269
1270	83	100				226	$self->in_element('hsp')
1271							&& $self->end_element( { 'Name' => 'Hsp' } );
1272	83	100				300	$self->in_element('hit')
1273							&& $self->end_element( { 'Name' => 'Hit' } );
1274
1275							# This is for the case when we specify -b 0 (or B=0 for WU-BLAST)
1276							# and still want to construct minimal Hit objects
1277	83	100				294	$self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1278	83	100				289	$self->within_element('iteration')
1279							&& $self->end_element( { 'Name' => 'Iteration' } );
1280
1281	83	50				290	next if /^\s+Subset/;
1282	83	100				400	my $blast = (/^(\s+Database\:)\|(\s*Lambda)/) ? 'ncbi' : 'wublast';
1283	83	50				242	if (/^\s*Histogram/) {
1284	0					0	$blast = 'btk';
1285							}
1286
1287	83					132	my $last = '';
1288
1289							# default is that gaps are allowed
1290	83					525	$self->element(
1291							{
1292							'Name' => 'Parameters_allowgaps',
1293							'Data' => 'yes'
1294							}
1295							);
1296	83					312	while ( defined( $_ = $self->_readline ) ) {
1297							# If Lambda/Kappa/Entropy numbers appear first at this point,
1298							# pushback and add the header line to process it correctly
1299	2819	100	100			16207	if (/^\s+[\d+\.]+\s+[\d+\.]+\s+[\d+\.]/ and $last eq '') {
		100	66
		100
1300	15					66	$self->_pushback($_);
1301	15					46	$self->_pushback("Lambda K H\n");
1302	15					31	next;
1303							}
1304							elsif (
1305							/^((?:\S+)?BLAST[NPX]?)\s+(.+)$/i # NCBI BLAST, PSIBLAST
1306							# RPSBLAST, MEGABLAST
1307							\|\| /^(P?GENEWISE\|HFRAME\|SWN\|TSWN)\s+(.+)/i #Paracel BTK
1308							)
1309							{
1310	5					32	$self->_pushback($_);
1311
1312							# let's handle this in the loop
1313	5					27	last;
1314							}
1315							elsif (/^Query=/) {
1316	4					15	$self->_pushback($_);
1317	4	50				22	$self->_pushback($reportline) if $reportline;
1318	4					12	last PARSER;
1319							}
1320
1321							# here is where difference between wublast and ncbiblast
1322							# is better handled by different logic
1323	2795	100	100			11275	if ( /Number of Sequences:\s+([\d\,]+)/i
		100
		50
		100
		50
1324							\|\| /of sequences in database:\s+(\-?[\d,]+)/i )
1325							{
1326	125					253	my $c = $1;
1327	125					229	$c =~ s/\,//g;
1328	125					451	$self->element(
1329							{
1330							'Name' => 'Statistics_db-len',
1331							'Data' => $c
1332							}
1333							);
1334							}
1335							elsif (/letters in database:\s+(\-?[\d,]+)/i) {
1336	73					188	my $s = $1;
1337	73					272	$s =~ s/,//g;
1338	73					263	$self->element(
1339							{
1340							'Name' => 'Statistics_db-let',
1341							'Data' => $s
1342							}
1343							);
1344							}
1345							elsif ( $blast eq 'btk' ) {
1346	0					0	next;
1347							}
1348							elsif ( $blast eq 'wublast' ) {
1349
1350							# warn($_);
1351	985	100				8245	if (/E=(\S+)/) {
		100
		100
		100
		100
		100
		100
		100
		50
		100
		100
		100
		100
		100
1352	29					139	$self->element(
1353							{
1354							'Name' => 'Parameters_expect',
1355							'Data' => $1
1356							}
1357							);
1358							}
1359							elsif (/nogaps/) {
1360	2					8	$self->element(
1361							{
1362							'Name' => 'Parameters_allowgaps',
1363							'Data' => 'no'
1364							}
1365							);
1366							}
1367							elsif (/ctxfactor=(\S+)/) {
1368	26					121	$self->element(
1369							{
1370							'Name' => 'Statistics_ctxfactor',
1371							'Data' => $1
1372							}
1373							);
1374							}
1375							elsif (
1376							/(postsw\|links\|span[12]?\|warnings\|notes\|gi\|noseqs\|qres\|qype)/
1377							)
1378							{
1379	34					172	$self->element(
1380							{
1381							'Name' => "Parameters_$1",
1382							'Data' => 'yes'
1383							}
1384							);
1385							}
1386							elsif (/(\S+)=(\S+)/) {
1387	141					590	$self->element(
1388							{
1389							'Name' => "Parameters_$1",
1390							'Data' => $2
1391							}
1392							);
1393							}
1394							elsif ( $last =~ /(Frame\|Strand)\s+MatID\s+Matrix name/i ) {
1395	26					42	my $firstgapinfo = 1;
1396	26					47	my $frame = undef;
1397	26		66			140	while ( defined($_) && !/^\s+$/ ) {
1398	108					283	s/^\s+//;
1399	108					420	s/\s+$//;
1400	108	100	100			371	if ( $firstgapinfo
1401							&& s/Q=(\d+),R=(\d+)\s+//x )
1402							{
1403	24					68	$firstgapinfo = 0;
1404
1405	24					125	$self->element(
1406							{
1407							'Name' => 'Parameters_gap-open',
1408							'Data' => $1
1409							}
1410							);
1411	24					120	$self->element(
1412							{
1413							'Name' => 'Parameters_gap-extend',
1414							'Data' => $2
1415							}
1416							);
1417	24					100	my @fields = split;
1418
1419	24					52	for my $type (
1420							qw(lambda_gapped
1421							kappa_gapped
1422							entropy_gapped)
1423							)
1424							{
1425	72	50				130	next if $type eq 'n/a';
1426	72	50				118	if ( !@fields ) {
1427	0					0	warn "fields is empty for $type\n";
1428	0					0	next;
1429							}
1430							$self->element(
1431							{
1432	72					233	'Name' =>
1433							"Statistics_frame$frame\_$type",
1434							'Data' => shift @fields
1435							}
1436							);
1437							}
1438							}
1439							else {
1440	84					350	my ( $frameo, $matid, $matrix, @fields ) =
1441							split;
1442	84	100				173	if ( !defined $frame ) {
1443
1444							# keep some sort of default feature I guess
1445							# even though this is sort of wrong
1446	26					112	$self->element(
1447							{
1448							'Name' => 'Parameters_matrix',
1449							'Data' => $matrix
1450							}
1451							);
1452	26					121	$self->element(
1453							{
1454							'Name' => 'Statistics_lambda',
1455							'Data' => $fields[0]
1456							}
1457							);
1458	26					105	$self->element(
1459							{
1460							'Name' => 'Statistics_kappa',
1461							'Data' => $fields[1]
1462							}
1463							);
1464	26					84	$self->element(
1465							{
1466							'Name' => 'Statistics_entropy',
1467							'Data' => $fields[2]
1468							}
1469							);
1470							}
1471	84					115	$frame = $frameo;
1472	84					92	my $ii = 0;
1473	84					120	for my $type (
1474							qw(lambda_used
1475							kappa_used
1476							entropy_used
1477							lambda_computed
1478							kappa_computed
1479							entropy_computed)
1480							)
1481							{
1482	504					544	my $f = $fields[$ii];
1483	504	100				681	next unless defined $f; # deal with n/a
1484	456	100				601	if ( $f eq 'same' ) {
1485	72					97	$f = $fields[ $ii - 3 ];
1486							}
1487	456					356	$ii++;
1488	456					1152	$self->element(
1489							{
1490							'Name' =>
1491							"Statistics_frame$frame\_$type",
1492							'Data' => $f
1493							}
1494							);
1495
1496							}
1497							}
1498
1499							# get the next line
1500	108					245	$_ = $self->_readline;
1501							}
1502	26					59	$last = $_;
1503							}
1504							elsif ( $last =~ /(Frame\|Strand)\s+MatID\s+Length/i ) {
1505	26					46	my $frame = undef;
1506	26		33			171	while ( defined($_) && !/^\s+/ ) {
1507	0					0	s/^\s+//;
1508	0					0	s/\s+$//;
1509	0					0	my @fields = split;
1510	0	0				0	if ( @fields <= 3 ) {
1511	0					0	for my $type (qw(X_gapped E2_gapped S2)) {
1512	0	0				0	last unless @fields;
1513	0					0	$self->element(
1514							{
1515							'Name' =>
1516							"Statistics_frame$frame\_$type",
1517							'Data' => shift @fields
1518							}
1519							);
1520							}
1521							}
1522							else {
1523
1524	0					0	for my $type (
1525							qw(length
1526							efflength
1527							E S W T X E2 S2)
1528							)
1529							{
1530	0					0	$self->element(
1531							{
1532							'Name' =>
1533							"Statistics_frame$frame\_$type",
1534							'Data' => shift @fields
1535							}
1536							);
1537							}
1538							}
1539	0					0	$_ = $self->_readline;
1540							}
1541	26					49	$last = $_;
1542							}
1543							elsif (/(\S+\s+\S+)\s+DFA:\s+(\S+)\s+$(.+)$/) {
1544	23	50				99	if ( $1 eq 'states in' ) {
		0
1545	23					145	$self->element(
1546							{
1547							'Name' => 'Statistics_DFA_states',
1548							'Data' => "$2 $3"
1549							}
1550							);
1551							}
1552							elsif ( $1 eq 'size of' ) {
1553	0					0	$self->element(
1554							{
1555							'Name' => 'Statistics_DFA_size',
1556							'Data' => "$2 $3"
1557							}
1558							);
1559							}
1560							}
1561							elsif (
1562							m/^\s+Time to generate neighborhood:\s+
1563							(\S+\s+\S+\s+\S+)/x
1564							)
1565							{
1566	0					0	$self->element(
1567							{
1568							'Name' => 'Statistics_neighbortime',
1569							'Data' => $1
1570							}
1571							);
1572							}
1573							elsif (/processors\s+used:\s+(\d+)/) {
1574	23					111	$self->element(
1575							{
1576							'Name' => 'Statistics_noprocessors',
1577							'Data' => $1
1578							}
1579							);
1580							}
1581							elsif (
1582							m/^\s+(\S+)\s+cpu\s+time:\s+ # cputype
1583							(\S+\s+\S+\s+\S+) # cputime
1584							\s+Elapsed:\s+(\S+)/x
1585							)
1586							{
1587	46					126	my $cputype = lc($1);
1588	46					207	$self->element(
1589							{
1590							'Name' => "Statistics_$cputype\_cputime",
1591							'Data' => $2
1592							}
1593							);
1594	46					190	$self->element(
1595							{
1596							'Name' => "Statistics_$cputype\_actualtime",
1597							'Data' => $3
1598							}
1599							);
1600							}
1601							elsif (/^\s+Start:/) {
1602	23					221	my ( $junk, $start, $stime, $end, $etime ) =
1603							split( /\s+(Start\|End)\:\s+/, $_ );
1604	23					53	chomp($stime);
1605	23					99	$self->element(
1606							{
1607							'Name' => 'Statistics_starttime',
1608							'Data' => $stime
1609							}
1610							);
1611	23					51	chomp($etime);
1612	23					70	$self->element(
1613							{
1614							'Name' => 'Statistics_endtime',
1615							'Data' => $etime
1616							}
1617							);
1618							}
1619							elsif (/^\s+Database:\s+(.+)$/) {
1620	21					125	$self->element(
1621							{
1622							'Name' => 'Parameters_full_dbpath',
1623							'Data' => $1
1624							}
1625							);
1626
1627							}
1628							elsif (/^\s+Posted:\s+(.+)/) {
1629	21					61	my $d = $1;
1630	21					35	chomp($d);
1631	21					89	$self->element(
1632							{
1633							'Name' => 'Statistics_posted_date',
1634							'Data' => $d
1635							}
1636							);
1637							}
1638							}
1639							elsif ( $blast eq 'ncbi' ) {
1640
1641	1612	100				11880	if (m/^Matrix:\s+(.+)\s*$/oxi) {
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
		100
1642	55					281	$self->element(
1643							{
1644							'Name' => 'Parameters_matrix',
1645							'Data' => $1
1646							}
1647							);
1648							}
1649							elsif (/^Gapped/) {
1650	49					99	$gapped_stats = 1;
1651							}
1652							elsif (/^Lambda/) {
1653	107					263	$_ = $self->_readline;
1654	107					340	s/^\s+//;
1655	107					377	my ( $lambda, $kappa, $entropy ) = split;
1656	107	100				243	if ($gapped_stats) {
1657	49					215	$self->element(
1658							{
1659							'Name' => "Statistics_gapped_lambda",
1660							'Data' => $lambda
1661							}
1662							);
1663	49					190	$self->element(
1664							{
1665							'Name' => "Statistics_gapped_kappa",
1666							'Data' => $kappa
1667							}
1668							);
1669	49					141	$self->element(
1670							{
1671							'Name' => "Statistics_gapped_entropy",
1672							'Data' => $entropy
1673							}
1674							);
1675							}
1676							else {
1677	58					235	$self->element(
1678							{
1679							'Name' => "Statistics_lambda",
1680							'Data' => $lambda
1681							}
1682							);
1683	58					272	$self->element(
1684							{
1685							'Name' => "Statistics_kappa",
1686							'Data' => $kappa
1687							}
1688							);
1689	58					237	$self->element(
1690							{
1691							'Name' => "Statistics_entropy",
1692							'Data' => $entropy
1693							}
1694							);
1695							}
1696							}
1697							elsif (m/effective\s+search\s+space\s+used:\s+(\d+)/oxi) {
1698	54					262	$self->element(
1699							{
1700							'Name' => 'Statistics_eff-spaceused',
1701							'Data' => $1
1702							}
1703							);
1704							}
1705							elsif (m/effective\s+search\s+space:\s+(\d+)/oxi) {
1706	48					249	$self->element(
1707							{
1708							'Name' => 'Statistics_eff-space',
1709							'Data' => $1
1710							}
1711							);
1712							}
1713							elsif (
1714							m/Gap\s+Penalties:\s+Existence:\s+(\d+)\,
1715							\s+Extension:\s+(\d+)/ox
1716							)
1717							{
1718	48					323	$self->element(
1719							{
1720							'Name' => 'Parameters_gap-open',
1721							'Data' => $1
1722							}
1723							);
1724	48					195	$self->element(
1725							{
1726							'Name' => 'Parameters_gap-extend',
1727							'Data' => $2
1728							}
1729							);
1730							}
1731							elsif (/effective\s+HSP\s+length:\s+(\d+)/) {
1732	45					231	$self->element(
1733							{
1734							'Name' => 'Statistics_hsp-len',
1735							'Data' => $1
1736							}
1737							);
1738							}
1739							elsif (/Number\s+of\s+HSP's\s+better\s+than\s+(\S+)\s+without\s+gapping:\s+(\d+)/) {
1740	34					161	$self->element(
1741							{
1742							'Name' => 'Statistics_number_of_hsps_better_than_expect_value_cutoff_without_gapping',
1743							'Data' => $2
1744							}
1745							);
1746							}
1747							elsif (/Number\s+of\s+HSP's\s+gapped:\s+(\d+)/) {
1748	7					35	$self->element(
1749							{
1750							'Name' => 'Statistics_number_of_hsps_gapped',
1751							'Data' => $1
1752							}
1753							);
1754							}
1755							elsif (/Number\s+of\s+HSP's\s+successfully\s+gapped:\s+(\d+)/) {
1756	7					30	$self->element(
1757							{
1758							'Name' => 'Statistics_number_of_hsps_successfully_gapped',
1759							'Data' => $1
1760							}
1761							);
1762							}
1763							elsif (/Length\s+adjustment:\s+(\d+)/) {
1764	6					33	$self->element(
1765							{
1766							'Name' => 'Statistics_length_adjustment',
1767							'Data' => $1
1768							}
1769							);
1770							}
1771							elsif (/effective\s+length\s+of\s+query:\s+([\d\,]+)/i) {
1772	48					141	my $c = $1;
1773	48					100	$c =~ s/\,//g;
1774	48					174	$self->element(
1775							{
1776							'Name' => 'Statistics_query-len',
1777							'Data' => $c
1778							}
1779							);
1780							}
1781							elsif (/effective\s+length\s+of\s+database:\s+([\d\,]+)/i) {
1782	51					117	my $c = $1;
1783	51					164	$c =~ s/\,//g;
1784	51					216	$self->element(
1785							{
1786							'Name' => 'Statistics_eff-dblen',
1787							'Data' => $c
1788							}
1789							);
1790							}
1791							elsif (
1792							/^(T\|A\|X1\|X2\|X3\|S1\|S2):\s+(\d+(\.\d+)?)\s+(?:$\s*(\d+\.\d+) bits$)?/
1793							)
1794							{
1795	325					608	my $v = $2;
1796	325					376	chomp($v);
1797	325					1241	$self->element(
1798							{
1799							'Name' => "Statistics_$1",
1800							'Data' => $v
1801							}
1802							);
1803	325	100				781	if ( defined $4 ) {
1804	235					731	$self->element(
1805							{
1806							'Name' => "Statistics_$1_bits",
1807							'Data' => $4
1808							}
1809							);
1810							}
1811							}
1812							elsif (
1813							m/frameshift\s+window\,
1814							\s+decay\s+const:\s+(\d+)\,\s+([\.\d]+)/x
1815							)
1816							{
1817	14					83	$self->element(
1818							{
1819							'Name' => 'Statistics_framewindow',
1820							'Data' => $1
1821							}
1822							);
1823	14					67	$self->element(
1824							{
1825							'Name' => 'Statistics_decay',
1826							'Data' => $2
1827							}
1828							);
1829							}
1830							elsif (m/^Number\s+of\s+Hits\s+to\s+DB:\s+(\S+)/ox) {
1831	52					244	$self->element(
1832							{
1833							'Name' => 'Statistics_hit_to_db',
1834							'Data' => $1
1835							}
1836							);
1837							}
1838							elsif (m/^Number\s+of\s+extensions:\s+(\S+)/ox) {
1839	51					222	$self->element(
1840							{
1841							'Name' => 'Statistics_num_extensions',
1842							'Data' => $1
1843							}
1844							);
1845							}
1846							elsif (
1847							m/^Number\s+of\s+successful\s+extensions:\s+
1848							(\S+)/ox
1849							)
1850							{
1851	51					276	$self->element(
1852							{
1853							'Name' => 'Statistics_num_suc_extensions',
1854							'Data' => $1
1855							}
1856							);
1857							}
1858							elsif (
1859							m/^Number\s+of\s+sequences\s+better\s+than\s+
1860							(\S+):\s+(\d+)/ox
1861							)
1862							{
1863	52					269	$self->element(
1864							{
1865							'Name' => 'Parameters_expect',
1866							'Data' => $1
1867							}
1868							);
1869	52					218	$self->element(
1870							{
1871							'Name' => 'Statistics_seqs_better_than_cutoff',
1872							'Data' => $2
1873							}
1874							);
1875							}
1876							elsif (/^\s+Posted\s+date:\s+(.+)/) {
1877	50					181	my $d = $1;
1878	50					101	chomp($d);
1879	50					175	$self->element(
1880							{
1881							'Name' => 'Statistics_posted_date',
1882							'Data' => $d
1883							}
1884							);
1885							}
1886							elsif ( !/^\s+$/ ) {
1887							#$self->debug( "unmatched stat $_");
1888							}
1889							}
1890	2795					6084	$last = $_;
1891							}
1892							} elsif ( $self->in_element('hsp') ) {
1893	9058					20376	$self->debug("blast.pm: Processing HSP\n");
1894							# let's read 3 lines at a time;
1895							# bl2seq hackiness... Not sure I like
1896	9058		66			14461	$self->{'_reporttype'} \|\|= $DEFAULTREPORTTYPE;
1897	9058					19671	my %data = (
1898							'Query' => '',
1899							'Mid' => '',
1900							'Hit' => ''
1901							);
1902	9058					8199	my $len;
1903	9058		66			24252	for ( my $i = 0 ; defined($_) && $i < 3 ; $i++ ) {
1904							# $self->debug("$i: $_") if $v;
1905	27150	50	66			104251	if ( ( $i == 0 && /^\s+$/) \|\|
			33
1906							/^\s*(?:Lambda\|Minus\|Plus\|Score)/i
1907							) {
1908	0	0				0	$self->_pushback($_) if defined $_;
1909	0					0	$self->end_element( { 'Name' => 'Hsp' } );
1910	0					0	last;
1911							}
1912	27150					29340	chomp;
1913	27150	100				77259	if (/^((Query\|Sbjct):?\s+(\-?\d+)?\s*)(\S+)\s+(\-?\d+)?/) {
		50
1914	18116					53822	my ( $full, $type, $start, $str, $end ) =
1915							( $1, $2, $3, $4, $5 );
1916
1917	18116	100				22174	if ( $str eq '-' ) {
1918	48	100				94	$i = 3 if $type eq 'Sbjct';
1919							}
1920							else {
1921	18068					25773	$data{$type} = $str;
1922							}
1923	18116					18091	$len = length($full);
1924							$self->{"\_$type"}->{'begin'} = $start
1925	18116	100				41738	unless $self->{"_$type"}->{'begin'};
1926	18116					27796	$self->{"\_$type"}->{'end'} = $end;
1927							} elsif (/^((Query\|Sbjct):?\s+(\-?0+)\s*)/) {
1928							# Bug from NCBI's BLAST: unaligned output
1929	0					0	$_ = $self->_readline() for 0..1;
1930	0					0	last;
1931							} else {
1932	9034	50	33			21382	$self->throw("no data for midline $_")
1933							unless ( defined $_ && defined $len );
1934	9034					15521	$data{'Mid'} = substr( $_, $len );
1935							}
1936	27150					43915	$_ = $self->_readline();
1937							}
1938							$self->characters(
1939							{
1940							'Name' => 'Hsp_qseq',
1941	9058					29692	'Data' => $data{'Query'}
1942							}
1943							);
1944							$self->characters(
1945							{
1946							'Name' => 'Hsp_hseq',
1947	9058					25373	'Data' => $data{'Sbjct'}
1948							}
1949							);
1950							$self->characters(
1951							{
1952							'Name' => 'Hsp_midline',
1953	9058					20854	'Data' => $data{'Mid'}
1954							}
1955							);
1956							}
1957							else {
1958							#$self->debug("blast.pm: unrecognized line $_");
1959							}
1960							}
1961
1962	102					442	$self->debug("blast.pm: End of BlastOutput\n");
1963	102	100				313	if ( $self->{'_seentop'} ) {
1964	96	100				253	$self->within_element('hsp')
1965							&& $self->end_element( { 'Name' => 'Hsp' } );
1966	96	100				223	$self->within_element('hit')
1967							&& $self->end_element( { 'Name' => 'Hit' } );
1968							# cleanup extra hits
1969	96	100				238	$self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1970	96	100				238	$self->within_element('iteration')
1971							&& $self->end_element( { 'Name' => 'Iteration' } );
1972	96	100				260	if ($bl2seq_fix) {
1973	6					156	$self->element(
1974							{
1975							'Name' => 'BlastOutput_program',
1976							'Data' => $reporttype
1977							}
1978							);
1979							}
1980	96					331	$self->end_element( { 'Name' => 'BlastOutput' } );
1981							}
1982	102					381	return $self->end_document();
1983							}
1984
1985							# Private method for internal use only.
1986							sub _start_blastoutput {
1987	96			96		167	my $self = shift;
1988	96					476	$self->start_element( { 'Name' => 'BlastOutput' } );
1989	96					187	$self->{'_seentop'} = 1;
1990	96					255	$self->{'_result_count'}++;
1991	96					171	$self->{'_handler_rc'} = undef;
1992							}
1993
1994							=head2 _will_handle
1995
1996							Title : _will_handle
1997							Usage : Private method. For internal use only.
1998							if( $self->_will_handle($type) ) { ... }
1999							Function: Provides an optimized way to check whether or not an element of a
2000							given type is to be handled.
2001							Returns : Reference to EventHandler object if the element type is to be handled.
2002							undef if the element type is not to be handled.
2003							Args : string containing type of element.
2004
2005							Optimizations:
2006
2007							=over 2
2008
2009							=item 1
2010
2011							Using the cached pointer to the EventHandler to minimize repeated
2012							lookups.
2013
2014							=item 2
2015
2016							Caching the will_handle status for each type that is encountered so
2017							that it only need be checked by calling
2018							handler-Ewill_handle($type) once.
2019
2020							=back
2021
2022							This does not lead to a major savings by itself (only 5-10%). In
2023							combination with other optimizations, or for large parse jobs, the
2024							savings good be significant.
2025
2026							To test against the unoptimized version, remove the parentheses from
2027							around the third term in the ternary " ? : " operator and add two
2028							calls to $self-E_eventHandler().
2029
2030							=cut
2031
2032							sub _will_handle {
2033	8516			8516		10504	my ( $self, $type ) = @_;
2034	8516					10278	my $handler = $self->{'_handler_cache'};
2035							my $will_handle =
2036							defined( $self->{'_will_handle_cache'}->{$type} )
2037							? $self->{'_will_handle_cache'}->{$type}
2038	8516	100				17376	: ( $self->{'_will_handle_cache'}->{$type} =
2039							$handler->will_handle($type) );
2040
2041	8516	50				13170	return $will_handle ? $handler : undef;
2042							}
2043
2044							=head2 start_element
2045
2046							Title : start_element
2047							Usage : $eventgenerator->start_element
2048							Function: Handles a start element event
2049							Returns : none
2050							Args : hashref with at least 2 keys 'Data' and 'Name'
2051
2052							=cut
2053
2054							sub start_element {
2055	4258			4258	1	5467	my ( $self, $data ) = @_;
2056
2057							# we currently don't care about attributes
2058	4258					4791	my $nm = $data->{'Name'};
2059	4258					5591	my $type = $MODEMAP{$nm};
2060	4258	50				6491	if ($type) {
2061	4258					6474	my $handler = $self->_will_handle($type);
2062	4258	50				6754	if ($handler) {
2063	4258					10762	my $func = sprintf( "start_%s", lc $type );
2064	4258					15160	$self->{'_handler_rc'} = $handler->$func( $data->{'Attributes'} );
2065							}
2066							#else {
2067							#$self->debug( # changed 4/29/2006 to play nice with other event handlers
2068							# "Bio::SearchIO::InternalParserError ".
2069							# "\nCan't handle elements of type \'$type.\'"
2070							#);
2071							#}
2072	4258					6515	unshift @{ $self->{'_elements'} }, $type;
	4258					7428
2073	4258	100				6431	if ( $type eq 'result' ) {
2074	96					229	$self->{'_values'} = {};
2075	96					221	$self->{'_result'} = undef;
2076							} else {
2077							# cleanup some things
2078	4162	50				7540	if ( defined $self->{'_values'} ) {
2079	4162					3945	foreach my $k (
2080	117000					258080	grep { /^\U$type\-/ }
2081	4162					18281	keys %{ $self->{'_values'} }
2082							)
2083							{
2084	38633					44156	delete $self->{'_values'}->{$k};
2085							}
2086							}
2087							}
2088							}
2089							}
2090
2091							=head2 end_element
2092
2093							Title : end_element
2094							Usage : $eventgenerator->end_element
2095							Function: Handles an end element event
2096							Returns : hashref with an element's worth of data
2097							Args : hashref with at least 2 keys 'Data' and 'Name'
2098
2099
2100							=cut
2101
2102							sub end_element {
2103	43007			43007	1	47030	my ( $self, $data ) = @_;
2104
2105	43007					42463	my $nm = $data->{'Name'};
2106	43007					39075	my $type;
2107							my $rc;
2108							# cache these (TODO: we should probably cache all cross-report data)
2109	43007	100				53546	if ( $nm eq 'BlastOutput_program' ) {
2110	86	100				452	if ( $self->{'_last_data'} =~ /(t?blast[npx])/i ) {
2111	83					268	$self->{'_reporttype'} = uc $1;
2112							}
2113	86		66			202	$self->{'_reporttype'} \|\|= $DEFAULTREPORTTYPE;
2114							}
2115
2116	43007	100				49338	if ( $nm eq 'BlastOutput_version' ) {
2117	80					155	$self->{'_reportversion'} = $self->{'_last_data'};
2118							}
2119
2120							# Hsps are sort of weird, in that they end when another
2121							# object begins so have to detect this in end_element for now
2122	43007	100				50529	if ( $nm eq 'Hsp' ) {
2123	1623					2412	foreach (qw(Hsp_qseq Hsp_midline Hsp_hseq Hsp_features)) {
2124							$self->element(
2125							{
2126							'Name' => $_,
2127							'Data' => $self->{'_last_hspdata'}->{$_}
2128							}
2129	6492	100				17549	) if defined $self->{'_last_hspdata'}->{$_};
2130							}
2131	1623					3431	$self->{'_last_hspdata'} = {};
2132							$self->element(
2133							{
2134							'Name' => 'Hsp_query-from',
2135	1623					4922	'Data' => $self->{'_Query'}->{'begin'}
2136							}
2137							);
2138							$self->element(
2139							{
2140							'Name' => 'Hsp_query-to',
2141	1623					5077	'Data' => $self->{'_Query'}->{'end'}
2142							}
2143							);
2144
2145							$self->element(
2146							{
2147							'Name' => 'Hsp_hit-from',
2148	1623					4995	'Data' => $self->{'_Sbjct'}->{'begin'}
2149							}
2150							);
2151							$self->element(
2152							{
2153							'Name' => 'Hsp_hit-to',
2154	1623					4247	'Data' => $self->{'_Sbjct'}->{'end'}
2155							}
2156							);
2157
2158							# } elsif( $nm eq 'Iteration' ) {
2159							# Nothing special needs to be done here.
2160							}
2161	43007	100				79360	if ( $type = $MODEMAP{$nm} ) {
		100
2162	4258					6883	my $handler = $self->_will_handle($type);
2163	4258	50				6708	if ($handler) {
2164	4258					11156	my $func = sprintf( "end_%s", lc $type );
2165	4258					13889	$rc = $handler->$func( $self->{'_reporttype'}, $self->{'_values'} );
2166							}
2167	4258					4495	shift @{ $self->{'_elements'} };
	4258					6292
2168
2169							}
2170							elsif ( $MAPPING{$nm} ) {
2171	38562	100				52162	if ( ref( $MAPPING{$nm} ) =~ /hash/i ) {
2172
2173							# this is where we shove in the data from the
2174							# hashref info about params or statistics
2175	2900					2933	my $key = ( keys %{ $MAPPING{$nm} } )[0];
	2900					6243
2176							$self->{'_values'}->{$key}->{ $MAPPING{$nm}->{$key} } =
2177	2900					8168	$self->{'_last_data'};
2178							}
2179							else {
2180	35662					61898	$self->{'_values'}->{ $MAPPING{$nm} } = $self->{'_last_data'};
2181							}
2182							}
2183							else {
2184							#$self->debug("blast.pm: unknown nm $nm, ignoring\n");
2185							}
2186	43007					44495	$self->{'_last_data'} = ''; # remove read data if we are at
2187							# end of an element
2188	43007	100	100			62218	$self->{'_result'} = $rc if ( defined $type && $type eq 'result' );
2189	43007					40358	$self->{'_seen_hsp_features'} = 0;
2190	43007					57287	return $rc;
2191							}
2192
2193							=head2 element
2194
2195							Title : element
2196							Usage : $eventhandler->element({'Name' => $name, 'Data' => $str});
2197							Function: Convenience method that calls start_element, characters, end_element
2198							Returns : none
2199							Args : Hash ref with the keys 'Name' and 'Data'
2200
2201
2202							=cut
2203
2204							sub element {
2205	38749			38749	1	44881	my ( $self, $data ) = @_;
2206							# Note that start element isn't needed for character data
2207							# Not too SAX-y, though
2208							#$self->start_element($data);
2209	38749					54889	$self->characters($data);
2210	38749					51034	$self->end_element($data);
2211							}
2212
2213							=head2 characters
2214
2215							Title : characters
2216							Usage : $eventgenerator->characters($str)
2217							Function: Send a character events
2218							Returns : none
2219							Args : string
2220
2221
2222							=cut
2223
2224							sub characters {
2225	65923			65923	1	71288	my ( $self, $data ) = @_;
2226	65923	100	100			72078	if ( $self->in_element('hsp')
2227							&& $data->{'Name'} =~ /^Hsp\_(qseq\|hseq\|midline)$/ )
2228							{
2229							$self->{'_last_hspdata'}->{ $data->{'Name'} } .= $data->{'Data'}
2230	32019	100				85014	if defined $data->{'Data'};
2231							}
2232	65923	100	100			194531	return unless ( defined $data->{'Data'} && $data->{'Data'} !~ /^\s+$/ );
2233	65853					109782	$self->{'_last_data'} = $data->{'Data'};
2234							}
2235
2236							=head2 within_element
2237
2238							Title : within_element
2239							Usage : if( $eventgenerator->within_element($element) ) {}
2240							Function: Test if we are within a particular element
2241							This is different than 'in' because within can be tested
2242							for a whole block.
2243							Returns : boolean
2244							Args : string element name
2245
2246							See Also: L
2247
2248							=cut
2249
2250							sub within_element {
2251	2214			2214	1	3087	my ( $self, $name ) = @_;
2252							return 0
2253							if ( !defined $name && !defined $self->{'_elements'}
2254	2214	100	33			4369	\|\| scalar @{ $self->{'_elements'} } == 0 );
	2214		66			4685
2255	2212					2168	foreach ( @{ $self->{'_elements'} } ) {
	2212					3091
2256	3134	100				4501	if ( $_ eq $name ) {
2257	1949					4592	return 1;
2258							}
2259							}
2260	263					531	return 0;
2261							}
2262
2263							=head2 in_element
2264
2265							Title : in_element
2266							Usage : if( $eventgenerator->in_element($element) ) {}
2267							Function: Test if we are in a particular element
2268							This is different than 'within_element' because within
2269							can be tested for a whole block.
2270							Returns : boolean
2271							Args : string element name
2272
2273							See Also: L
2274
2275							=cut
2276
2277							sub in_element {
2278	264029			264029	1	305149	my ( $self, $name ) = @_;
2279	264029	100				337271	return 0 if !defined $self->{'_elements'}->[0];
2280	263910					847437	return ( $self->{'_elements'}->[0] eq $name );
2281							}
2282
2283							=head2 start_document
2284
2285							Title : start_document
2286							Usage : $eventgenerator->start_document
2287							Function: Handle a start document event
2288							Returns : none
2289							Args : none
2290
2291
2292							=cut
2293
2294							sub start_document {
2295	102			102	1	198	my ($self) = @_;
2296	102					203	$self->{'_lasttype'} = '';
2297	102					424	$self->{'_values'} = {};
2298	102					217	$self->{'_result'} = undef;
2299	102					221	$self->{'_elements'} = [];
2300							}
2301
2302							=head2 end_document
2303
2304							Title : end_document
2305							Usage : $eventgenerator->end_document
2306							Function: Handles an end document event
2307							Returns : Bio::Search::Result::ResultI object
2308							Args : none
2309
2310
2311							=cut
2312
2313							sub end_document {
2314	102			102	1	225	my ( $self, @args ) = @_;
2315
2316							#$self->debug("blast.pm: end_document\n");
2317	102					745	return $self->{'_result'};
2318							}
2319
2320							sub write_result {
2321	5			5	1	37	my ( $self, $blast, @args ) = @_;
2322
2323	5	50				27	if ( not defined( $self->writer ) ) {
2324	0					0	$self->warn("Writer not defined. Using a $DEFAULT_BLAST_WRITER_CLASS");
2325	0					0	$self->writer( $DEFAULT_BLAST_WRITER_CLASS->new() );
2326							}
2327	5					34	$self->SUPER::write_result( $blast, @args );
2328							}
2329
2330							sub result_count {
2331	0			0	1	0	my $self = shift;
2332	0					0	return $self->{'_result_count'};
2333							}
2334
2335	0			0	0	0	sub report_count { shift->result_count }
2336
2337							=head2 inclusion_threshold
2338
2339							Title : inclusion_threshold
2340							Usage : my $incl_thresh = $isreb->inclusion_threshold;
2341							: $isreb->inclusion_threshold(1e-5);
2342							Function: Get/Set the e-value threshold for inclusion in the PSI-BLAST
2343							score matrix model (blastpgp) that was used for generating the reports
2344							being parsed.
2345							Returns : number (real)
2346							Default value: $Bio::SearchIO::IteratedSearchResultEventBuilder::DEFAULT_INCLUSION_THRESHOLD
2347							Args : number (real) (e.g., 0.0001 or 1e-4 )
2348
2349							=cut
2350
2351							=head2 max_significance
2352
2353							Usage : $obj->max_significance();
2354							Purpose : Set/Get the P or Expect value used as significance screening cutoff.
2355							This is the value of the -signif parameter supplied to new().
2356							Hits with P or E-value above this are skipped.
2357							Returns : Scientific notation number with this format: 1.0e-05.
2358							Argument : Scientific notation number or float (when setting)
2359							Comments : Screening of significant hits uses the data provided on the
2360							: description line. For NCBI BLAST1 and WU-BLAST, this data
2361							: is P-value. for NCBI BLAST2 it is an Expect value.
2362
2363							=cut
2364
2365							=head2 signif
2366
2367							Synonym for L
2368
2369							=cut
2370
2371							=head2 min_score
2372
2373							Usage : $obj->min_score();
2374							Purpose : Set/Get the Blast score used as screening cutoff.
2375							This is the value of the -score parameter supplied to new().
2376							Hits with scores below this are skipped.
2377							Returns : Integer or scientific notation number.
2378							Argument : Integer or scientific notation number (when setting)
2379							Comments : Screening of significant hits uses the data provided on the
2380							: description line.
2381
2382							=cut
2383
2384							=head2 min_query_length
2385
2386							Usage : $obj->min_query_length();
2387							Purpose : Gets the query sequence length used as screening criteria.
2388							This is the value of the -min_query_len parameter supplied to new().
2389							Hits with sequence length below this are skipped.
2390							Returns : Integer
2391							Argument : n/a
2392
2393							=cut
2394
2395							=head2 best_hit_only
2396
2397							Title : best_hit_only
2398							Usage : print "only getting best hit.\n" if $obj->best_hit_only;
2399							Purpose : Set/Get the indicator for whether or not to process only
2400							: the best BlastHit.
2401							Returns : Boolean (1 \| 0)
2402							Argument : Boolean (1 \| 0) (when setting)
2403
2404							=cut
2405
2406							=head2 check_all_hits
2407
2408							Title : check_all_hits
2409							Usage : print "checking all hits.\n" if $obj->check_all_hits;
2410							Purpose : Set/Get the indicator for whether or not to process all hits.
2411							: If false, the parser will stop processing hits after the
2412							: the first non-significance hit or the first hit that fails
2413							: any hit filter.
2414							Returns : Boolean (1 \| 0)
2415							Argument : Boolean (1 \| 0) (when setting)
2416
2417							=cut
2418
2419							# general private method used to make minimal hits from leftover
2420							# data in the hit table
2421
2422							sub _cleanup_hits {
2423	8			8		18	my ($self, $hits) = @_;
2424	8					13	while ( my $v = shift @{ $hits }) {
	1532					2871
2425	1524	50				2113	next unless defined $v;
2426	1524					3945	$self->start_element( { 'Name' => 'Hit' } );
2427	1524					4156	my $id = $v->[2];
2428	1524					1818	my $desc = $v->[3];
2429	1524					4370	$self->element(
2430							{
2431							'Name' => 'Hit_id',
2432							'Data' => $id
2433							}
2434							);
2435	1524					4315	my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
2436	1524					3929	$self->element(
2437							{
2438							'Name' => 'Hit_accession',
2439							'Data' => $acc
2440							}
2441							);
2442	1524	50				2892	if ( defined $v ) {
2443	1524					4004	$self->element(
2444							{
2445							'Name' => 'Hit_signif',
2446							'Data' => $v->[0]
2447							}
2448							);
2449	1524	50				2846	if (exists $self->{'_wublast'}) {
2450	0					0	$self->element(
2451							{
2452							'Name' => 'Hit_score',
2453							'Data' => $v->[1]
2454							}
2455							);
2456							} else {
2457	1524					3524	$self->element(
2458							{
2459							'Name' => 'Hit_bits',
2460							'Data' => $v->[1]
2461							}
2462							);
2463							}
2464
2465							}
2466							$self->element(
2467							{
2468	1524					4096	'Name' => 'Hit_def',
2469							'Data' => $desc
2470							}
2471							);
2472	1524					3235	$self->end_element( { 'Name' => 'Hit' } );
2473							}
2474							}
2475
2476
2477							1;
2478
2479							__END__