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

Annotation of /FigKernelPackages/gjoparseblast.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3