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

Annotation of /FigKernelPackages/gjoparseblast.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3