File Coverage

Bio/SearchIO/blast.pm

Criterion	Covered	Total	%
statement	570	623	91.4
branch	339	390	86.9
condition	120	154	77.9
subroutine	22	24	91.6
pod	12	13	92.3
total	1063	1204	88.2

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		4266	use Bio::SearchIO::IteratedSearchResultEventBuilder;
	12					18
	12					327
145	12			12		51	use strict;
	12					13
	12					255
146	12					849	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		36	);
	12					15
153
154
155	12			12		37	use base qw(Bio::SearchIO);
	12					16
	12					700
156	12			12		44	use Data::Dumper;
	12					12
	12					5042
157
158							BEGIN {
159
160							# mapping of NCBI Blast terms to Bioperl hash keys
161	12			12		57	%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					538	%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					47	for my $frame ( 0 .. 3 ) {
278	48					44	for my $strand ( '+', '-' ) {
279	96					84	for my $ind (
280							qw(length efflength E S W T X X_gapped E2
281							E2_gapped S2)
282							)
283							{
284	1056					2080	$MAPPING{"Statistics_frame$strand$frame\_$ind"} =
285							{ 'RESULT-statistics' => "Frame$strand$frame\_$ind" };
286							}
287	96					81	for my $val (qw(lambda kappa entropy )) {
288	288					211	for my $type (qw(used computed gapped)) {
289	864					802	my $key = "Statistics_frame$strand$frame\_$val\_$type";
290	864					1148	my $val =
291							{ 'RESULT-statistics' =>
292							"Frame$strand$frame\_$val\_$type" };
293	864					1051	$MAPPING{$key} = $val;
294							}
295							}
296							}
297							}
298
299							# add Statistics
300	12					19	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					203	my $key = "Statistics_$stats";
312	264					236	my $val = { 'RESULT-statistics' => $stats };
313	264					366	$MAPPING{$key} = $val;
314							}
315
316							# add WU-BLAST Parameters
317	12					17	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					256	my $key = "Parameters_$param";
326	348					356	my $val = { 'RESULT-parameters' => $param };
327	348					398	$MAPPING{$key} = $val;
328							}
329
330	12					15	$DEFAULT_BLAST_WRITER_CLASS = 'Bio::SearchIO::Writer::HitTableWriter';
331	12					15	$MAX_HSP_OVERLAP = 2; # Used when tiling multiple HSPs.
332	12					35647	$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		240	my ( $self, @args ) = @_;
395	92					352	$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					690	my $handler = Bio::SearchIO::IteratedSearchResultEventBuilder->new(@args);
403	92					221	$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					333	my ($rpttype ) = $self->_rearrange(
413							[
414							qw(
415							REPORT_TYPE)
416							],
417							@args
418							);
419	92	100				406	defined $rpttype && ( $self->{'_reporttype'} = $rpttype );
420							}
421
422							sub attach_EventHandler {
423	184			184	1	209	my ($self,$handler) = @_;
424
425	184					595	$self->SUPER::attach_EventHandler($handler);
426
427							# Optimization: caching the EventHandler since it is used a lot
428							# during the parse.
429
430	184					278	$self->{'_handler_cache'} = $handler;
431	184					282	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	2271	my ($self) = @_;
446	102					289	my $v = $self->verbose;
447	102					186	my $data = '';
448	102					140	my $flavor = '';
449	102					194	$self->{'_seentop'} = 0; # start next report at top
450
451	102					126	my ( $reporttype, $seenquery, $reportline, $reportversion );
452	0					0	my ( $seeniteration, $found_again );
453	102					385	my $incl_threshold = $self->inclusion_threshold;
454	102					131	my $bl2seq_fix;
455	102					291	$self->start_document(); # let the fun begin...
456	102					118	my (@hit_signifs);
457	102					182	my $gapped_stats = 0; # for switching between gapped/ungapped
458							# lambda, K, H
459	102					204	local $_ = "\n"; #consistency
460							PARSER:
461	102					418	while ( defined( $_ = $self->_readline ) ) {
462	20670	100				45889	next if (/^\s+$/); # skip empty lines
463	15331	100				22472	next if (/CPU time:/);
464	15329	50				18904	next if (/^>\s*$/);
465	15329	100				19636	next if (/[]+\s+No hits found\s+[]+/);
466	15325	100	66			186180	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					332	($reporttype, $reportversion) = ($1, $2);
473							# need to keep track of whether this is WU-BLAST
474	85	50	33			449	if ($reportversion && $reportversion =~ m{WashU$}) {
475	0					0	$self->{'_wublast'}++;
476							}
477	85					568	$self->debug("blast.pm: Start of new report: $reporttype, $reportversion\n");
478	85	100				264	if ( $self->{'_seentop'} ) {
479							# This handles multi-result input streams
480	5					15	$self->_pushback($_);
481	5					12	last PARSER;
482							}
483	80					278	$self->_start_blastoutput;
484	80	100				232	if ($reporttype =~ /RPS-BLAST/) {
485	4					10	$reporttype .= '(BLASTP)'; # default RPS-BLAST type
486							}
487	80					119	$reportline = $_; # to fix the fact that RPS-BLAST output is wrong
488	80					397	$self->element(
489							{
490							'Name' => 'BlastOutput_program',
491							'Data' => $reporttype
492							}
493							);
494
495	80					336	$self->element(
496							{
497							'Name' => 'BlastOutput_version',
498							'Data' => $reportversion
499							}
500							);
501	80					311	$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					273	my $algorithm_reference = "$1\n";
512	73					206	$_ = $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		66			609	while($_ !~ /^$/ && $_ !~ /^RID:/ && $_ !~ /^Database:/ && $_ !~ /^Query=/) {
			66
			66
518	146					331	$algorithm_reference .= "$_";
519	146					256	$_ = $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					266	$self->_pushback($_);
524	73					275	$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					25	my $rid = $1;
534	9					39	$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				362	next unless $1 =~ /Results from round/;
544	2					6	$self->debug("blast.pm: Possible psi blast iterations found...\n");
545
546	2	100				7	$self->in_element('hsp')
547							&& $self->end_element( { 'Name' => 'Hsp' } );
548	2	100				6	$self->in_element('hit')
549							&& $self->end_element( { 'Name' => 'Hit' } );
550	2	100				7	if ( defined $seeniteration ) {
551	1	50				4	$self->within_element('iteration')
552							&& $self->end_element( { 'Name' => 'Iteration' } );
553	1					5	$self->start_element( { 'Name' => 'Iteration' } );
554							}
555							else {
556	1					5	$self->start_element( { 'Name' => 'Iteration' } );
557							}
558	2					10	$seeniteration = 1;
559							}
560							elsif (/^Query=\s(.)$/) {
561	102					262	my $q = $1;
562	102					495	$self->debug("blast.pm: Query= found...$_\n");
563	102					165	my $size = 0;
564	102	100				253	if ( defined $seenquery ) {
565	8					24	$self->_pushback($_);
566	8	100				24	$self->_pushback($reportline) if $reportline;
567	8					13	last PARSER;
568							}
569	94	100				215	if ( !defined $reporttype ) {
570	14					50	$self->_start_blastoutput;
571	14	50				36	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					48	$self->start_element( { 'Name' => 'Iteration' } );
578							}
579	14					28	$seeniteration = 1;
580							}
581	94					151	$seenquery = $q;
582	94					223	$_ = $self->_readline;
583	94					270	while ( defined($_) ) {
584	142	50				298	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			866	if ( /$(\-?[\d,]+)\s+letters.*$/ \|\| /^Length=(\-?[\d,]+)/ ) {
591	94					210	$size = $1;
592	94					174	$size =~ s/,//g;
593	94					150	last;
594							}
595							else {
596							# bug 2391
597	48	100	100			371	$q .= ($q =~ /\w$/ && $_ =~ /^\w/) ? " $_" : $_;
598	48					366	$q =~ s/\s+/ /g; # this catches the newline as well
599	48					496	$q =~ s/^ \| $//g;
600							}
601
602	48					133	$_ = $self->_readline;
603							}
604	94					206	chomp($q);
605	94					414	my ( $nm, $desc ) = split( /\s+/, $q, 2 );
606	94	100				557	$self->element(
607							{
608							'Name' => 'BlastOutput_query-def',
609							'Data' => $nm
610							}
611							) if $nm;
612	94					458	$self->element(
613							{
614							'Name' => 'BlastOutput_query-len',
615							'Data' => $size
616							}
617							);
618	94	100				402	defined $desc && $desc =~ s/\s+$//;
619	94					332	$self->element(
620							{
621							'Name' => 'BlastOutput_querydesc',
622							'Data' => $desc
623							}
624							);
625	94					408	my ( $gi, $acc, $version ) = $self->_get_seq_identifiers($nm);
626	94	100	66			461	$version = defined($version) && length($version) ? ".$version" : "";
627	94	100				523	$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				337	) if $self->{'_blsdb_length'};
642							$self->element(
643							{
644							'Name' => 'BlastOutput_db-let',
645							'Data' => $self->{'_blsdb_letters'}
646							}
647	94	100				345	) if $self->{'_blsdb_letters'};
648							$self->element(
649							{
650							'Name' => 'BlastOutput_db',
651							'Data' => $self->{'_blsdb'}
652							}
653	94	100				419	) 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					5	while($_ !~ /^Database:/) {
659	69					113	$self->debug("Bypassing features line: $_");
660	69					82	$_ = $self->_readline;
661							}
662	1					4	$self->_pushback($_);
663							}
664							elsif (/Sequences producing significant alignments:/) {
665	56					210	$self->debug("blast.pm: Processing NCBI-BLAST descriptions\n");
666	56					91	$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				148	if ( !$self->in_element('iteration') ) {
677	50					190	$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					152	my $h_regex;
683							my $seen_block;
684							DESCLINE:
685	56					195	while ( defined( my $descline = $self->_readline() ) ) {
686	1354	100				2524	if ($descline =~ m{^\s*$}) {
687	105	100				283	last DESCLINE if $seen_block;
688	55					184	next DESCLINE;
689							}
690							# any text match is part of block...
691	1249					771	$seen_block++;
692							# GCG multiline oddness...
693	1249	100				1633	if ($descline =~ /^(\S+)\s+Begin:\s\d+\s+End:\s+\d+/xms) {
694	3					7	my ($id, $nextline) = ($1, $self->_readline);
695	3					9	$nextline =~ s{^!}{};
696	3					6	$descline = "$id $nextline";
697							}
698							# NCBI style hit table (no N)
699	1249	100				204025	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					1713	my ( $score, $evalue ) = ($1, $2);
706
707							# Some data clean-up so e-value will appear numeric to perl
708	1243					1225	$evalue =~ s/^e/1e/i;
709
710							# This to handle no-HSP case
711	1243					2833	my @line = split ' ',$descline;
712
713							# we want to throw away the score, evalue
714	1243					922	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				2009	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					4994	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					27	$self->_pushback($descline); # Catch leading > (end of section)
732	6					18	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					84	$self->debug("blast.pm: Processing WU-BLAST descriptions\n");
741	26					61	$_ = $self->_readline();
742	26					40	$flavor = 'wu';
743
744	26	100				58	if ( !$self->in_element('iteration') ) {
745	23					86	$self->start_element( { 'Name' => 'Iteration' } );
746							}
747
748	26		66			99	while ( defined( $_ = $self->_readline() )
749							&& !/^\s+$/ )
750							{
751	923					2818	my @line = split;
752	923					782	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					4689	push @hit_signifs,
757							[ pop @line, pop @line, shift @line, join( ' ', @line ) ];
758							}
759
760							}
761							elsif (/^Database:\s(.+?)\s$/) {
762
763	77					428	$self->debug("blast.pm: Database: $1\n");
764	77					162	my $db = $1;
765	77					208	while ( defined( $_ = $self->_readline ) ) {
766	155	100				582	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					216	my ( $s, $l ) = ( $1, $2 );
773	75					180	$s =~ s/,//g;
774	75					203	$l =~ s/,//g;
775	75					318	$self->element(
776							{
777							'Name' => 'BlastOutput_db-len',
778							'Data' => $s
779							}
780							);
781	75					297	$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					201	$self->{'_blsdb'} = $db;
789	75					150	$self->{'_blsdb_length'} = $s;
790	75					145	$self->{'_blsdb_letters'} = $l;
791	75					133	last;
792							}
793							else {
794	80					69	chomp;
795	80					158	$db .= $_;
796							}
797							}
798							$self->element(
799							{
800	77					272	'Name' => 'BlastOutput_db',
801							'Data' => $db
802							}
803							);
804							}
805
806							# move inside of a hit
807							elsif (/^(?:Subject=\|>)\s(\S+)\s(.*)?/) {
808	922					1335	chomp;
809
810	922					3661	$self->debug("blast.pm: Hit: $1\n");
811	922	100				1586	$self->in_element('hsp')
812							&& $self->end_element( { 'Name' => 'Hsp' } );
813	922	100				2178	$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				2192	if ( !$self->within_element('result') ) {
		100
819	2					13	$self->_start_blastoutput;
820	2					7	$self->start_element( { 'Name' => 'Iteration' } );
821							}
822							elsif ( !$self->within_element('iteration') ) {
823	2					8	$self->start_element( { 'Name' => 'Iteration' } );
824							}
825	922					2338	$self->start_element( { 'Name' => 'Hit' } );
826	922					3299	my $id = $1;
827	922					1224	my $restofline = $2;
828
829	922					3746	$self->debug("Starting a hit: $1 $2\n");
830	922					2845	$self->element(
831							{
832							'Name' => 'Hit_id',
833							'Data' => $id
834							}
835							);
836	922					2856	my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
837	922					2250	$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					2470	while (my $v = shift @hit_signifs) {
851	642					778	my $tableid = $v->[2];
852	642	100				6860	if ($tableid !~ m{\Q$id\E}) {
853	14					48	$self->debug("Hit table ID $tableid doesn't match current hit id $id, checking next hit table entry...\n");
854	14					46	next HITTABLE;
855							} else {
856	628					1104	last HITTABLE;
857							}
858							}
859	922					3079	while ( defined( $_ = $self->_readline() ) ) {
860	1461	100				3790	next if (/^\s+$/);
861	1458					1375	chomp;
862	1458	100				3956	if (/Length\s=\s([\d,]+)/) {
863	922					1560	my $l = $1;
864	922					1006	$l =~ s/\,//g;
865	922					2610	$self->element(
866							{
867							'Name' => 'Hit_len',
868							'Data' => $l
869							}
870							);
871	922					1464	last;
872							}
873							else {
874	536	100				1179	if ($restofline !~ /\s$/) { # bug #3235
875	519					806	s/^\s(?!\s)/\x01/; #new line to concatenate desc lines with
876							}
877	536	100	100			2487	$restofline .= ($restofline =~ /\w$/ && $_ =~ /^\w/) ? " $_" : $_;
878	536					4188	$restofline =~ s/\s+/ /g; # this catches the newline as well
879	536					6433	$restofline =~ s/^ \| $//g;
880							}
881							}
882	922					5135	$restofline =~ s/\s+/ /g;
883	922					2205	$self->element(
884							{
885							'Name' => 'Hit_def',
886							'Data' => $restofline
887							}
888							);
889							}
890							elsif (/\s+(Plus\|Minus) Strand HSPs:/i) {
891	19					41	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				559	$self->in_element('hsp')
982							&& $self->end_element( { 'Name' => 'Hsp' } );
983	423					1111	$self->start_element( { 'Name' => 'Hsp' } );
984
985							# Some data clean-up so e-value will appear numeric to perl
986	423					1835	my ( $score, $bits, $evalue, $pvalue, $group ) =
987							( $1, $2, $3, $4, $5 );
988	423					627	$evalue =~ s/^e/1e/i;
989	423					372	$pvalue =~ s/^e/1e/i;
990
991	423					1052	$self->element(
992							{
993							'Name' => 'Hsp_score',
994							'Data' => $score
995							}
996							);
997	423					1008	$self->element(
998							{
999							'Name' => 'Hsp_bit-score',
1000							'Data' => $bits
1001							}
1002							);
1003	423					957	$self->element(
1004							{
1005							'Name' => 'Hsp_evalue',
1006							'Data' => $evalue
1007							}
1008							);
1009	423					865	$self->element(
1010							{
1011							'Name' => 'Hsp_pvalue',
1012							'Data' => $pvalue
1013							}
1014							);
1015
1016	423	100				1560	if ( defined $group ) {
1017	29					80	$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				180	$self->in_element('hsp')
1034							&& $self->end_element( { 'Name' => 'Hsp' } );
1035	70					143	my $featline;
1036	70					243	$_ = $self->_readline;
1037	70					365	while($_ !~ /^\s*$/) {
1038	119					200	chomp;
1039	119					255	$featline .= $_;
1040	119					265	$_ = $self->_readline;
1041							}
1042	70					288	$self->_pushback($_);
1043	70					879	$featline =~ s{(?:^\s+\|\s+^)}{}g;
1044	70					122	$featline =~ s{\n}{;}g;
1045	70					283	$self->start_element( { 'Name' => 'Hsp' } );
1046	70					424	$self->element(
1047							{
1048							'Name' => 'Hsp_features',
1049							'Data' => $featline
1050							}
1051							);
1052	70					365	$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				2397	if( !$self->{'_seen_hsp_features'} ) {
1064	1130	100				1459	$self->in_element('hsp')
1065							&& $self->end_element( { 'Name' => 'Hsp' } );
1066	1130					2952	$self->start_element( { 'Name' => 'Hsp' } );
1067							}
1068
1069							# Some data clean-up so e-value will appear numeric to perl
1070	1200					5015	my ( $bits, $score, $n, $evalue ) = ( $1, $2, $3, $4 );
1071	1200					1854	$evalue =~ s/^e/1e/i;
1072	1200					3599	$self->element(
1073							{
1074							'Name' => 'Hsp_score',
1075							'Data' => $score
1076							}
1077							);
1078	1200					3376	$self->element(
1079							{
1080							'Name' => 'Hsp_bit-score',
1081							'Data' => $bits
1082							}
1083							);
1084	1200					3046	$self->element(
1085							{
1086							'Name' => 'Hsp_evalue',
1087							'Data' => $evalue
1088							}
1089							);
1090	1200	100				2246	$self->element(
1091							{
1092							'Name' => 'Hsp_n',
1093							'Data' => $n
1094							}
1095							) if defined $n;
1096	1200	50				2160	$score = '' unless defined $score; # deal with BLAT which
1097							# has no score only bits
1098	1200					4422	$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					5816	$self->element(
1109							{
1110							'Name' => 'Hsp_identity',
1111							'Data' => $1
1112							}
1113							);
1114	1623					5123	$self->element(
1115							{
1116							'Name' => 'Hsp_align-len',
1117							'Data' => $2
1118							}
1119							);
1120	1623	100				2895	if ( defined $3 ) {
1121	1251					2867	$self->element(
1122							{
1123							'Name' => 'Hsp_positive',
1124							'Data' => $3
1125							}
1126							);
1127							}
1128							else {
1129	372					1099	$self->element(
1130							{
1131							'Name' => 'Hsp_positive',
1132							'Data' => $1
1133							}
1134							);
1135							}
1136	1623	100				3521	if ( defined $6 ) {
1137	787					2053	$self->element(
1138							{
1139							'Name' => 'Hsp_gaps',
1140							'Data' => $5
1141							}
1142							);
1143							}
1144
1145	1623					3536	$self->{'_Query'} = { 'begin' => 0, 'end' => 0 };
1146	1623					3018	$self->{'_Sbjct'} = { 'begin' => 0, 'end' => 0 };
1147
1148	1623	100				7393	if (/(Frame\s=\s.+)$/) {
1149
1150							# handle wu-blast Frame listing on same line
1151	44					100	$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				771	if (!defined($reporttype)) {
1160	2					6	$self->{'_reporttype'} = $reporttype = 'BLASTN';
1161	2					4	$bl2seq_fix = 1; # special case to resubmit the algorithm
1162							# reporttype
1163							}
1164	372					1103	next;
1165							}
1166							elsif ( $self->in_element('hsp')
1167							&& /Links\s=\s(\S+)/ox )
1168							{
1169	342					918	$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			787	my $frame1 = $1 \|\| 0;
1180	308		100			699	my $frame2 = $2 \|\| 0;
1181							# this is for bl2seq only
1182	308	100				466	if ( not defined $reporttype ) {
1183	4					6	$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					100	my $fh = $self->_fh;
1195	3					33	my $file_pos = tell $fh;
1196
1197	3					11	my $a_position = '';
1198	3					4	my $ali_length = '';
1199	3					6	my $b_position = '';
1200	3					17	while (my $line = <$fh>) {
1201	6	100				36	if ($line =~ m/^(?:Query\|Sbjct):?\s+(\-?\d+)?\s*(\S+)\s+(\-?\d+)?/) {
1202	3					9	$a_position = $1;
1203	3					8	my $alignment = $2;
1204	3					7	$b_position = $3;
1205
1206	12			12		4447	use Bio::LocatableSeq;
	12					21
	12					53110
1207	3					7	my $gap_symbols = $Bio::LocatableSeq::GAP_SYMBOLS;
1208	3					44	$alignment =~ s/[$gap_symbols]//g;
1209	3					6	$ali_length = length($alignment);
1210	3					8	last;
1211							}
1212							}
1213	3	100				17	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					16	seek $fh, $file_pos, 0;
1219							}
1220	4					11	$self->{'_reporttype'} = $reporttype;
1221							}
1222
1223	308					268	my ( $queryframe, $hitframe );
1224	308	100	66			698	if ( $reporttype eq 'TBLASTX' ) {
		100	33
		50
1225	187					204	( $queryframe, $hitframe ) = ( $frame1, $frame2 );
1226	187					488	$hitframe =~ s/\/\s*//g;
1227							}
1228							elsif ( $reporttype eq 'TBLASTN' \|\| $reporttype eq 'PSITBLASTN') {
1229	95					125	( $hitframe, $queryframe ) = ( $frame1, 0 );
1230							}
1231							elsif ( $reporttype eq 'BLASTX' \|\| $reporttype eq 'RPS-BLAST(BLASTP)') {
1232	26					34	( $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				49	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					828	'Name' => 'Hsp_query-frame',
1247							'Data' => $queryframe
1248							}
1249							);
1250
1251	308					776	$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					250	$self->debug("blast.pm: found parameters section \n");
1269
1270	83	100				205	$self->in_element('hsp')
1271							&& $self->end_element( { 'Name' => 'Hsp' } );
1272	83	100				281	$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				304	$self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1278	83	100				240	$self->within_element('iteration')
1279							&& $self->end_element( { 'Name' => 'Iteration' } );
1280
1281	83	50				299	next if /^\s+Subset/;
1282	83	100				420	my $blast = (/^(\s+Database\:)\|(\s*Lambda)/) ? 'ncbi' : 'wublast';
1283	83	50				289	if (/^\s*Histogram/) {
1284	0					0	$blast = 'btk';
1285							}
1286
1287	83					138	my $last = '';
1288
1289							# default is that gaps are allowed
1290	83					792	$self->element(
1291							{
1292							'Name' => 'Parameters_allowgaps',
1293							'Data' => 'yes'
1294							}
1295							);
1296	83					320	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			17811	if (/^\s+[\d+\.]+\s+[\d+\.]+\s+[\d+\.]/ and $last eq '') {
		100	66
		100
1300	15					60	$self->_pushback($_);
1301	15					43	$self->_pushback("Lambda K H\n");
1302	15					36	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					20	$self->_pushback($_);
1311
1312							# let's handle this in the loop
1313	5					24	last;
1314							}
1315							elsif (/^Query=/) {
1316	4					13	$self->_pushback($_);
1317	4	50				16	$self->_pushback($reportline) if $reportline;
1318	4					9	last PARSER;
1319							}
1320
1321							# here is where difference between wublast and ncbiblast
1322							# is better handled by different logic
1323	2795	100	100			12434	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					190	$c =~ s/\,//g;
1328	125					407	$self->element(
1329							{
1330							'Name' => 'Statistics_db-len',
1331							'Data' => $c
1332							}
1333							);
1334							}
1335							elsif (/letters in database:\s+(\-?[\d,]+)/i) {
1336	73					170	my $s = $1;
1337	73					203	$s =~ s/,//g;
1338	73					255	$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				8675	if (/E=(\S+)/) {
		100
		100
		100
		100
		100
		100
		100
		50
		100
		100
		100
		100
		100
1352	29					121	$self->element(
1353							{
1354							'Name' => 'Parameters_expect',
1355							'Data' => $1
1356							}
1357							);
1358							}
1359							elsif (/nogaps/) {
1360	2					9	$self->element(
1361							{
1362							'Name' => 'Parameters_allowgaps',
1363							'Data' => 'no'
1364							}
1365							);
1366							}
1367							elsif (/ctxfactor=(\S+)/) {
1368	26					110	$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					156	$self->element(
1380							{
1381							'Name' => "Parameters_$1",
1382							'Data' => 'yes'
1383							}
1384							);
1385							}
1386							elsif (/(\S+)=(\S+)/) {
1387	141					506	$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					31	my $firstgapinfo = 1;
1396	26					35	my $frame = undef;
1397	26		66			163	while ( defined($_) && !/^\s+$/ ) {
1398	108					224	s/^\s+//;
1399	108					341	s/\s+$//;
1400	108	100	100			385	if ( $firstgapinfo
1401							&& s/Q=(\d+),R=(\d+)\s+//x )
1402							{
1403	24					28	$firstgapinfo = 0;
1404
1405	24					95	$self->element(
1406							{
1407							'Name' => 'Parameters_gap-open',
1408							'Data' => $1
1409							}
1410							);
1411	24					102	$self->element(
1412							{
1413							'Name' => 'Parameters_gap-extend',
1414							'Data' => $2
1415							}
1416							);
1417	24					112	my @fields = split;
1418
1419	24					51	for my $type (
1420							qw(lambda_gapped
1421							kappa_gapped
1422							entropy_gapped)
1423							)
1424							{
1425	72	50				106	next if $type eq 'n/a';
1426	72	50				116	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					321	my ( $frameo, $matid, $matrix, @fields ) =
1441							split;
1442	84	100				161	if ( !defined $frame ) {
1443
1444							# keep some sort of default feature I guess
1445							# even though this is sort of wrong
1446	26					102	$self->element(
1447							{
1448							'Name' => 'Parameters_matrix',
1449							'Data' => $matrix
1450							}
1451							);
1452	26					107	$self->element(
1453							{
1454							'Name' => 'Statistics_lambda',
1455							'Data' => $fields[0]
1456							}
1457							);
1458	26					111	$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					83	$frame = $frameo;
1472	84					74	my $ii = 0;
1473	84					114	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					407	my $f = $fields[$ii];
1483	504	100				604	next unless defined $f; # deal with n/a
1484	456	100				520	if ( $f eq 'same' ) {
1485	72					98	$f = $fields[ $ii - 3 ];
1486							}
1487	456					282	$ii++;
1488	456					1044	$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					212	$_ = $self->_readline;
1501							}
1502	26					38	$last = $_;
1503							}
1504							elsif ( $last =~ /(Frame\|Strand)\s+MatID\s+Length/i ) {
1505	26					33	my $frame = undef;
1506	26		33			180	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					41	$last = $_;
1542							}
1543							elsif (/(\S+\s+\S+)\s+DFA:\s+(\S+)\s+$(.+)$/) {
1544	23	50				89	if ( $1 eq 'states in' ) {
		0
1545	23					128	$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					103	$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					104	my $cputype = lc($1);
1588	46					180	$self->element(
1589							{
1590							'Name' => "Statistics_$cputype\_cputime",
1591							'Data' => $2
1592							}
1593							);
1594	46					171	$self->element(
1595							{
1596							'Name' => "Statistics_$cputype\_actualtime",
1597							'Data' => $3
1598							}
1599							);
1600							}
1601							elsif (/^\s+Start:/) {
1602	23					204	my ( $junk, $start, $stime, $end, $etime ) =
1603							split( /\s+(Start\|End)\:\s+/, $_ );
1604	23					51	chomp($stime);
1605	23					81	$self->element(
1606							{
1607							'Name' => 'Statistics_starttime',
1608							'Data' => $stime
1609							}
1610							);
1611	23					34	chomp($etime);
1612	23					68	$self->element(
1613							{
1614							'Name' => 'Statistics_endtime',
1615							'Data' => $etime
1616							}
1617							);
1618							}
1619							elsif (/^\s+Database:\s+(.+)$/) {
1620	21					113	$self->element(
1621							{
1622							'Name' => 'Parameters_full_dbpath',
1623							'Data' => $1
1624							}
1625							);
1626
1627							}
1628							elsif (/^\s+Posted:\s+(.+)/) {
1629	21					47	my $d = $1;
1630	21					33	chomp($d);
1631	21					81	$self->element(
1632							{
1633							'Name' => 'Statistics_posted_date',
1634							'Data' => $d
1635							}
1636							);
1637							}
1638							}
1639							elsif ( $blast eq 'ncbi' ) {
1640
1641	1612	100				14404	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					284	$self->element(
1643							{
1644							'Name' => 'Parameters_matrix',
1645							'Data' => $1
1646							}
1647							);
1648							}
1649							elsif (/^Gapped/) {
1650	49					85	$gapped_stats = 1;
1651							}
1652							elsif (/^Lambda/) {
1653	107					241	$_ = $self->_readline;
1654	107					303	s/^\s+//;
1655	107					334	my ( $lambda, $kappa, $entropy ) = split;
1656	107	100				209	if ($gapped_stats) {
1657	49					177	$self->element(
1658							{
1659							'Name' => "Statistics_gapped_lambda",
1660							'Data' => $lambda
1661							}
1662							);
1663	49					182	$self->element(
1664							{
1665							'Name' => "Statistics_gapped_kappa",
1666							'Data' => $kappa
1667							}
1668							);
1669	49					134	$self->element(
1670							{
1671							'Name' => "Statistics_gapped_entropy",
1672							'Data' => $entropy
1673							}
1674							);
1675							}
1676							else {
1677	58					212	$self->element(
1678							{
1679							'Name' => "Statistics_lambda",
1680							'Data' => $lambda
1681							}
1682							);
1683	58					241	$self->element(
1684							{
1685							'Name' => "Statistics_kappa",
1686							'Data' => $kappa
1687							}
1688							);
1689	58					198	$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					419	$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					311	$self->element(
1719							{
1720							'Name' => 'Parameters_gap-open',
1721							'Data' => $1
1722							}
1723							);
1724	48					189	$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					203	$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					168	$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					32	$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					36	$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					29	$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					119	my $c = $1;
1773	48					84	$c =~ s/\,//g;
1774	48					203	$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					108	my $c = $1;
1783	51					150	$c =~ s/\,//g;
1784	51					222	$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					471	my $v = $2;
1796	325					341	chomp($v);
1797	325					1144	$self->element(
1798							{
1799							'Name' => "Statistics_$1",
1800							'Data' => $v
1801							}
1802							);
1803	325	100				760	if ( defined $4 ) {
1804	235					735	$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					76	$self->element(
1818							{
1819							'Name' => 'Statistics_framewindow',
1820							'Data' => $1
1821							}
1822							);
1823	14					53	$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					306	$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					265	$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					250	$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					276	$self->element(
1864							{
1865							'Name' => 'Parameters_expect',
1866							'Data' => $1
1867							}
1868							);
1869	52					193	$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					136	my $d = $1;
1878	50					67	chomp($d);
1879	50					208	$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					5706	$last = $_;
1891							}
1892							} elsif ( $self->in_element('hsp') ) {
1893	9058					18614	$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			15815	$self->{'_reporttype'} \|\|= $DEFAULTREPORTTYPE;
1897	9058					19013	my %data = (
1898							'Query' => '',
1899							'Mid' => '',
1900							'Hit' => ''
1901							);
1902	9058					6079	my $len;
1903	9058		66			28779	for ( my $i = 0 ; defined($_) && $i < 3 ; $i++ ) {
1904							# $self->debug("$i: $_") if $v;
1905	27150	50	66			127007	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					22638	chomp;
1913	27150	100				76834	if (/^((Query\|Sbjct):?\s+(\-?\d+)?\s*)(\S+)\s+(\-?\d+)?/) {
		50
1914	18116					45675	my ( $full, $type, $start, $str, $end ) =
1915							( $1, $2, $3, $4, $5 );
1916
1917	18116	100				20643	if ( $str eq '-' ) {
1918	48	100				83	$i = 3 if $type eq 'Sbjct';
1919							}
1920							else {
1921	18068					23406	$data{$type} = $str;
1922							}
1923	18116					12787	$len = length($full);
1924							$self->{"\_$type"}->{'begin'} = $start
1925	18116	100				38222	unless $self->{"_$type"}->{'begin'};
1926	18116					24305	$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			26807	$self->throw("no data for midline $_")
1933							unless ( defined $_ && defined $len );
1934	9034					14839	$data{'Mid'} = substr( $_, $len );
1935							}
1936	27150					45573	$_ = $self->_readline();
1937							}
1938							$self->characters(
1939							{
1940							'Name' => 'Hsp_qseq',
1941	9058					26814	'Data' => $data{'Query'}
1942							}
1943							);
1944							$self->characters(
1945							{
1946							'Name' => 'Hsp_hseq',
1947	9058					23007	'Data' => $data{'Sbjct'}
1948							}
1949							);
1950							$self->characters(
1951							{
1952							'Name' => 'Hsp_midline',
1953	9058					20298	'Data' => $data{'Mid'}
1954							}
1955							);
1956							}
1957							else {
1958							#$self->debug("blast.pm: unrecognized line $_");
1959							}
1960							}
1961
1962	102					399	$self->debug("blast.pm: End of BlastOutput\n");
1963	102	100				354	if ( $self->{'_seentop'} ) {
1964	96	100				231	$self->within_element('hsp')
1965							&& $self->end_element( { 'Name' => 'Hsp' } );
1966	96	100				237	$self->within_element('hit')
1967							&& $self->end_element( { 'Name' => 'Hit' } );
1968							# cleanup extra hits
1969	96	100				257	$self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1970	96	100				222	$self->within_element('iteration')
1971							&& $self->end_element( { 'Name' => 'Iteration' } );
1972	96	100				216	if ($bl2seq_fix) {
1973	6					24	$self->element(
1974							{
1975							'Name' => 'BlastOutput_program',
1976							'Data' => $reporttype
1977							}
1978							);
1979							}
1980	96					313	$self->end_element( { 'Name' => 'BlastOutput' } );
1981							}
1982	102					335	return $self->end_document();
1983							}
1984
1985							# Private method for internal use only.
1986							sub _start_blastoutput {
1987	96			96		149	my $self = shift;
1988	96					514	$self->start_element( { 'Name' => 'BlastOutput' } );
1989	96					192	$self->{'_seentop'} = 1;
1990	96					159	$self->{'_result_count'}++;
1991	96					133	$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		7078	my ( $self, $type ) = @_;
2034	8516					8600	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				16875	: ( $self->{'_will_handle_cache'}->{$type} =
2039							$handler->will_handle($type) );
2040
2041	8516	50				12917	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	3983	my ( $self, $data ) = @_;
2056
2057							# we currently don't care about attributes
2058	4258					3927	my $nm = $data->{'Name'};
2059	4258					4610	my $type = $MODEMAP{$nm};
2060	4258	50				6698	if ($type) {
2061	4258					5730	my $handler = $self->_will_handle($type);
2062	4258	50				7439	if ($handler) {
2063	4258					9642	my $func = sprintf( "start_%s", lc $type );
2064	4258					15355	$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					5240	unshift @{ $self->{'_elements'} }, $type;
	4258					6677
2073	4258	100				6074	if ( $type eq 'result' ) {
2074	96					200	$self->{'_values'} = {};
2075	96					202	$self->{'_result'} = undef;
2076							} else {
2077							# cleanup some things
2078	4162	50				7337	if ( defined $self->{'_values'} ) {
2079	4162					2931	foreach my $k (
2080	117000					201776	grep { /^\U$type\-/ }
2081	4162					16558	keys %{ $self->{'_values'} }
2082							)
2083							{
2084	38633					36655	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	30352	my ( $self, $data ) = @_;
2104
2105	43007					34088	my $nm = $data->{'Name'};
2106	43007					28339	my $type;
2107							my $rc;
2108							# cache these (TODO: we should probably cache all cross-report data)
2109	43007	100				53919	if ( $nm eq 'BlastOutput_program' ) {
2110	86	100				433	if ( $self->{'_last_data'} =~ /(t?blast[npx])/i ) {
2111	83					244	$self->{'_reporttype'} = uc $1;
2112							}
2113	86		66			219	$self->{'_reporttype'} \|\|= $DEFAULTREPORTTYPE;
2114							}
2115
2116	43007	100				47090	if ( $nm eq 'BlastOutput_version' ) {
2117	80					148	$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				47606	if ( $nm eq 'Hsp' ) {
2123	1623					2745	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				18448	) if defined $self->{'_last_hspdata'}->{$_};
2130							}
2131	1623					1892	$self->{'_last_hspdata'} = {};
2132							$self->element(
2133							{
2134							'Name' => 'Hsp_query-from',
2135	1623					5719	'Data' => $self->{'_Query'}->{'begin'}
2136							}
2137							);
2138							$self->element(
2139							{
2140							'Name' => 'Hsp_query-to',
2141	1623					4261	'Data' => $self->{'_Query'}->{'end'}
2142							}
2143							);
2144
2145							$self->element(
2146							{
2147							'Name' => 'Hsp_hit-from',
2148	1623					4622	'Data' => $self->{'_Sbjct'}->{'begin'}
2149							}
2150							);
2151							$self->element(
2152							{
2153							'Name' => 'Hsp_hit-to',
2154	1623					4259	'Data' => $self->{'_Sbjct'}->{'end'}
2155							}
2156							);
2157
2158							# } elsif( $nm eq 'Iteration' ) {
2159							# Nothing special needs to be done here.
2160							}
2161	43007	100				78985	if ( $type = $MODEMAP{$nm} ) {
		100
2162	4258					5990	my $handler = $self->_will_handle($type);
2163	4258	50				6575	if ($handler) {
2164	4258					10604	my $func = sprintf( "end_%s", lc $type );
2165	4258					13568	$rc = $handler->$func( $self->{'_reporttype'}, $self->{'_values'} );
2166							}
2167	4258					3144	shift @{ $self->{'_elements'} };
	4258					6180
2168
2169							}
2170							elsif ( $MAPPING{$nm} ) {
2171	38562	100				45950	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					2058	my $key = ( keys %{ $MAPPING{$nm} } )[0];
	2900					5766
2176							$self->{'_values'}->{$key}->{ $MAPPING{$nm}->{$key} } =
2177	2900					7579	$self->{'_last_data'};
2178							}
2179							else {
2180	35662					60113	$self->{'_values'}->{ $MAPPING{$nm} } = $self->{'_last_data'};
2181							}
2182							}
2183							else {
2184							#$self->debug("blast.pm: unknown nm $nm, ignoring\n");
2185							}
2186	43007					34121	$self->{'_last_data'} = ''; # remove read data if we are at
2187							# end of an element
2188	43007	100	100			70170	$self->{'_result'} = $rc if ( defined $type && $type eq 'result' );
2189	43007					31676	$self->{'_seen_hsp_features'} = 0;
2190	43007					49519	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	30593	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					37922	$self->characters($data);
2210	38749					43311	$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	50430	my ( $self, $data ) = @_;
2226	65923	100	100			64941	if ( $self->in_element('hsp')
2227							&& $data->{'Name'} =~ /^Hsp\_(qseq\|hseq\|midline)$/ )
2228							{
2229							$self->{'_last_hspdata'}->{ $data->{'Name'} } .= $data->{'Data'}
2230	32019	100				82882	if defined $data->{'Data'};
2231							}
2232	65923	100	100			218951	return unless ( defined $data->{'Data'} && $data->{'Data'} !~ /^\s+$/ );
2233	65853					96099	$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	2008	my ( $self, $name ) = @_;
2252							return 0
2253							if ( !defined $name && !defined $self->{'_elements'}
2254	2214	100	33			4850	\|\| scalar @{ $self->{'_elements'} } == 0 );
	2214		66			5512
2255	2212					1702	foreach ( @{ $self->{'_elements'} } ) {
	2212					3471
2256	3134	100				4751	if ( $_ eq $name ) {
2257	1949					4562	return 1;
2258							}
2259							}
2260	263					500	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	194503	my ( $self, $name ) = @_;
2279	264029	100				320818	return 0 if !defined $self->{'_elements'}->[0];
2280	263910					989808	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	138	my ($self) = @_;
2296	102					205	$self->{'_lasttype'} = '';
2297	102					199	$self->{'_values'} = {};
2298	102					414	$self->{'_result'} = undef;
2299	102					216	$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	154	my ( $self, @args ) = @_;
2315
2316							#$self->debug("blast.pm: end_document\n");
2317	102					575	return $self->{'_result'};
2318							}
2319
2320							sub write_result {
2321	5			5	1	24	my ( $self, $blast, @args ) = @_;
2322
2323	5	50				14	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					31	$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		43	my ($self, $hits) = @_;
2424	8					21	while ( my $v = shift @{ $hits }) {
	1532					3043
2425	1524	50				2137	next unless defined $v;
2426	1524					3601	$self->start_element( { 'Name' => 'Hit' } );
2427	1524					3586	my $id = $v->[2];
2428	1524					1437	my $desc = $v->[3];
2429	1524					4055	$self->element(
2430							{
2431							'Name' => 'Hit_id',
2432							'Data' => $id
2433							}
2434							);
2435	1524					3719	my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
2436	1524					3135	$self->element(
2437							{
2438							'Name' => 'Hit_accession',
2439							'Data' => $acc
2440							}
2441							);
2442	1524	50				2725	if ( defined $v ) {
2443	1524					3343	$self->element(
2444							{
2445							'Name' => 'Hit_signif',
2446							'Data' => $v->[0]
2447							}
2448							);
2449	1524	50				2360	if (exists $self->{'_wublast'}) {
2450	0					0	$self->element(
2451							{
2452							'Name' => 'Hit_score',
2453							'Data' => $v->[1]
2454							}
2455							);
2456							} else {
2457	1524					3175	$self->element(
2458							{
2459							'Name' => 'Hit_bits',
2460							'Data' => $v->[1]
2461							}
2462							);
2463							}
2464
2465							}
2466							$self->element(
2467							{
2468	1524					3300	'Name' => 'Hit_def',
2469							'Data' => $desc
2470							}
2471							);
2472	1524					3017	$self->end_element( { 'Name' => 'Hit' } );
2473							}
2474							}
2475
2476
2477							1;
2478
2479							__END__