[Bio] / FigKernelPackages / gjoparseblast.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/gjoparseblast.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : golsen 1.1 package gjoparseblast;
2 :    
3 :     #===============================================================================
4 :     # This is a set of functions for reading blast output from a file, a pipe or
5 :     # an array reference (D = \*STDIN), and providing various perl interfaces.
6 :     #
7 :     # BEWARE: When used with an array reference, the array is consumed as it is
8 :     # read. If you need the original, make a copy.
9 :     #
10 :     # The output data are equivalent to those provided by the rationalize_blast
11 :     # script, but can read directly from blastall, without an intermediate file
12 :     # or a shell. It is possible to gather all the output at once, to gather
13 :     # it a query at a time, to gather it a subject sequence at a time, to
14 :     # gather it an HSP at a time, or to read it record-by-record (as it is
15 :     # provided by the rationalize_blast script). The flexible granulatity of
16 :     # access is intended to faciliate the prossessing of large streams of
17 :     # output without loading it all into memory.
18 :     #
19 :     # The 'self' option enables returning matches with query_id eq subject_id.
20 :     # These are normally discarded (matching the behavior of ratonalize_blast).
21 :     #
22 :     #===============================================================================
23 :     #
24 :     # Structured collection of all blast output:
25 :     #
26 :     # @output = structured_blast_output( $input, $self )
27 :     # \@output = structured_blast_output( $input, $self )
28 :     #
29 :     # Output is clustered heirarchically by query, by subject and by hsp. The
30 :     # highest level is query records:
31 :     #
32 :     # [ qid, qdef, qlen, [ [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
33 :     # [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
34 :     # ...
35 :     # ]
36 :     # ]
37 :     #
38 :     # hsp_data:
39 :     #
40 :     # [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
41 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
42 :     #
43 :     #-------------------------------------------------------------------------------
44 :     #
45 :     # Flattened collection of all blast output, one record per HSP:
46 :     #
47 :     # @hsps = blast_hsp_list( $input, $self )
48 :     # \@hsps = blast_hsp_list( $input, $self )
49 :     #
50 :     # Output records are all of the form:
51 :     #
52 :     # [ qid qdef qlen sid sdef slen scr e_val p_n p_val n_mat n_id n_pos n_gap dir q1 q2 qseq s1 s2 sseq ]
53 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
54 :     #
55 :     #-------------------------------------------------------------------------------
56 :     #
57 :     # Collection of all blast output in the record types of rationalize_blast
58 :     # (Query=, > and HSP):
59 :     #
60 :     # @records = blast_record_list( $input, $self )
61 :     # \@records = blast_record_list( $input, $self )
62 :     #
63 :     # There are 3 record types: 'Query=', '>' and 'HSP', with fields:
64 :     #
65 :     # [ 'Query=' query_id query_def query_len ]
66 :     # 0 1 2 3
67 :     #
68 :     # [ '>' sbjct_id sbjct_def sbjct_len ]
69 :     # 0 1 2 3
70 :     #
71 :     # [ 'HSP' scr exp p_n p_val n_mat n_id n_sim n_gap dir q1 q2 qseq s1 s2 sseq ]
72 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
73 :     #
74 :     #-------------------------------------------------------------------------------
75 :     #
76 :     # Blast output one query at a time:
77 :     #
78 :     # $query_results = next_blast_query( $input, $self )
79 :     #
80 :     # Query record structure is defined above (see structured_blast_output)
81 :     #
82 :     #-------------------------------------------------------------------------------
83 :     #
84 :     # Blast output one subject sequence at a time:
85 :     #
86 :     # $subject_results = next_blast_subject( $input, $self )
87 :     #
88 :     # Output fields are:
89 :     #
90 :     # [ qid, qdef, qlen, sid, sdef, slen, [ hsp_data, hsp_data, ... ] ]
91 :     #
92 :     # hsp_data is defined above (see structured_blast_output)
93 :     #
94 :     #-------------------------------------------------------------------------------
95 :     #
96 :     # Blast output one HSP at a time:
97 :     #
98 :     # $hsp = next_blast_hsp( $input, $self )
99 :     #
100 :     # HSP record fields are defined above (see blast_hsp_list)
101 :     #
102 :     #-------------------------------------------------------------------------------
103 :     #
104 :     # Blast output one record (Query=, > and HSP) at a time:
105 :     #
106 :     # $record = next_blast_record( $input, $self )
107 :     #
108 :     # Record types and fields are defined above (see blast_record_list)
109 :     #
110 :     #===============================================================================
111 :     #
112 :     # The following code fragment would read blast output from STDIN, one
113 :     # HSP at a time, process that HSP and move on. It is the extreme form
114 :     # of "memory-light".
115 :     #
116 :     # my $hsp;
117 :     # while ( defined( $hsp = next_blast_hsp() )
118 :     # {
119 :     # # Process the HSP
120 :     # }
121 :     #
122 :     #-------------------------------------------------------------------------------
123 :     #
124 :     # The following code fragment would launch blastall in a forked process,
125 :     # read blast output one subject at a time, process those HSPs and move on.
126 :     #
127 :     # my( $prog, $query, $db ) = @_;
128 :     #
129 :     # # This tedious approach allows file names with blanks and tabs, not to
130 :     # # mention newlines (sorry, but it was too much to resist):
131 :     #
132 :     # my @command = ('blastall', '-p', $prog, '-d', $db, '-i', $query, '-e', '1e-5');
133 :     #
134 :     # # Of course we could (ab)use perl to make it more (or less) opaque:
135 :     # #
136 :     # # my @command = qw( blastall -p prog -d db -i query -e 1e-5 );
137 :     # # @command[2,6,4] = @_;
138 :     #
139 :     # my $bfh;
140 :     # my $pid = open( $bfh, '-|' );
141 :     # if ( $pid == 0 )
142 :     # {
143 :     # exec( @command );
144 :     # die "'" . join(" ", @command) . "' failed: $!\n";
145 :     # }
146 :     #
147 :     # my $subj_matches;
148 :     # while ( defined( $subj_matches = next_blast_subject( $bfh ) )
149 :     # {
150 :     # my ( $qid, $qdef, $qlen, $sid, $sdef, $slen, $hsps ) = @$subj_matches;
151 :     # foreach ( sort { $b->[5]/$b->[4] <=> $a->[5]/$a->[4] } @$hsps )
152 :     # {
153 :     # # Process the HSPs, sorted by percent identity
154 :     # }
155 :     # }
156 :     #
157 :     #===============================================================================
158 :    
159 :     use strict;
160 :    
161 :     require Exporter;
162 :    
163 :     our @ISA = qw(Exporter);
164 :     our @EXPORT = qw(
165 :     structured_blast_output
166 :     blast_hsp_list
167 :     blast_record_list
168 :    
169 :     next_blast_query
170 :     next_blast_subject
171 :     next_blast_hsp
172 :     next_blast_record
173 :     );
174 :    
175 :    
176 :     #===============================================================================
177 :     # Collect BLAST program output into parsed records. This version returns
178 :     # the entire output in one list, hence care should be taken with multiple
179 :     # query searches. The record types correspond to those returned by the
180 :     # ratonalize_blast script.
181 :     #
182 :     # @records = blast_record_list( $input, $self )
183 :     #
184 :     # or
185 :     #
186 :     # \@records = blast_record_list( $input, $self )
187 :     #
188 :     #
189 :     # $input An input file or pipe handle, or array ref. If it is undef,
190 :     # \*STDIN will be used.
191 :     #
192 :     # $self normally matches of a query to itself are discarded. This can be
193 :     # overriden by setting $self to true.
194 :     #
195 :     # There are 3 output record types: 'Query=', '>' and 'HSP'. Their fields are:
196 :     #
197 :     # [ 'Query=' query_id query_def query_len ]
198 :     # 0 1 2 3
199 :     #
200 :     # [ '>' sbjct_id sbjct_def sbjct_len ]
201 :     # 0 1 2 3
202 :     #
203 :     # [ 'HSP' scr exp p_n p_val n_mat n_id n_sim n_gap dir q1 q2 qseq s1 s2 sseq ]
204 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
205 :     #
206 :     #===============================================================================
207 :    
208 :     sub blast_record_list
209 :     {
210 :     my @out = ();
211 :    
212 :     while ( $_ = next_blast_record( @_ ) ) { push @out, $_; }
213 :    
214 :     wantarray ? @out : \@out;
215 :     }
216 :    
217 :    
218 :     #===============================================================================
219 :     # Collect BLAST program output into parsed records. This version returns
220 :     # the entire output in one list, hence care should be taken with multiple
221 :     # query searches:
222 :     #
223 :     # @hsps = blast_hsp_list( $input, $self )
224 :     #
225 :     # or
226 :     #
227 :     # \@hsps = blast_hsp_list( $input, $self )
228 :     #
229 :     # $input An input file or pipe handle, or array ref. If it is undef,
230 :     # \*STDIN will be used.
231 :     #
232 :     # $self normally matches of a query to itself are discarded. This can be
233 :     # overriden by setting $self to true.
234 :     #
235 :     # There is one record per HSP, and all output records are stand alone, having
236 :     # the query and subject sequence data:
237 :     #
238 :     # qid qdef qlen sid sdef slen scr e_val p_n p_val n_mat n_id n_pos n_gap dir q1 q2 qseq s1 s2 sseq
239 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
240 :     #
241 :     #===============================================================================
242 :    
243 :     sub blast_hsp_list
244 :     {
245 :     my @out = ();
246 :    
247 :     while ( $_ = next_blast_hsp( @_ ) ) { push @out, $_; }
248 :    
249 :     wantarray ? @out : \@out;
250 :     }
251 :    
252 :    
253 :     #===============================================================================
254 :     # Collect BLAST program output into perl structures. This returns
255 :     # the entire output in one list, hence care should be taken with multiple
256 :     # query searches:
257 :     #
258 :     # @output = structured_blast_output( $input, $self )
259 :     #
260 :     # or
261 :     #
262 :     # \@output = structured_blast_output( $input, $self )
263 :     #
264 :     # $input An input file or pipe handle, or array ref. If it is undef,
265 :     # \*STDIN will be used.
266 :     #
267 :     # $self normally matches of a query to itself are discarded. This can be
268 :     # overriden by setting $self to true.
269 :     #
270 :     # Output is clustered heirarchically:
271 :     #
272 :     # ( [ qid, qdef, qlen, [ [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
273 :     # [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
274 :     # ...
275 :     # ]
276 :     # ],
277 :     # [ qid, qdef, qlen, [ [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
278 :     # [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
279 :     # ...
280 :     # ]
281 :     # ],
282 :     # ...
283 :     # )
284 :     #
285 :     # hsp_data = [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
286 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
287 :     #
288 :     # Each query will be reported even if it does not have hits.
289 :     #
290 :     #===============================================================================
291 :    
292 :     sub structured_blast_output
293 :     {
294 :     my ( $subj_list, $hsp_list );
295 :     my @out = ();
296 :    
297 :     while ( $_ = next_blast_record( @_ ) )
298 :     {
299 :     my $type = shift @$_;
300 :    
301 :     if ( $type eq 'Query=' )
302 :     {
303 :     $subj_list = [];
304 :     $hsp_list = undef;
305 :     push @out, [ @$_, $subj_list ];
306 :     }
307 :    
308 :     elsif ( $subj_list and ( $type eq '>' ) )
309 :     {
310 :     $hsp_list = [];
311 :     push @$subj_list, [ @$_, $hsp_list ];
312 :     }
313 :    
314 :     elsif ( $hsp_list && ( $type eq 'HSP' ) )
315 :     {
316 :     push @$hsp_list, $_;
317 :     }
318 :     }
319 :    
320 :     wantarray ? @out : \@out;
321 :     }
322 :    
323 :    
324 :     #===============================================================================
325 :     # Collect BLAST program output into perl structures. This returns
326 :     # the output for one query sequence:
327 :     #
328 :     # $query_results = next_blast_query( $input, $self )
329 :     #
330 :     # $input An input file or pipe handle, or array ref. If it is undef,
331 :     # \*STDIN will be used.
332 :     #
333 :     # $self normally matches of a query to itself are discarded. This can be
334 :     # overriden by setting $self to true.
335 :     #
336 :     # Output structure:
337 :     #
338 :     # [ qid, qdef, qlen, [ [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
339 :     # [ sid, sdef, slen, [ hsp_data, hsp_data, ... ] ],
340 :     # ...
341 :     # ]
342 :     # ]
343 :     #
344 :     # hsp_data = [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
345 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
346 :     #
347 :     # Each query will be reported even if it does not have hits.
348 :     #
349 :     #===============================================================================
350 :    
351 :     {
352 :     my %query_info = ();
353 :    
354 :     sub next_blast_query
355 :     { my ( $input, $self ) = @_;
356 :     $input ||= \*STDIN;
357 :    
358 :     my ( $query_info, $subj_list, $hsp_list );
359 :    
360 :     if ( $query_info = $query_info{ $input } )
361 :     {
362 :     $subj_list = [];
363 :     $query_info{ $input } = undef;
364 :     }
365 :    
366 :     while ( $_ = next_blast_record( @_ ) )
367 :     {
368 :     my $type = shift @$_;
369 :    
370 :     if ( $type eq 'Query=' )
371 :     {
372 :     if ( $subj_list )
373 :     {
374 :     $query_info{ $input } = [ @$_ ];
375 :     return [ @$query_info, $subj_list ];
376 :     }
377 :    
378 :     $query_info = [ @$_ ];
379 :     $subj_list = [];
380 :     $hsp_list = undef;
381 :     }
382 :    
383 :     elsif ( $subj_list and ( $type eq '>' ) )
384 :     {
385 :     $hsp_list = [];
386 :     push @$subj_list, [ @$_, $hsp_list ];
387 :     }
388 :    
389 :     elsif ( $hsp_list && ( $type eq 'HSP' ) )
390 :     {
391 :     push @$hsp_list, $_;
392 :     }
393 :     }
394 :    
395 :     $query_info{ $input } = undef;
396 :     $subj_list ? [ @$query_info, $subj_list ] : undef;
397 :     }
398 :     }
399 :    
400 :    
401 :     #===============================================================================
402 :     # Collect BLAST program output into perl structures. This returns
403 :     # the output for one subject sequence:
404 :     #
405 :     # $subject_results = next_blast_subject( $input, $self )
406 :     #
407 :     # $input An input file or pipe handle, or array ref. If it is undef,
408 :     # \*STDIN will be used.
409 :     #
410 :     # $self normally matches of a query to itself are discarded. This can be
411 :     # overriden by setting $self to true.
412 :     #
413 :     # Output structure:
414 :     #
415 :     # [ qid, qdef, qlen, sid, sdef, slen, [ hsp_data, hsp_data, ... ] ]
416 :     #
417 :     # hsp_data = [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
418 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
419 :     #
420 :     #===============================================================================
421 :    
422 :     {
423 :     my %q_and_s_info = ();
424 :    
425 :     sub next_blast_subject
426 :     { my ( $input, $self ) = @_;
427 :     $input ||= \*STDIN;
428 :    
429 :     my $q_and_s_info = $q_and_s_info{ $input } || [ (undef) x 6 ];
430 :     my $hsp_list = defined( $q_and_s_info->[3] ) ? [] : undef;
431 :    
432 :     while ( $_ = next_blast_record( @_ ) )
433 :     {
434 :     my $type = shift @$_;
435 :    
436 :     if ( $type eq 'Query=' )
437 :     {
438 :     if ( $hsp_list && @$hsp_list )
439 :     {
440 :     $q_and_s_info{ $input } = [ @$_, undef, undef, undef ];
441 :     return [ @$q_and_s_info, $hsp_list ];
442 :     }
443 :    
444 :     @$q_and_s_info = ( @$_, undef, undef, undef );
445 :     $hsp_list = undef;
446 :     }
447 :    
448 :     elsif ( $type eq '>' )
449 :     {
450 :     if ( $hsp_list && @$hsp_list )
451 :     {
452 :     $q_and_s_info{ $input } = [ @$q_and_s_info[0..2], @$_ ];
453 :     return [ @$q_and_s_info, $hsp_list ];
454 :     }
455 :    
456 :     if ( defined( $q_and_s_info->[0] ) )
457 :     {
458 :     @$q_and_s_info[ 3 .. 5 ] = @$_;
459 :     $hsp_list = [];
460 :     }
461 :     }
462 :    
463 :     elsif ( ( $type eq 'HSP' ) && $hsp_list )
464 :     {
465 :     push @$hsp_list, $_;
466 :     }
467 :     }
468 :    
469 :     $q_and_s_info{ $input } = undef;
470 :     $hsp_list ? [ @$q_and_s_info, $hsp_list ] : undef;
471 :     }
472 :     }
473 :    
474 :    
475 :     #===============================================================================
476 :     # Collect BLAST program output into parsed records. This version returns
477 :     # the output one HSP at a time, hence is the memory light version:
478 :     #
479 :     # $hsp = next_blast_hsp( $input, $self )
480 :     #
481 :     # $input An input file or pipe handle, or array ref. If it is undef,
482 :     # \*STDIN will be used.
483 :     #
484 :     # $self normally matches of a query to itself are discarded. This can be
485 :     # overriden by setting $self to true.
486 :     #
487 :     # Output record fields are:
488 :     #
489 :     # qid qdef qlen sid sdef slen scr e_val p_n p_val n_mat n_id n_pos n_gap dir q1 q2 qseq s1 s2 sseq
490 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
491 :     #
492 :     #===============================================================================
493 :    
494 :     {
495 :     my %q_and_s_info = (); # Saving query and subject info between calls
496 :    
497 :     sub next_blast_hsp
498 :     { my ( $input, $self ) = @_;
499 :     $input ||= \*STDIN;
500 :    
501 :     my $q_and_s_info = $q_and_s_info{ $input } || [ (undef) x 6 ];
502 :    
503 :     while ( $_ = next_blast_record( @_ ) )
504 :     {
505 :     my $type = shift @$_;
506 :     if ( $type eq 'Query=' )
507 :     {
508 :     @$q_and_s_info = ( @$_, undef, undef, undef );
509 :     }
510 :    
511 :     elsif ( $type eq '>' )
512 :     {
513 :     @$q_and_s_info[ 3 .. 5 ] = @$_ if defined( $q_and_s_info->[0] );
514 :     }
515 :    
516 :     elsif ( $type eq 'HSP' && defined( $q_and_s_info->[3] ) )
517 :     {
518 :     $q_and_s_info{ $input } = $q_and_s_info;
519 :     return [ @$q_and_s_info, @$_ ];
520 :     }
521 :     }
522 :    
523 :     $q_and_s_info{ $input } = undef;
524 :     return undef;
525 :     }
526 :     }
527 :    
528 :    
529 :     #===============================================================================
530 :     # Collect BLAST program output into parsed records. Each call returns one
531 :     # record. Record types correspond to those of the rationalize_blast script.
532 :     # This can be used to progressively read blast output from a file or pipe,
533 :     # without putting it all in memory.
534 :     #
535 :     # $record = next_blast_record( $input, $self )
536 :     #
537 :     # $input An input file or pipe handle, or array ref. If it is undef,
538 :     # \*STDIN will be used.
539 :     #
540 :     # $self Normally matches of a query to itself are discarded. This can be
541 :     # overriden by setting $self to true.
542 :     #
543 :     # There are 3 output record types: 'Query=', '>' and 'HSP'. Their fields are:
544 :     #
545 :     # [ 'Query=' query_id query_def query_len ]
546 :     # 0 1 2 3
547 :     #
548 :     # [ '>' sbjct_id sbjct_def sbjct_len ]
549 :     # 0 1 2 3
550 :     #
551 :     # [ 'HSP' scr exp p_n p_val n_mat n_id n_sim n_gap dir q1 q2 qseq s1 s2 sseq ]
552 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
553 :     #
554 :     #===============================================================================
555 :    
556 :     {
557 :     my %blast_state = ();
558 :    
559 :     sub next_blast_record
560 :     { my ( $input, $self ) = @_;
561 :     $input ||= \*STDIN;
562 :    
563 :     my ( $qid, $qdef, $qlen, $q1, $q2, $qseq );
564 :     my ( $sid, $sdef, $slen, $s1, $s2, $sseq );
565 :     my ( $s, $e, $n, $p, $ident, $outof, $posit, $ngap, $frame );
566 :     my $q_ok = 0;
567 :     my $s_ok = 0;
568 :     my $hsp_ok = 0;
569 :    
570 :     my $saved_line;
571 :     my $state = $blast_state{ $input };
572 :     ( $saved_line, $q_ok, $qid, $s_ok, $sid ) = @$state if $state;
573 :    
574 :     while ( 1 )
575 :     {
576 :     $_ = defined( $saved_line ) ? $saved_line : nextline( $input );
577 :     $saved_line = undef;
578 :    
579 :     s/\s+$//; # trim trailing spaces (including the newline)
580 :    
581 :     # These all signify the end of an HSP. Report it. ---------------------
582 :    
583 :     if ( $hsp_ok && ( /^Query=/
584 :     || /^>/
585 :     || /^Parameters:/
586 :     || /^ +Score = /
587 :     || /^ +Plus Strand HSPs:/
588 :     || /^ +Minus Strand HSPs:/
589 :     || /^ +Database:/
590 :     || ! defined( $_ )
591 :     )
592 :     )
593 :     {
594 :     if ( $s_ok && $q1 && $qseq && $s1 && $sseq )
595 :     {
596 :     $e =~ s/^e-/1.0e-/i; # Fix missing digits in exponential
597 :     $p =~ s/^e-/1.0e-/i;
598 :    
599 :     $frame ||= (seqdir($q1,$q2) * seqdir($s1,$s2) > 0) ? "+" : "-";
600 :    
601 :     $blast_state{ $input } = [ $_, $q_ok, $qid, $s_ok, $sid ];
602 :     return [ 'HSP', $s, $e, $n, $p,
603 :     $outof, $ident, $posit, $ngap, $frame,
604 :     $q1, $q2, $qseq,
605 :     $s1, $s2, $sseq
606 :     ];
607 :     }
608 :    
609 :     $hsp_ok = 0;
610 :     }
611 :    
612 :     # This is the end condition -------------------------------------------
613 :    
614 :     last if ! defined( $_ );
615 :    
616 :     # Now that reporting is up-to-date, process the line -----------------
617 :    
618 :     # Query sequence description -----------------------------------------
619 :    
620 :     if ( s/^Query=\s+// )
621 :     {
622 :     $q_ok = 0;
623 :     $s_ok = 0;
624 :     $hsp_ok = 0;
625 :    
626 :     ( $qid, $qdef ) = split_id( $_ ); # query description
627 :    
628 :     # Continue reading to query length
629 :    
630 :     while ( defined( $_ = nextline( $input ) ) )
631 :     {
632 :     s/\s+$//; # trailing space
633 :    
634 :     # Query length marks end of description
635 :    
636 :     if ( /^ +\(([1-9][\d,]*) letters.*\)/ )
637 :     {
638 :     $qlen = $1; # grab query length
639 :     $qlen =~ s/,//g; # get rid of commas
640 :     $q_ok = 1;
641 :     $blast_state{ $input } = [ undef, $q_ok, $qid, $s_ok, $sid ];
642 :     return [ 'Query=', $qid, $qdef, $qlen ]; # return query description
643 :     }
644 :    
645 :     # Database before length is an error
646 :    
647 :     elsif ( /^Database:/ ) { last; }
648 :    
649 :     # Otherwise this is a continuation of the query description.
650 :     # Remove leading white space.
651 :    
652 :     s/^\s+//;
653 :     $qdef .= ( $qdef =~ /-$/ ) ? $_ : " $_";
654 :     }
655 :    
656 :     # Failed to get query length?
657 :    
658 :     if ( ! $qlen )
659 :     {
660 :     $_ || ( $_ = defined($_) ? "" : "[EOF]" );
661 :     Print STDERR "Error parsing query definition for sequence length:\n";
662 :     Print STDERR "$qid $qdef\n<<<here>>>\n$_\n\n";
663 :     print STDERR "Flushing this query sequence\n\n";
664 :     next;
665 :     }
666 :     }
667 :    
668 :     # Subject sequence description -----------------------------------------
669 :    
670 :     elsif ( $q_ok && s/^>\s*// )
671 :     {
672 :     $s_ok = 0;
673 :     $hsp_ok = 0;
674 :    
675 :     ( $sid, $sdef ) = split_id( $_ ); # subject description
676 :    
677 :     while ( defined( $_ = nextline( $input ) ) )
678 :     {
679 :     chomp;
680 :    
681 :     # Length marks end of subject sequence description
682 :    
683 :     if ( /^ +Length = ([1-9][\d,]*)/ )
684 :     {
685 :     $slen = $1;
686 :     $slen =~ s/,//g;
687 :     if ( $self || ( $sid ne $qid ) )
688 :     {
689 :     $s_ok = 1;
690 :     $blast_state{ $input } = [ undef, $q_ok, $qid, $s_ok, $sid ];
691 :     return [ '>', $sid, $sdef, $slen ]; # return subject description
692 :     }
693 :     last;
694 :     }
695 :    
696 :     # Multiple spaces marks a continuation
697 :    
698 :     elsif ( s/^ +// )
699 :     {
700 :     $sdef .= ( $sdef =~ /-$/ ) ? $_ : " $_";
701 :     }
702 :    
703 :     # Merged nr entries start with one space
704 :    
705 :     elsif ( s/^ /\001/ ) {
706 :     $sdef .= $_;
707 :     }
708 :    
709 :     # Anything else is an error
710 :    
711 :     else { last; }
712 :     }
713 :    
714 :     if ( ! $slen )
715 :     {
716 :     $_ ||= defined($_) ? "" : "[EOF]";
717 :     print STDERR "Error parsing subject definition for sequence length:\n";
718 :     print STDERR "$sid $sdef\n<<<here>>>\n$_\n\n";
719 :     print STDERR "Flushing this subject sequence\n\n";
720 :     }
721 :     }
722 :    
723 :     # Score marks the start of an HSP description --------------------------
724 :    
725 :     elsif ( $s_ok && /^ Score = +([\d.e+-]+) / )
726 :     {
727 :     $hsp_ok = 0;
728 :    
729 :     $s = $1;
730 :     if ( ! /Expect = +([^ ,]+)/ && ! /Expect[(]\d+[)] = +([^ ,]+)/ )
731 :     {
732 :     print STDERR "Error parsing Score line for Expect:\n";
733 :     print STDERR "Query = $qid; Subject = $sid\n<<<here>>>\n$_\n";
734 :     print STDERR "Flushing this HSP\n\n";
735 :     next;
736 :     }
737 :    
738 :     $e = $1;
739 :     $n = /Expect\((\d+)\)/ ? $1
740 :     : /P\((\d+)\)/ ? $1
741 :     : 1;
742 :     $p = / = +(\S+)$/ ? $1 : $e;
743 :    
744 :     if ( ! defined( $_ = nextline( $input ) ) )
745 :     {
746 :     print STDERR "End-of-file while looking for Identities line:\n";
747 :     print STDERR "Query = $qid; Subject = $sid\n<<<here>>>\n";
748 :     print STDERR "Flushing this HSP\n\n";
749 :     next;
750 :     }
751 :     chomp;
752 :    
753 :     if ( ! /^ Identities = +(\d+)\/(\d+)/ )
754 :     {
755 :     print STDERR "Error parsing Identities line:\n";
756 :     print STDERR "Query = $qid; Subject = $sid\n<<<here>>>\n$_\n";
757 :     print STDERR "Flushing this HSP\n\n";
758 :     next;
759 :     }
760 :     $ident = $1;
761 :     $outof = $2;
762 :    
763 :     $posit = /Positives = +(\d+)/ ? $1 : $ident;
764 :     $ngap = /Gaps = +(\d+)/ ? $1 : 0;
765 :    
766 :     $q1 = $s1 = undef;
767 :     $qseq = $sseq = "";
768 :     $hsp_ok = 1;
769 :     }
770 :    
771 :     # Frame of a translated blast -----------------------------------------
772 :    
773 :     elsif ( $hsp_ok && /^ Frame = +(\S+)/ )
774 :     {
775 :     $frame = $1;
776 :     }
777 :    
778 :     # Query sequence data -------------------------------------------------
779 :    
780 :     elsif ( $hsp_ok && s/^Query: +// )
781 :     {
782 :     my ($t1, $t2, $t3) = /(\d+)\s*([^\d\s]*)\s*(\d+)/;
783 :    
784 :     # First fragment of alignment?
785 :    
786 :     if ( ! $q1 )
787 :     {
788 :     $q1 = $t1;
789 :     }
790 :    
791 :     # Additional fragment of alignment:
792 :    
793 :     elsif ( ( $t1 != $q2 + 1 ) && ( $t1 != $q2 - 1 ) )
794 :     {
795 :     print STDERR "Warning: Query position $t1 follows $q2\n";
796 :     print STDERR "Query = $qid; Subject = $sid\n";
797 :     print STDERR "$_\nFlushing this HSP\n\n";
798 :     $hsp_ok = 0;
799 :     next;
800 :     }
801 :    
802 :     $qseq .= $t2; # Append sequence
803 :     $q2 = $t3; # New last residue number
804 :    
805 :     # Flush the alignment match symbol line
806 :    
807 :     if ( ! defined( $_ = nextline( $input ) ) )
808 :     {
809 :     print STDERR "End-of-file while reading alignment:\n";
810 :     print STDERR "Query = $qid; Subject = $sid\n";
811 :     print STDERR "Flushing this HSP\n\n";
812 :     $hsp_ok = 0;
813 :     }
814 :     }
815 :    
816 :     # Subject sequence data -----------------------------------------------
817 :    
818 :     elsif ( $hsp_ok && s/^Sbjct: +// )
819 :     {
820 :     my ( $t1, $t2, $t3 ) = /(\d+)\s*([^\d\s]*)\s*(\d+)/;
821 :    
822 :     # First fragment of alignment?
823 :    
824 :     if ( ! $s1 )
825 :     {
826 :     $s1 = $t1;
827 :     }
828 :    
829 :     # Additional fragment of alignment:
830 :    
831 :     elsif ( ( $t1 != $s2 + 1 ) && ( $t1 != $s2 - 1 ) )
832 :     {
833 :     print STDERR "Warning: Subject position $t1 follows $q2\n";
834 :     print STDERR "Query = $qid; Subject = $sid\n";
835 :     print STDERR "$_\nFlushing this HSP\n\n";
836 :     $hsp_ok = 0;
837 :     next;
838 :     }
839 :    
840 :     $sseq .= $t2; # Append sequence
841 :     $s2 = $t3; # New last residue number
842 :    
843 :     # Flush the blank line
844 :    
845 :     nextline( $input );
846 :     }
847 :    
848 :     # Back up to the top to read another line -----------------------------
849 :     }
850 :    
851 :     # Breaking of the loop when an undefined input line is encountered,
852 :     # we are done.
853 :    
854 :     $blast_state{ $input } = undef;
855 :     return undef;
856 :     }
857 :     } # End of bare block with state information
858 :    
859 :    
860 :     #===============================================================================
861 :     # Useful functions:
862 :     #===============================================================================
863 :    
864 :     # This consumes the array. It should use keep a counter, but there are issues
865 :     # of getting the same array twice, confusing the issue.
866 :    
867 :     sub nextline
868 :     {
869 :     my $input = shift;
870 :     ref( $input ) eq "GLOB" ? <$input> :
871 :     ref( $input ) eq "ARRAY" ? shift @$input :
872 :     undef
873 :     }
874 :    
875 :    
876 :     sub split_id
877 :     {
878 :     ( $_[0] =~ m/(\S+)(\s+(\S.*))?$/ ) ? ( $1, defined($3) ? $3 : "" ) : ();
879 :     }
880 :    
881 :    
882 :     sub seqdir { ( $_[0] <= $_[1] ) ? 1 : -1 }
883 :    
884 :    
885 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3