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		4259	use Bio::SearchIO::IteratedSearchResultEventBuilder;
	12					13
	12					283
145	12			12		52	use strict;
	12					11
	12					233
146	12					806	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		35	);
	12					11
153
154
155	12			12		38	use base qw(Bio::SearchIO);
	12					11
	12					689
156	12			12		48	use Data::Dumper;
	12					15
	12					4889
157
158							BEGIN {
159
160							# mapping of NCBI Blast terms to Bioperl hash keys
161	12			12		50	%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					512	%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					45	for my $frame ( 0 .. 3 ) {
278	48					46	for my $strand ( '+', '-' ) {
279	96					81	for my $ind (
280							qw(length efflength E S W T X X_gapped E2
281							E2_gapped S2)
282							)
283							{
284	1056					2003	$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					213	for my $type (qw(used computed gapped)) {
289	864					792	my $key = "Statistics_frame$strand$frame\_$val\_$type";
290	864					1150	my $val =
291							{ 'RESULT-statistics' =>
292							"Frame$strand$frame\_$val\_$type" };
293	864					1029	$MAPPING{$key} = $val;
294							}
295							}
296							}
297							}
298
299							# add Statistics
300	12					20	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					207	my $key = "Statistics_$stats";
312	264					231	my $val = { 'RESULT-statistics' => $stats };
313	264					390	$MAPPING{$key} = $val;
314							}
315
316							# add WU-BLAST Parameters
317	12					18	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					279	my $key = "Parameters_$param";
326	348					315	my $val = { 'RESULT-parameters' => $param };
327	348					398	$MAPPING{$key} = $val;
328							}
329
330	12					13	$DEFAULT_BLAST_WRITER_CLASS = 'Bio::SearchIO::Writer::HitTableWriter';
331	12					15	$MAX_HSP_OVERLAP = 2; # Used when tiling multiple HSPs.
332	12					34545	$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		243	my ( $self, @args ) = @_;
395	92					356	$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					646	my $handler = Bio::SearchIO::IteratedSearchResultEventBuilder->new(@args);
403	92					219	$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					370	my ($rpttype ) = $self->_rearrange(
413							[
414							qw(
415							REPORT_TYPE)
416							],
417							@args
418							);
419	92	100				336	defined $rpttype && ( $self->{'_reporttype'} = $rpttype );
420							}
421
422							sub attach_EventHandler {
423	184			184	1	259	my ($self,$handler) = @_;
424
425	184					604	$self->SUPER::attach_EventHandler($handler);
426
427							# Optimization: caching the EventHandler since it is used a lot
428							# during the parse.
429
430	184					252	$self->{'_handler_cache'} = $handler;
431	184					297	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	3017	my ($self) = @_;
446	102					285	my $v = $self->verbose;
447	102					182	my $data = '';
448	102					143	my $flavor = '';
449	102					176	$self->{'_seentop'} = 0; # start next report at top
450
451	102					132	my ( $reporttype, $seenquery, $reportline, $reportversion );
452	0					0	my ( $seeniteration, $found_again );
453	102					365	my $incl_threshold = $self->inclusion_threshold;
454	102					141	my $bl2seq_fix;
455	102					309	$self->start_document(); # let the fun begin...
456	102					111	my (@hit_signifs);
457	102					141	my $gapped_stats = 0; # for switching between gapped/ungapped
458							# lambda, K, H
459	102					183	local $_ = "\n"; #consistency
460							PARSER:
461	102					384	while ( defined( $_ = $self->_readline ) ) {
462	20670	100				47813	next if (/^\s+$/); # skip empty lines
463	15331	100				21566	next if (/CPU time:/);
464	15329	50				18153	next if (/^>\s*$/);
465	15329	100				19866	next if (/[]+\s+No hits found\s+[]+/);
466	15325	100	66			183891	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					317	($reporttype, $reportversion) = ($1, $2);
473							# need to keep track of whether this is WU-BLAST
474	85	50	33			413	if ($reportversion && $reportversion =~ m{WashU$}) {
475	0					0	$self->{'_wublast'}++;
476							}
477	85					536	$self->debug("blast.pm: Start of new report: $reporttype, $reportversion\n");
478	85	100				262	if ( $self->{'_seentop'} ) {
479							# This handles multi-result input streams
480	5					15	$self->_pushback($_);
481	5					11	last PARSER;
482							}
483	80					285	$self->_start_blastoutput;
484	80	100				242	if ($reporttype =~ /RPS-BLAST/) {
485	4					6	$reporttype .= '(BLASTP)'; # default RPS-BLAST type
486							}
487	80					112	$reportline = $_; # to fix the fact that RPS-BLAST output is wrong
488	80					412	$self->element(
489							{
490							'Name' => 'BlastOutput_program',
491							'Data' => $reporttype
492							}
493							);
494
495	80					321	$self->element(
496							{
497							'Name' => 'BlastOutput_version',
498							'Data' => $reportversion
499							}
500							);
501	80					303	$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					252	my $algorithm_reference = "$1\n";
512	73					185	$_ = $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			601	while($_ !~ /^$/ && $_ !~ /^RID:/ && $_ !~ /^Database:/ && $_ !~ /^Query=/) {
			66
			66
518	146					313	$algorithm_reference .= "$_";
519	146					270	$_ = $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					258	$self->_pushback($_);
524	73					244	$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					27	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				311	next unless $1 =~ /Results from round/;
544	2					9	$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				8	$self->in_element('hit')
549							&& $self->end_element( { 'Name' => 'Hit' } );
550	2	100				7	if ( defined $seeniteration ) {
551	1	50				3	$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					11	$seeniteration = 1;
559							}
560							elsif (/^Query=\s(.)$/) {
561	102					265	my $q = $1;
562	102					474	$self->debug("blast.pm: Query= found...$_\n");
563	102					155	my $size = 0;
564	102	100				249	if ( defined $seenquery ) {
565	8					24	$self->_pushback($_);
566	8	100				26	$self->_pushback($reportline) if $reportline;
567	8					16	last PARSER;
568							}
569	94	100				254	if ( !defined $reporttype ) {
570	14					52	$self->_start_blastoutput;
571	14	50				41	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					44	$self->start_element( { 'Name' => 'Iteration' } );
578							}
579	14					31	$seeniteration = 1;
580							}
581	94					149	$seenquery = $q;
582	94					212	$_ = $self->_readline;
583	94					223	while ( defined($_) ) {
584	142	50				322	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			870	if ( /$(\-?[\d,]+)\s+letters.*$/ \|\| /^Length=(\-?[\d,]+)/ ) {
591	94					191	$size = $1;
592	94					160	$size =~ s/,//g;
593	94					149	last;
594							}
595							else {
596							# bug 2391
597	48	100	100			352	$q .= ($q =~ /\w$/ && $_ =~ /^\w/) ? " $_" : $_;
598	48					335	$q =~ s/\s+/ /g; # this catches the newline as well
599	48					482	$q =~ s/^ \| $//g;
600							}
601
602	48					121	$_ = $self->_readline;
603							}
604	94					207	chomp($q);
605	94					404	my ( $nm, $desc ) = split( /\s+/, $q, 2 );
606	94	100				451	$self->element(
607							{
608							'Name' => 'BlastOutput_query-def',
609							'Data' => $nm
610							}
611							) if $nm;
612	94					359	$self->element(
613							{
614							'Name' => 'BlastOutput_query-len',
615							'Data' => $size
616							}
617							);
618	94	100				413	defined $desc && $desc =~ s/\s+$//;
619	94					309	$self->element(
620							{
621							'Name' => 'BlastOutput_querydesc',
622							'Data' => $desc
623							}
624							);
625	94					490	my ( $gi, $acc, $version ) = $self->_get_seq_identifiers($nm);
626	94	100	66			407	$version = defined($version) && length($version) ? ".$version" : "";
627	94	100				474	$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				328	) if $self->{'_blsdb_length'};
642							$self->element(
643							{
644							'Name' => 'BlastOutput_db-let',
645							'Data' => $self->{'_blsdb_letters'}
646							}
647	94	100				305	) if $self->{'_blsdb_letters'};
648							$self->element(
649							{
650							'Name' => 'BlastOutput_db',
651							'Data' => $self->{'_blsdb'}
652							}
653	94	100				417	) 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					128	$self->debug("Bypassing features line: $_");
660	69					96	$_ = $self->_readline;
661							}
662	1					4	$self->_pushback($_);
663							}
664							elsif (/Sequences producing significant alignments:/) {
665	56					211	$self->debug("blast.pm: Processing NCBI-BLAST descriptions\n");
666	56					86	$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				117	if ( !$self->in_element('iteration') ) {
677	50					168	$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					144	my $h_regex;
683							my $seen_block;
684							DESCLINE:
685	56					169	while ( defined( my $descline = $self->_readline() ) ) {
686	1354	100				2699	if ($descline =~ m{^\s*$}) {
687	105	100				299	last DESCLINE if $seen_block;
688	55					177	next DESCLINE;
689							}
690							# any text match is part of block...
691	1249					914	$seen_block++;
692							# GCG multiline oddness...
693	1249	100				1737	if ($descline =~ /^(\S+)\s+Begin:\s\d+\s+End:\s+\d+/xms) {
694	3					7	my ($id, $nextline) = ($1, $self->_readline);
695	3					8	$nextline =~ s{^!}{};
696	3					5	$descline = "$id $nextline";
697							}
698							# NCBI style hit table (no N)
699	1249	100				208017	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					1916	my ( $score, $evalue ) = ($1, $2);
706
707							# Some data clean-up so e-value will appear numeric to perl
708	1243					1397	$evalue =~ s/^e/1e/i;
709
710							# This to handle no-HSP case
711	1243					3053	my @line = split ' ',$descline;
712
713							# we want to throw away the score, evalue
714	1243					945	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				2238	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					5613	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					28	$self->_pushback($descline); # Catch leading > (end of section)
732	6					21	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					62	$_ = $self->_readline();
742	26					42	$flavor = 'wu';
743
744	26	100				75	if ( !$self->in_element('iteration') ) {
745	23					74	$self->start_element( { 'Name' => 'Iteration' } );
746							}
747
748	26		66			98	while ( defined( $_ = $self->_readline() )
749							&& !/^\s+$/ )
750							{
751	923					2089	my @line = split;
752	923					583	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					3269	push @hit_signifs,
757							[ pop @line, pop @line, shift @line, join( ' ', @line ) ];
758							}
759
760							}
761							elsif (/^Database:\s(.+?)\s$/) {
762
763	77					429	$self->debug("blast.pm: Database: $1\n");
764	77					171	my $db = $1;
765	77					184	while ( defined( $_ = $self->_readline ) ) {
766	155	100				585	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					195	my ( $s, $l ) = ( $1, $2 );
773	75					195	$s =~ s/,//g;
774	75					207	$l =~ s/,//g;
775	75					318	$self->element(
776							{
777							'Name' => 'BlastOutput_db-len',
778							'Data' => $s
779							}
780							);
781	75					281	$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					180	$self->{'_blsdb'} = $db;
789	75					143	$self->{'_blsdb_length'} = $s;
790	75					169	$self->{'_blsdb_letters'} = $l;
791	75					137	last;
792							}
793							else {
794	80					71	chomp;
795	80					178	$db .= $_;
796							}
797							}
798							$self->element(
799							{
800	77					267	'Name' => 'BlastOutput_db',
801							'Data' => $db
802							}
803							);
804							}
805
806							# move inside of a hit
807							elsif (/^(?:Subject=\|>)\s(\S+)\s(.*)?/) {
808	922					1402	chomp;
809
810	922					3736	$self->debug("blast.pm: Hit: $1\n");
811	922	100				1680	$self->in_element('hsp')
812							&& $self->end_element( { 'Name' => 'Hsp' } );
813	922	100				2188	$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				2209	if ( !$self->within_element('result') ) {
		100
819	2					10	$self->_start_blastoutput;
820	2					8	$self->start_element( { 'Name' => 'Iteration' } );
821							}
822							elsif ( !$self->within_element('iteration') ) {
823	2					11	$self->start_element( { 'Name' => 'Iteration' } );
824							}
825	922					2388	$self->start_element( { 'Name' => 'Hit' } );
826	922					2814	my $id = $1;
827	922					1129	my $restofline = $2;
828
829	922					3791	$self->debug("Starting a hit: $1 $2\n");
830	922					2777	$self->element(
831							{
832							'Name' => 'Hit_id',
833							'Data' => $id
834							}
835							);
836	922					3113	my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
837	922					2358	$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					2703	while (my $v = shift @hit_signifs) {
851	642					748	my $tableid = $v->[2];
852	642	100				7106	if ($tableid !~ m{\Q$id\E}) {
853	14					49	$self->debug("Hit table ID $tableid doesn't match current hit id $id, checking next hit table entry...\n");
854	14					54	next HITTABLE;
855							} else {
856	628					1163	last HITTABLE;
857							}
858							}
859	922					2925	while ( defined( $_ = $self->_readline() ) ) {
860	1461	100				4155	next if (/^\s+$/);
861	1458					1417	chomp;
862	1458	100				3828	if (/Length\s=\s([\d,]+)/) {
863	922					1409	my $l = $1;
864	922					1166	$l =~ s/\,//g;
865	922					2365	$self->element(
866							{
867							'Name' => 'Hit_len',
868							'Data' => $l
869							}
870							);
871	922					1594	last;
872							}
873							else {
874	536	100				1272	if ($restofline !~ /\s$/) { # bug #3235
875	519					832	s/^\s(?!\s)/\x01/; #new line to concatenate desc lines with
876							}
877	536	100	100			2571	$restofline .= ($restofline =~ /\w$/ && $_ =~ /^\w/) ? " $_" : $_;
878	536					4303	$restofline =~ s/\s+/ /g; # this catches the newline as well
879	536					6468	$restofline =~ s/^ \| $//g;
880							}
881							}
882	922					5164	$restofline =~ s/\s+/ /g;
883	922					2196	$self->element(
884							{
885							'Name' => 'Hit_def',
886							'Data' => $restofline
887							}
888							);
889							}
890							elsif (/\s+(Plus\|Minus) Strand HSPs:/i) {
891	19					44	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				549	$self->in_element('hsp')
982							&& $self->end_element( { 'Name' => 'Hsp' } );
983	423					1538	$self->start_element( { 'Name' => 'Hsp' } );
984
985							# Some data clean-up so e-value will appear numeric to perl
986	423					2012	my ( $score, $bits, $evalue, $pvalue, $group ) =
987							( $1, $2, $3, $4, $5 );
988	423					662	$evalue =~ s/^e/1e/i;
989	423					439	$pvalue =~ s/^e/1e/i;
990
991	423					1090	$self->element(
992							{
993							'Name' => 'Hsp_score',
994							'Data' => $score
995							}
996							);
997	423					1072	$self->element(
998							{
999							'Name' => 'Hsp_bit-score',
1000							'Data' => $bits
1001							}
1002							);
1003	423					935	$self->element(
1004							{
1005							'Name' => 'Hsp_evalue',
1006							'Data' => $evalue
1007							}
1008							);
1009	423					995	$self->element(
1010							{
1011							'Name' => 'Hsp_pvalue',
1012							'Data' => $pvalue
1013							}
1014							);
1015
1016	423	100				1549	if ( defined $group ) {
1017	29					69	$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				154	$self->in_element('hsp')
1034							&& $self->end_element( { 'Name' => 'Hsp' } );
1035	70					120	my $featline;
1036	70					197	$_ = $self->_readline;
1037	70					305	while($_ !~ /^\s*$/) {
1038	119					120	chomp;
1039	119					198	$featline .= $_;
1040	119					218	$_ = $self->_readline;
1041							}
1042	70					173	$self->_pushback($_);
1043	70					740	$featline =~ s{(?:^\s+\|\s+^)}{}g;
1044	70					112	$featline =~ s{\n}{;}g;
1045	70					210	$self->start_element( { 'Name' => 'Hsp' } );
1046	70					307	$self->element(
1047							{
1048							'Name' => 'Hsp_features',
1049							'Data' => $featline
1050							}
1051							);
1052	70					276	$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				2327	if( !$self->{'_seen_hsp_features'} ) {
1064	1130	100				1617	$self->in_element('hsp')
1065							&& $self->end_element( { 'Name' => 'Hsp' } );
1066	1130					3096	$self->start_element( { 'Name' => 'Hsp' } );
1067							}
1068
1069							# Some data clean-up so e-value will appear numeric to perl
1070	1200					4957	my ( $bits, $score, $n, $evalue ) = ( $1, $2, $3, $4 );
1071	1200					1883	$evalue =~ s/^e/1e/i;
1072	1200					3785	$self->element(
1073							{
1074							'Name' => 'Hsp_score',
1075							'Data' => $score
1076							}
1077							);
1078	1200					3304	$self->element(
1079							{
1080							'Name' => 'Hsp_bit-score',
1081							'Data' => $bits
1082							}
1083							);
1084	1200					2985	$self->element(
1085							{
1086							'Name' => 'Hsp_evalue',
1087							'Data' => $evalue
1088							}
1089							);
1090	1200	100				2316	$self->element(
1091							{
1092							'Name' => 'Hsp_n',
1093							'Data' => $n
1094							}
1095							) if defined $n;
1096	1200	50				1647	$score = '' unless defined $score; # deal with BLAT which
1097							# has no score only bits
1098	1200					3864	$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					5696	$self->element(
1109							{
1110							'Name' => 'Hsp_identity',
1111							'Data' => $1
1112							}
1113							);
1114	1623					4961	$self->element(
1115							{
1116							'Name' => 'Hsp_align-len',
1117							'Data' => $2
1118							}
1119							);
1120	1623	100				3126	if ( defined $3 ) {
1121	1251					2815	$self->element(
1122							{
1123							'Name' => 'Hsp_positive',
1124							'Data' => $3
1125							}
1126							);
1127							}
1128							else {
1129	372					917	$self->element(
1130							{
1131							'Name' => 'Hsp_positive',
1132							'Data' => $1
1133							}
1134							);
1135							}
1136	1623	100				3302	if ( defined $6 ) {
1137	787					1984	$self->element(
1138							{
1139							'Name' => 'Hsp_gaps',
1140							'Data' => $5
1141							}
1142							);
1143							}
1144
1145	1623					3757	$self->{'_Query'} = { 'begin' => 0, 'end' => 0 };
1146	1623					3260	$self->{'_Sbjct'} = { 'begin' => 0, 'end' => 0 };
1147
1148	1623	100				6799	if (/(Frame\s=\s.+)$/) {
1149
1150							# handle wu-blast Frame listing on same line
1151	44					111	$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				678	if (!defined($reporttype)) {
1160	2					6	$self->{'_reporttype'} = $reporttype = 'BLASTN';
1161	2					3	$bl2seq_fix = 1; # special case to resubmit the algorithm
1162							# reporttype
1163							}
1164	372					1077	next;
1165							}
1166							elsif ( $self->in_element('hsp')
1167							&& /Links\s=\s(\S+)/ox )
1168							{
1169	342					975	$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			798	my $frame1 = $1 \|\| 0;
1180	308		100			700	my $frame2 = $2 \|\| 0;
1181							# this is for bl2seq only
1182	308	100				465	if ( not defined $reporttype ) {
1183	4					7	$bl2seq_fix = 1;
1184	4	100	66			25	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					101	my $fh = $self->_fh;
1195	3					30	my $file_pos = tell $fh;
1196
1197	3					5	my $a_position = '';
1198	3					6	my $ali_length = '';
1199	3					5	my $b_position = '';
1200	3					14	while (my $line = <$fh>) {
1201	6	100				27	if ($line =~ m/^(?:Query\|Sbjct):?\s+(\-?\d+)?\s*(\S+)\s+(\-?\d+)?/) {
1202	3					8	$a_position = $1;
1203	3					7	my $alignment = $2;
1204	3					6	$b_position = $3;
1205
1206	12			12		4327	use Bio::LocatableSeq;
	12					20
	12					53331
1207	3					4	my $gap_symbols = $Bio::LocatableSeq::GAP_SYMBOLS;
1208	3					42	$alignment =~ s/[$gap_symbols]//g;
1209	3					5	$ali_length = length($alignment);
1210	3					6	last;
1211							}
1212							}
1213	3	100				15	my $coord_length = ($a_position < $b_position) ? ($b_position - $a_position + 1)
1214							: ($a_position - $b_position + 1);
1215	3	100				12	($coord_length == ($ali_length * 3)) ? ($reporttype = 'BLASTX') : ($reporttype = 'TBLASTN');
1216
1217							# Rewind filehandle to its original position to continue parsing
1218	3					15	seek $fh, $file_pos, 0;
1219							}
1220	4					12	$self->{'_reporttype'} = $reporttype;
1221							}
1222
1223	308					300	my ( $queryframe, $hitframe );
1224	308	100	66			688	if ( $reporttype eq 'TBLASTX' ) {
		100	33
		50
1225	187					195	( $queryframe, $hitframe ) = ( $frame1, $frame2 );
1226	187					512	$hitframe =~ s/\/\s*//g;
1227							}
1228							elsif ( $reporttype eq 'TBLASTN' \|\| $reporttype eq 'PSITBLASTN') {
1229	95					123	( $hitframe, $queryframe ) = ( $frame1, 0 );
1230							}
1231							elsif ( $reporttype eq 'BLASTX' \|\| $reporttype eq 'RPS-BLAST(BLASTP)') {
1232	26					36	( $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				48	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					848	'Name' => 'Hsp_query-frame',
1247							'Data' => $queryframe
1248							}
1249							);
1250
1251	308					851	$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					252	$self->debug("blast.pm: found parameters section \n");
1269
1270	83	100				200	$self->in_element('hsp')
1271							&& $self->end_element( { 'Name' => 'Hsp' } );
1272	83	100				264	$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				308	$self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1278	83	100				226	$self->within_element('iteration')
1279							&& $self->end_element( { 'Name' => 'Iteration' } );
1280
1281	83	50				296	next if /^\s+Subset/;
1282	83	100				393	my $blast = (/^(\s+Database\:)\|(\s*Lambda)/) ? 'ncbi' : 'wublast';
1283	83	50				266	if (/^\s*Histogram/) {
1284	0					0	$blast = 'btk';
1285							}
1286
1287	83					130	my $last = '';
1288
1289							# default is that gaps are allowed
1290	83					560	$self->element(
1291							{
1292							'Name' => 'Parameters_allowgaps',
1293							'Data' => 'yes'
1294							}
1295							);
1296	83					304	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			17961	if (/^\s+[\d+\.]+\s+[\d+\.]+\s+[\d+\.]/ and $last eq '') {
		100	66
		100
1300	15					50	$self->_pushback($_);
1301	15					59	$self->_pushback("Lambda K H\n");
1302	15					29	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					16	$self->_pushback($_);
1311
1312							# let's handle this in the loop
1313	5					20	last;
1314							}
1315							elsif (/^Query=/) {
1316	4					10	$self->_pushback($_);
1317	4	50				15	$self->_pushback($reportline) if $reportline;
1318	4					8	last PARSER;
1319							}
1320
1321							# here is where difference between wublast and ncbiblast
1322							# is better handled by different logic
1323	2795	100	100			12697	if ( /Number of Sequences:\s+([\d\,]+)/i
		100
		50
		100
		50
1324							\|\| /of sequences in database:\s+(\-?[\d,]+)/i )
1325							{
1326	125					216	my $c = $1;
1327	125					218	$c =~ s/\,//g;
1328	125					427	$self->element(
1329							{
1330							'Name' => 'Statistics_db-len',
1331							'Data' => $c
1332							}
1333							);
1334							}
1335							elsif (/letters in database:\s+(\-?[\d,]+)/i) {
1336	73					144	my $s = $1;
1337	73					206	$s =~ s/,//g;
1338	73					245	$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				8700	if (/E=(\S+)/) {
		100
		100
		100
		100
		100
		100
		100
		50
		100
		100
		100
		100
		100
1352	29					124	$self->element(
1353							{
1354							'Name' => 'Parameters_expect',
1355							'Data' => $1
1356							}
1357							);
1358							}
1359							elsif (/nogaps/) {
1360	2					7	$self->element(
1361							{
1362							'Name' => 'Parameters_allowgaps',
1363							'Data' => 'no'
1364							}
1365							);
1366							}
1367							elsif (/ctxfactor=(\S+)/) {
1368	26					100	$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					155	$self->element(
1380							{
1381							'Name' => "Parameters_$1",
1382							'Data' => 'yes'
1383							}
1384							);
1385							}
1386							elsif (/(\S+)=(\S+)/) {
1387	141					501	$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					41	my $firstgapinfo = 1;
1396	26					35	my $frame = undef;
1397	26		66			177	while ( defined($_) && !/^\s+$/ ) {
1398	108					248	s/^\s+//;
1399	108					360	s/\s+$//;
1400	108	100	100			404	if ( $firstgapinfo
1401							&& s/Q=(\d+),R=(\d+)\s+//x )
1402							{
1403	24					32	$firstgapinfo = 0;
1404
1405	24					102	$self->element(
1406							{
1407							'Name' => 'Parameters_gap-open',
1408							'Data' => $1
1409							}
1410							);
1411	24					96	$self->element(
1412							{
1413							'Name' => 'Parameters_gap-extend',
1414							'Data' => $2
1415							}
1416							);
1417	24					111	my @fields = split;
1418
1419	24					47	for my $type (
1420							qw(lambda_gapped
1421							kappa_gapped
1422							entropy_gapped)
1423							)
1424							{
1425	72	50				116	next if $type eq 'n/a';
1426	72	50				144	if ( !@fields ) {
1427	0					0	warn "fields is empty for $type\n";
1428	0					0	next;
1429							}
1430							$self->element(
1431							{
1432	72					211	'Name' =>
1433							"Statistics_frame$frame\_$type",
1434							'Data' => shift @fields
1435							}
1436							);
1437							}
1438							}
1439							else {
1440	84					364	my ( $frameo, $matid, $matrix, @fields ) =
1441							split;
1442	84	100				158	if ( !defined $frame ) {
1443
1444							# keep some sort of default feature I guess
1445							# even though this is sort of wrong
1446	26					97	$self->element(
1447							{
1448							'Name' => 'Parameters_matrix',
1449							'Data' => $matrix
1450							}
1451							);
1452	26					101	$self->element(
1453							{
1454							'Name' => 'Statistics_lambda',
1455							'Data' => $fields[0]
1456							}
1457							);
1458	26					103	$self->element(
1459							{
1460							'Name' => 'Statistics_kappa',
1461							'Data' => $fields[1]
1462							}
1463							);
1464	26					88	$self->element(
1465							{
1466							'Name' => 'Statistics_entropy',
1467							'Data' => $fields[2]
1468							}
1469							);
1470							}
1471	84					88	$frame = $frameo;
1472	84					63	my $ii = 0;
1473	84					119	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					403	my $f = $fields[$ii];
1483	504	100				626	next unless defined $f; # deal with n/a
1484	456	100				600	if ( $f eq 'same' ) {
1485	72					90	$f = $fields[ $ii - 3 ];
1486							}
1487	456					302	$ii++;
1488	456					1041	$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					240	$_ = $self->_readline;
1501							}
1502	26					41	$last = $_;
1503							}
1504							elsif ( $last =~ /(Frame\|Strand)\s+MatID\s+Length/i ) {
1505	26					39	my $frame = undef;
1506	26		33			150	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					38	$last = $_;
1542							}
1543							elsif (/(\S+\s+\S+)\s+DFA:\s+(\S+)\s+$(.+)$/) {
1544	23	50				73	if ( $1 eq 'states in' ) {
		0
1545	23					142	$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					122	$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					116	my $cputype = lc($1);
1588	46					181	$self->element(
1589							{
1590							'Name' => "Statistics_$cputype\_cputime",
1591							'Data' => $2
1592							}
1593							);
1594	46					178	$self->element(
1595							{
1596							'Name' => "Statistics_$cputype\_actualtime",
1597							'Data' => $3
1598							}
1599							);
1600							}
1601							elsif (/^\s+Start:/) {
1602	23					213	my ( $junk, $start, $stime, $end, $etime ) =
1603							split( /\s+(Start\|End)\:\s+/, $_ );
1604	23					45	chomp($stime);
1605	23					78	$self->element(
1606							{
1607							'Name' => 'Statistics_starttime',
1608							'Data' => $stime
1609							}
1610							);
1611	23					34	chomp($etime);
1612	23					66	$self->element(
1613							{
1614							'Name' => 'Statistics_endtime',
1615							'Data' => $etime
1616							}
1617							);
1618							}
1619							elsif (/^\s+Database:\s+(.+)$/) {
1620	21					99	$self->element(
1621							{
1622							'Name' => 'Parameters_full_dbpath',
1623							'Data' => $1
1624							}
1625							);
1626
1627							}
1628							elsif (/^\s+Posted:\s+(.+)/) {
1629	21					49	my $d = $1;
1630	21					34	chomp($d);
1631	21					74	$self->element(
1632							{
1633							'Name' => 'Statistics_posted_date',
1634							'Data' => $d
1635							}
1636							);
1637							}
1638							}
1639							elsif ( $blast eq 'ncbi' ) {
1640
1641	1612	100				13834	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					298	$self->element(
1643							{
1644							'Name' => 'Parameters_matrix',
1645							'Data' => $1
1646							}
1647							);
1648							}
1649							elsif (/^Gapped/) {
1650	49					74	$gapped_stats = 1;
1651							}
1652							elsif (/^Lambda/) {
1653	107					231	$_ = $self->_readline;
1654	107					330	s/^\s+//;
1655	107					332	my ( $lambda, $kappa, $entropy ) = split;
1656	107	100				220	if ($gapped_stats) {
1657	49					172	$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					145	$self->element(
1670							{
1671							'Name' => "Statistics_gapped_entropy",
1672							'Data' => $entropy
1673							}
1674							);
1675							}
1676							else {
1677	58					250	$self->element(
1678							{
1679							'Name' => "Statistics_lambda",
1680							'Data' => $lambda
1681							}
1682							);
1683	58					306	$self->element(
1684							{
1685							'Name' => "Statistics_kappa",
1686							'Data' => $kappa
1687							}
1688							);
1689	58					204	$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					260	$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					215	$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					278	$self->element(
1719							{
1720							'Name' => 'Parameters_gap-open',
1721							'Data' => $1
1722							}
1723							);
1724	48					194	$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					221	$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					160	$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					32	$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					24	$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					100	my $c = $1;
1773	48					81	$c =~ s/\,//g;
1774	48					169	$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					145	my $c = $1;
1783	51					132	$c =~ s/\,//g;
1784	51					190	$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					445	my $v = $2;
1796	325					326	chomp($v);
1797	325					1213	$self->element(
1798							{
1799							'Name' => "Statistics_$1",
1800							'Data' => $v
1801							}
1802							);
1803	325	100				797	if ( defined $4 ) {
1804	235					773	$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					74	$self->element(
1818							{
1819							'Name' => 'Statistics_framewindow',
1820							'Data' => $1
1821							}
1822							);
1823	14					60	$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					247	$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					221	$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					256	$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					226	$self->element(
1864							{
1865							'Name' => 'Parameters_expect',
1866							'Data' => $1
1867							}
1868							);
1869	52					208	$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					142	my $d = $1;
1878	50					76	chomp($d);
1879	50					190	$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					5811	$last = $_;
1891							}
1892							} elsif ( $self->in_element('hsp') ) {
1893	9058					18188	$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			16083	$self->{'_reporttype'} \|\|= $DEFAULTREPORTTYPE;
1897	9058					17693	my %data = (
1898							'Query' => '',
1899							'Mid' => '',
1900							'Hit' => ''
1901							);
1902	9058					5840	my $len;
1903	9058		66			28101	for ( my $i = 0 ; defined($_) && $i < 3 ; $i++ ) {
1904							# $self->debug("$i: $_") if $v;
1905	27150	50	66			127040	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					21168	chomp;
1913	27150	100				74979	if (/^((Query\|Sbjct):?\s+(\-?\d+)?\s*)(\S+)\s+(\-?\d+)?/) {
		50
1914	18116					43894	my ( $full, $type, $start, $str, $end ) =
1915							( $1, $2, $3, $4, $5 );
1916
1917	18116	100				21254	if ( $str eq '-' ) {
1918	48	100				83	$i = 3 if $type eq 'Sbjct';
1919							}
1920							else {
1921	18068					24806	$data{$type} = $str;
1922							}
1923	18116					11763	$len = length($full);
1924							$self->{"\_$type"}->{'begin'} = $start
1925	18116	100				37955	unless $self->{"_$type"}->{'begin'};
1926	18116					22762	$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			24692	$self->throw("no data for midline $_")
1933							unless ( defined $_ && defined $len );
1934	9034					14991	$data{'Mid'} = substr( $_, $len );
1935							}
1936	27150					44138	$_ = $self->_readline();
1937							}
1938							$self->characters(
1939							{
1940							'Name' => 'Hsp_qseq',
1941	9058					27473	'Data' => $data{'Query'}
1942							}
1943							);
1944							$self->characters(
1945							{
1946							'Name' => 'Hsp_hseq',
1947	9058					23302	'Data' => $data{'Sbjct'}
1948							}
1949							);
1950							$self->characters(
1951							{
1952							'Name' => 'Hsp_midline',
1953	9058					21235	'Data' => $data{'Mid'}
1954							}
1955							);
1956							}
1957							else {
1958							#$self->debug("blast.pm: unrecognized line $_");
1959							}
1960							}
1961
1962	102					434	$self->debug("blast.pm: End of BlastOutput\n");
1963	102	100				313	if ( $self->{'_seentop'} ) {
1964	96	100				249	$self->within_element('hsp')
1965							&& $self->end_element( { 'Name' => 'Hsp' } );
1966	96	100				501	$self->within_element('hit')
1967							&& $self->end_element( { 'Name' => 'Hit' } );
1968							# cleanup extra hits
1969	96	100				232	$self->_cleanup_hits(\@hit_signifs) if scalar(@hit_signifs);
1970	96	100				229	$self->within_element('iteration')
1971							&& $self->end_element( { 'Name' => 'Iteration' } );
1972	96	100				219	if ($bl2seq_fix) {
1973	6					30	$self->element(
1974							{
1975							'Name' => 'BlastOutput_program',
1976							'Data' => $reporttype
1977							}
1978							);
1979							}
1980	96					303	$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		172	my $self = shift;
1988	96					443	$self->start_element( { 'Name' => 'BlastOutput' } );
1989	96					205	$self->{'_seentop'} = 1;
1990	96					176	$self->{'_result_count'}++;
1991	96					131	$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		6965	my ( $self, $type ) = @_;
2034	8516					7986	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				15581	: ( $self->{'_will_handle_cache'}->{$type} =
2039							$handler->will_handle($type) );
2040
2041	8516	50				12473	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	4189	my ( $self, $data ) = @_;
2056
2057							# we currently don't care about attributes
2058	4258					3897	my $nm = $data->{'Name'};
2059	4258					4580	my $type = $MODEMAP{$nm};
2060	4258	50				6299	if ($type) {
2061	4258					5549	my $handler = $self->_will_handle($type);
2062	4258	50				6444	if ($handler) {
2063	4258					9212	my $func = sprintf( "start_%s", lc $type );
2064	4258					15425	$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					5420	unshift @{ $self->{'_elements'} }, $type;
	4258					6273
2073	4258	100				5556	if ( $type eq 'result' ) {
2074	96					186	$self->{'_values'} = {};
2075	96					258	$self->{'_result'} = undef;
2076							} else {
2077							# cleanup some things
2078	4162	50				7090	if ( defined $self->{'_values'} ) {
2079	4162					3062	foreach my $k (
2080	117000					190996	grep { /^\U$type\-/ }
2081	4162					16454	keys %{ $self->{'_values'} }
2082							)
2083							{
2084	38633					37315	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	31218	my ( $self, $data ) = @_;
2104
2105	43007					33451	my $nm = $data->{'Name'};
2106	43007					25374	my $type;
2107							my $rc;
2108							# cache these (TODO: we should probably cache all cross-report data)
2109	43007	100				51549	if ( $nm eq 'BlastOutput_program' ) {
2110	86	100				470	if ( $self->{'_last_data'} =~ /(t?blast[npx])/i ) {
2111	83					250	$self->{'_reporttype'} = uc $1;
2112							}
2113	86		66			216	$self->{'_reporttype'} \|\|= $DEFAULTREPORTTYPE;
2114							}
2115
2116	43007	100				44691	if ( $nm eq 'BlastOutput_version' ) {
2117	80					142	$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				46452	if ( $nm eq 'Hsp' ) {
2123	1623					2786	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				18725	) if defined $self->{'_last_hspdata'}->{$_};
2130							}
2131	1623					1976	$self->{'_last_hspdata'} = {};
2132							$self->element(
2133							{
2134							'Name' => 'Hsp_query-from',
2135	1623					5356	'Data' => $self->{'_Query'}->{'begin'}
2136							}
2137							);
2138							$self->element(
2139							{
2140							'Name' => 'Hsp_query-to',
2141	1623					4426	'Data' => $self->{'_Query'}->{'end'}
2142							}
2143							);
2144
2145							$self->element(
2146							{
2147							'Name' => 'Hsp_hit-from',
2148	1623					4361	'Data' => $self->{'_Sbjct'}->{'begin'}
2149							}
2150							);
2151							$self->element(
2152							{
2153							'Name' => 'Hsp_hit-to',
2154	1623					4311	'Data' => $self->{'_Sbjct'}->{'end'}
2155							}
2156							);
2157
2158							# } elsif( $nm eq 'Iteration' ) {
2159							# Nothing special needs to be done here.
2160							}
2161	43007	100				79183	if ( $type = $MODEMAP{$nm} ) {
		100
2162	4258					5889	my $handler = $self->_will_handle($type);
2163	4258	50				6361	if ($handler) {
2164	4258					9861	my $func = sprintf( "end_%s", lc $type );
2165	4258					14061	$rc = $handler->$func( $self->{'_reporttype'}, $self->{'_values'} );
2166							}
2167	4258					3285	shift @{ $self->{'_elements'} };
	4258					5418
2168
2169							}
2170							elsif ( $MAPPING{$nm} ) {
2171	38562	100				44885	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					2191	my $key = ( keys %{ $MAPPING{$nm} } )[0];
	2900					6172
2176							$self->{'_values'}->{$key}->{ $MAPPING{$nm}->{$key} } =
2177	2900					7511	$self->{'_last_data'};
2178							}
2179							else {
2180	35662					59938	$self->{'_values'}->{ $MAPPING{$nm} } = $self->{'_last_data'};
2181							}
2182							}
2183							else {
2184							#$self->debug("blast.pm: unknown nm $nm, ignoring\n");
2185							}
2186	43007					33795	$self->{'_last_data'} = ''; # remove read data if we are at
2187							# end of an element
2188	43007	100	100			69822	$self->{'_result'} = $rc if ( defined $type && $type eq 'result' );
2189	43007					34493	$self->{'_seen_hsp_features'} = 0;
2190	43007					51820	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	31266	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					39513	$self->characters($data);
2210	38749					47990	$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	48596	my ( $self, $data ) = @_;
2226	65923	100	100			68202	if ( $self->in_element('hsp')
2227							&& $data->{'Name'} =~ /^Hsp\_(qseq\|hseq\|midline)$/ )
2228							{
2229							$self->{'_last_hspdata'}->{ $data->{'Name'} } .= $data->{'Data'}
2230	32019	100				81936	if defined $data->{'Data'};
2231							}
2232	65923	100	100			225749	return unless ( defined $data->{'Data'} && $data->{'Data'} !~ /^\s+$/ );
2233	65853					97885	$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	2066	my ( $self, $name ) = @_;
2252							return 0
2253							if ( !defined $name && !defined $self->{'_elements'}
2254	2214	100	33			5009	\|\| scalar @{ $self->{'_elements'} } == 0 );
	2214		66			5282
2255	2212					1632	foreach ( @{ $self->{'_elements'} } ) {
	2212					3377
2256	3134	100				4620	if ( $_ eq $name ) {
2257	1949					4602	return 1;
2258							}
2259							}
2260	263					529	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	195460	my ( $self, $name ) = @_;
2279	264029	100				333173	return 0 if !defined $self->{'_elements'}->[0];
2280	263910					990920	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	129	my ($self) = @_;
2296	102					194	$self->{'_lasttype'} = '';
2297	102					248	$self->{'_values'} = {};
2298	102					418	$self->{'_result'} = undef;
2299	102					240	$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	172	my ( $self, @args ) = @_;
2315
2316							#$self->debug("blast.pm: end_document\n");
2317	102					618	return $self->{'_result'};
2318							}
2319
2320							sub write_result {
2321	5			5	1	27	my ( $self, $blast, @args ) = @_;
2322
2323	5	50				20	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					40	$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		17	my ($self, $hits) = @_;
2424	8					15	while ( my $v = shift @{ $hits }) {
	1532					2809
2425	1524	50				2016	next unless defined $v;
2426	1524					3108	$self->start_element( { 'Name' => 'Hit' } );
2427	1524					3209	my $id = $v->[2];
2428	1524					1383	my $desc = $v->[3];
2429	1524					3268	$self->element(
2430							{
2431							'Name' => 'Hit_id',
2432							'Data' => $id
2433							}
2434							);
2435	1524					3766	my ($gi, $acc, $version ) = $self->_get_seq_identifiers($id);
2436	1524					3255	$self->element(
2437							{
2438							'Name' => 'Hit_accession',
2439							'Data' => $acc
2440							}
2441							);
2442	1524	50				2715	if ( defined $v ) {
2443	1524					3169	$self->element(
2444							{
2445							'Name' => 'Hit_signif',
2446							'Data' => $v->[0]
2447							}
2448							);
2449	1524	50				2346	if (exists $self->{'_wublast'}) {
2450	0					0	$self->element(
2451							{
2452							'Name' => 'Hit_score',
2453							'Data' => $v->[1]
2454							}
2455							);
2456							} else {
2457	1524					3177	$self->element(
2458							{
2459							'Name' => 'Hit_bits',
2460							'Data' => $v->[1]
2461							}
2462							);
2463							}
2464
2465							}
2466							$self->element(
2467							{
2468	1524					3358	'Name' => 'Hit_def',
2469							'Data' => $desc
2470							}
2471							);
2472	1524					2948	$self->end_element( { 'Name' => 'Hit' } );
2473							}
2474							}
2475
2476
2477							1;
2478
2479							__END__