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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package gjoalignandtree;
2 :    
3 : golsen 1.3 # This is a SAS component.
4 :    
5 :     #-------------------------------------------------------------------------------
6 :     #
7 :     # Order a set of sequences:
8 :     #
9 :     # @seqs = order_sequences_by_length( \@seqs, \%opts )
10 :     # \@seqs = order_sequences_by_length( \@seqs, \%opts )
11 :     #
12 :     # Write representative seqeunce groups to a file
13 :     #
14 :     # write_representative_groups( \%rep, $file )
15 : overbeek 1.1 #
16 : golsen 1.3 # @align = trim_align_to_median_ends( \@align, \%opts )
17 :     # \@align = trim_align_to_median_ends( \@align, \%opts )
18 :     #
19 :     # print_alignment_as_pseudoclustal( $align, \%opts )
20 :     #
21 :     # $prof_rep = representative_for_profile( $align )
22 :     #
23 :     # \@trimmed_seq = extract_with_psiblast( $db, $profile, \%opts )
24 :     # ( \@trimmed_seq, \%report ) = extract_with_psiblast( $db, $profile, \%opts )
25 :     #
26 :     # $structured_blast = blastpgp( $dbfile, $profilefile, \%options )
27 :     # $structured_blast = blastpgp( \@dbseq, $profilefile, \%options )
28 :     # $structured_blast = blastpgp( $dbfile, \@profileseq, \%options )
29 :     # $structured_blast = blastpgp( \@dbseq, \@profileseq, \%options )
30 :     #
31 :     # print_blast_as_records( $file, \@queries )
32 :     # print_blast_as_records( \*FH, \@queries )
33 :     # print_blast_as_records( \@queries ) # STDOUT
34 :     #
35 :     # \@queries = read_blast_from_records( $file )
36 :     # \@queries = read_blast_from_records( \*FH )
37 :     # \@queries = read_blast_from_records( ) # STDIN
38 :     #
39 :     # $exists = verify_db( $db, $type );
40 :     #
41 :     #-------------------------------------------------------------------------------
42 : overbeek 1.1
43 :     use strict;
44 : golsen 1.3 use SeedAware;
45 : overbeek 1.1 use gjoseqlib;
46 : golsen 1.3 use gjoparseblast;
47 : overbeek 1.1
48 :     #-------------------------------------------------------------------------------
49 :     # Order a set of sequences:
50 :     #
51 :     # @seqs = order_sequences_by_length( \@seqs, \%opts )
52 :     # \@seqs = order_sequences_by_length( \@seqs, \%opts )
53 :     #
54 :     # Options:
55 :     #
56 : golsen 1.3 # order => $key # increasing (D), decreasing, median_outward,
57 :     # # closest_to_median, closest_to_length
58 :     # length => $prefered_length
59 : overbeek 1.1 #
60 :     #-------------------------------------------------------------------------------
61 :    
62 :     sub order_sequences_by_length
63 :     {
64 :     my ( $seqs, $opts ) = @_;
65 :     $seqs && ref $seqs eq 'ARRAY' && @$seqs
66 :     or print STDERR "order_sequences_by_length called with invalid sequence list.\n"
67 :     and return wantarray ? () : [];
68 :    
69 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
70 :     my $order = $opts->{ order } || 'increasing';
71 :    
72 :     my @seqs = ();
73 :     if ( $order =~ /^dec/i )
74 :     {
75 :     $opts->{ order } = 'decreaseing';
76 :     @seqs = sort { length( $b->[2] ) <=> length( $a->[2] ) } @$seqs;
77 :     }
78 :     elsif ( $order =~ /^med/i )
79 :     {
80 :     $opts->{ order } = 'median_outward';
81 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
82 :     while ( $_ = shift @seqs0 )
83 :     {
84 :     push @seqs, $_;
85 :     push @seqs, pop @seqs0 if @seqs0;
86 :     }
87 :     @seqs = reverse @seqs;
88 :     }
89 :     elsif ( $order =~ /^close/i )
90 :     {
91 :     my $pref_len;
92 :     if ( defined $opts->{ length } ) # caller supplies preferred length?
93 :     {
94 :     $opts->{ order } = 'closest_to_length';
95 :     $pref_len = $opts->{ length } + 0;
96 :     }
97 :     else # preferred length is median
98 :     {
99 :     $opts->{ order } = 'closest_to_median';
100 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
101 :     my $i = int( @seqs0 / 2 );
102 :     $pref_len = ( @seqs0 % 2 ) ? length( $seqs0[$i]->[2] )
103 :     : ( length( $seqs0[$i-1]->[2] ) + length( $seqs0[$i]->[2] ) ) / 2;
104 :     }
105 :    
106 :     @seqs = map { $_->[0] }
107 :     sort { $a->[1] <=> $b->[1] } # closest to preferred first
108 :     map { [ $_, abs( length( $_->[2] ) - $pref_len ) ] } # dist from pref?
109 :     @$seqs;
110 :     }
111 :     else
112 :     {
113 :     $opts->{ order } = 'increasing';
114 :     @seqs = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
115 :     }
116 :    
117 :     wantarray ? @seqs : \@seqs;
118 :     }
119 :    
120 :    
121 :     #-------------------------------------------------------------------------------
122 :     # Write representative seqeunce groups to a file
123 :     #
124 :     # write_representative_groups( \%rep, $file )
125 :     #
126 :     #-------------------------------------------------------------------------------
127 :    
128 :     sub write_representative_groups
129 :     {
130 :     my ( $repH, $file ) = @_;
131 :     $repH && (ref $repH eq 'HASH') && keys %$repH
132 :     or print STDERR "write_representative_groups called with invalid rep hash.\n"
133 :     and return;
134 :    
135 :     my ( $fh, $close ) = &output_file_handle( $file );
136 :    
137 :     # Order keys from largest to smallest set, then alphabetical
138 :    
139 :     my @keys = sort { @{ $repH->{$b} } <=> @{ $repH->{$a} } || lc $a cmp lc $b } keys %$repH;
140 :    
141 :     foreach ( @keys ) { print $fh join( "\t", ( $_, @{ $repH->{ $_ } } ) ), "\n" }
142 :    
143 :     close $fh if $close;
144 :     }
145 :    
146 :    
147 :     #-------------------------------------------------------------------------------
148 :     #
149 :     # @align = trim_align_to_median_ends( \@align, \%opts )
150 :     # \@align = trim_align_to_median_ends( \@align, \%opts )
151 :     #
152 :     # Options:
153 :     #
154 :     # begin => bool # Trim start (specifically)
155 :     # end => bool # Trim end (specifically)
156 :     #
157 :     #-------------------------------------------------------------------------------
158 :     sub trim_align_to_median_ends
159 :     {
160 :     my ( $align, $opts ) = @_;
161 :     $align && ref $align eq 'ARRAY' && @$align
162 :     or print STDERR "trim_align_to_median_ends called with invalid sequence list.\n"
163 :     and return wantarray ? () : [];
164 :    
165 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
166 :     my $tr_beg = $opts->{ begin } || $opts->{ beg } || $opts->{ start } ? 1 : 0;
167 :     my $tr_end = $opts->{ end } ? 1 : 0;
168 :     $tr_beg = $tr_end = 1 if ! ( $tr_beg || $tr_end );
169 :    
170 :     my( @pos1, @pos2, $i );
171 :     foreach my $seq ( @$align )
172 :     {
173 :     my( $b, $e ) = $seq->[2] =~ /^(-*).*[^-](-*)$/;
174 :     push @pos1, length( $b || '' );
175 :     push @pos2, length( $e || '' );
176 :     }
177 :    
178 :     @pos1 = sort { $a <=> $b } @pos1;
179 :     @pos2 = sort { $a <=> $b } @pos2;
180 :    
181 :     my $pos1 = $tr_beg ? $pos1[ int( @pos1/2 ) ] : 0;
182 :     my $pos2 = $tr_end ? $pos2[ int( @pos2/2 ) ] : 0;
183 :    
184 :     my $ori_len = length( $align->[0]->[2] );
185 :     my $new_len = $ori_len - ( $pos1 + $pos2 );
186 :     my @align2 = map { [ @$_[0,1], substr( $_->[2], $pos1, $new_len ) ] }
187 :     @$align;
188 :    
189 :     wantarray ? @align2 : \@align2;
190 :     }
191 :    
192 :    
193 :     #-------------------------------------------------------------------------------
194 :     #
195 : golsen 1.3 # print_alignment_as_pseudoclustal( $align, \%opts )
196 :     #
197 : overbeek 1.1 # Options:
198 :     #
199 :     # file => $filename # supply a file name to open and write
200 :     # file => \*FH # supply an open file handle (D = STDOUT)
201 :     # line => $linelen # residues per line (D = 60)
202 :     # lower => $bool # all lower case sequence
203 :     # upper => $bool # all upper case sequence
204 :     #
205 :     #-------------------------------------------------------------------------------
206 :    
207 :     sub print_alignment_as_pseudoclustal
208 :     {
209 :     my ( $align, $opts ) = @_;
210 :     $align && ref $align eq 'ARRAY' && @$align
211 :     or print STDERR "print_alignment_as_pseudoclustal called with invalid sequence list.\n"
212 :     and return wantarray ? () : [];
213 :    
214 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
215 :     my $line_len = $opts->{ line } || 60;
216 :     my $case = $opts->{ upper } ? 1 : $opts->{ lower } ? -1 : 0;
217 :    
218 :     my ( $fh, $close ) = &output_file_handle( $opts->{ file } );
219 :    
220 :     my $namelen = 0;
221 :     foreach ( @$align ) { $namelen = length $_->[0] if $namelen < length $_->[0] }
222 :     my $fmt = "%-${namelen}s %s\n";
223 :    
224 :     my $id;
225 :     my @lines = map { $id = $_->[0]; [ map { sprintf $fmt, $id, $_ }
226 : golsen 1.3 map { $case < 0 ? lc $_ : $case > 0 ? uc $_ : $_ } # map sequence only
227 : overbeek 1.1 $_->[2] =~ m/.{1,$line_len}/g
228 :     ] }
229 :    
230 :     @$align;
231 :    
232 :     my $ngroup = @{ $lines[0] };
233 :     for ( my $i = 0; $i < $ngroup; $i++ )
234 :     {
235 :     foreach ( @lines ) { print $fh $_->[$i] if $_->[$i] }
236 :     print $fh "\n";
237 :     }
238 :    
239 :     close $fh if $close;
240 :     }
241 :    
242 :    
243 :     #-------------------------------------------------------------------------------
244 : golsen 1.3 #
245 :     # @seqs = read_pseudoclustal( ) # D = STDIN
246 :     # \@seqs = read_pseudoclustal( ) # D = STDIN
247 :     # @seqs = read_pseudoclustal( $file_name )
248 :     # \@seqs = read_pseudoclustal( $file_name )
249 :     # @seqs = read_pseudoclustal( \*FH )
250 :     # \@seqs = read_pseudoclustal( \*FH )
251 :     #
252 :     #-------------------------------------------------------------------------------
253 :    
254 :     sub read_pseudoclustal
255 :     {
256 :     my ( $file ) = @_;
257 :     my ( $fh, $close ) = input_file_handle( $file );
258 :     my %seq;
259 :     my @ids;
260 :     while ( <$fh> )
261 :     {
262 :     chomp;
263 :     my ( $id, $data ) = /^(\S+)\s+(\S.*)$/;
264 :     if ( defined $id && defined $data )
265 :     {
266 :     push @ids, $id if ! $seq{ $id };
267 :     $data =~ s/\s+//g;
268 :     push @{ $seq{ $id } }, $data;
269 :     }
270 :     }
271 :     close $fh if $close;
272 :    
273 :     my @seq = map { [ $_, '', join( '', @{ $seq{ $_ } } ) ] } @ids;
274 :     wantarray ? @seq : \@seq;
275 :     }
276 :    
277 :    
278 :     #-------------------------------------------------------------------------------
279 : overbeek 1.1 # The profile 'query' sequence:
280 :     # 1. No terminal gaps, if possible
281 :     # 1.1 Otherwise, no gaps at start
282 :     # 2. Longest sequence passing above
283 :     #
284 :     # $prof_rep = representative_for_profile( $align )
285 :     #-------------------------------------------------------------------------------
286 :     sub representative_for_profile
287 :     {
288 :     my ( $align ) = @_;
289 :     $align && ref $align eq 'ARRAY' && @$align
290 :     or die "representative_for_profile called with invalid sequence list.\n";
291 :    
292 :     my @cand;
293 :     @cand = grep { $_->[2] =~ /^[^-].*[^-]$/ } @$align; # No end gaps
294 :     @cand = grep { $_->[2] =~ /^[^-]/ } @$align if ! @cand; # No start gaps
295 :    
296 : golsen 1.3 my ( $rep ) = map { $_->[0] } # sequence entry
297 :     sort { $b->[1] <=> $a->[1] } # max nongaps
298 :     map { [ $_, $_->[2] =~ tr/-//c ] } # count nongaps
299 : overbeek 1.1 @cand;
300 :    
301 :     $rep = [ @$rep ]; # Make a copy
302 :     $rep->[2] =~ s/-+//g; # Compress sequence
303 :    
304 :     return $rep;
305 :     }
306 :    
307 :    
308 :     #-------------------------------------------------------------------------------
309 :     #
310 : golsen 1.3 # \@seq = extract_with_psiblast( $db, $profile, \%opts )
311 :     # ( \@seq, \%report ) = extract_with_psiblast( $db, $profile, \%opts )
312 :     #
313 :     # $db can be $db_file_name or \@db_seqs
314 :     # $profile can be $profile_file_name or \@profile_seqs
315 :     #
316 :     # If supplied as file, profile is pseudoclustal
317 :     #
318 :     # Report records:
319 :     #
320 :     # [ $sid, $slen, $status, $frac_id, $frac_pos,
321 :     # $q_uncov_n_term, $q_uncov_c_term,
322 :     # $s_uncov_n_term, $s_uncov_c_term ]
323 :     #
324 :     # Options:
325 :     #
326 :     # e_value => $max_e_value # maximum blastpgp E-value (D = 0.01)
327 :     # max_e_val => $max_e_value # maximum blastpgp E-value (D = 0.01)
328 :     # max_q_uncov => $aa_per_end # maximum unmatched query (D = 20)
329 :     # max_q_uncov_c => $aa_per_end # maximum unmatched query, c-term (D = 20)
330 :     # max_q_uncov_n => $aa_per_end # maximum unmatched query, n-term (D = 20)
331 :     # min_ident => $frac_ident # minimum fraction identity (D = 0.15)
332 :     # min_positive => $frac_positive # minimum fraction positive scoring (D = 0.20)
333 :     # nresult => $max_seq # maximim matching sequences (D = 5000)
334 :     # query => $q_file_name # query sequence file (D = most complete)
335 :     # query => \@q_seq_entry # query sequence (D = most complete)
336 :     # stderr => $file # blastpgp stderr (D = /dev/stderr)
337 :     #
338 :     # If supplied, query must be identical to a profile sequence, but no gaps.
339 :     #
340 :     #-------------------------------------------------------------------------------
341 :    
342 :     sub extract_with_psiblast
343 :     {
344 :     my ( $profile, $seqs, $opts ) = @_;
345 :    
346 :     $opts && ref $opts eq 'HASH' or $opts = {};
347 :    
348 :     $opts->{ e_value } ||= $opts->{ max_e_val } ||= 0.01;
349 :     $opts->{ nresult } ||= 5000;
350 :     my $max_q_uncov_c = $opts->{ max_q_uncov_c } || $opts->{ max_q_uncov } || 20;
351 :     my $max_q_uncov_n = $opts->{ max_q_uncov_n } || $opts->{ max_q_uncov } || 20;
352 :     my $min_ident = $opts->{ min_ident } || 0.15;
353 :     my $min_pos = $opts->{ min_positive } || 0.20;
354 :    
355 :     my $blast = blastpgp( $profile, $seqs, $opts );
356 :     $blast && @$blast or return ();
357 :    
358 :     my ( $qid, $qdef, $qlen, $qhits ) = @{ $blast->[0] };
359 :    
360 :     my @trimmed;
361 :     my %report;
362 :    
363 :     foreach my $sdata ( @$qhits )
364 :     {
365 :     my ( $sid, $sdef, $slen, $hsps ) = @$sdata;
366 :    
367 :     if ( one_real_hsp( $hsps ) )
368 :     {
369 :     # [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
370 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
371 :    
372 :     my $hsp0 = $hsps->[0];
373 :     my ( $scr, $nmat, $nid, $npos, $q1, $q2, $qseq, $s1, $s2, $sseq ) = ( @$hsp0 )[ 0, 4, 5, 6, 9, 10, 11, 12, 13, 14 ];
374 :    
375 :     my $status;
376 :     if ( $q1-1 > $max_q_uncov_n ) { $status = 'missing start' }
377 :     elsif ( $qlen-$q2 > $max_q_uncov_c ) { $status = 'missing end' }
378 :     elsif ( $nid / $nmat < $min_ident ) { $status = 'low identity' }
379 :     elsif ( $npos / $nmat < $min_pos ) { $status = 'low positives' }
380 :     else
381 :     {
382 :     $sseq =~ s/-+//g;
383 :     $sdef .= " ($s1-$s2/$slen)" if ( ( $s1 > 1 ) || ( $s2 < $slen ) );
384 :     push @trimmed, [ $sid, $sdef, $sseq ];
385 :     $status = 'included';
386 :     }
387 :    
388 :     $report{ $sid } = [ $sid, $slen, $status, $nid/$nmat, $npos/$nmat, $q1-1, $qlen-$q2, $s1-1, $slen-$s2 ];
389 :     }
390 :     }
391 :    
392 :     wantarray ? ( \@trimmed, \%report ) : \@trimmed;
393 :     }
394 :    
395 :    
396 :     #-------------------------------------------------------------------------------
397 :     #
398 :     # Allow fragmentary matches inside of the highest-scoring hsp:
399 :     #
400 :     #-------------------------------------------------------------------------------
401 :    
402 :     sub one_real_hsp
403 :     {
404 :     my ( $hsps ) = @_;
405 :     return 0 if ! ( $hsps && ( ref( $hsps ) eq 'ARRAY' ) && @$hsps );
406 :     return 1 if @$hsps == 1;
407 :    
408 :     my ( $q1_0, $q2_0 ) = ( @{ $hsps->[0] } )[9, 10];
409 :     for ( my $i = 1; $i < @$hsps; $i++ )
410 :     {
411 :     my ($q1, $q2) = ( @{ $hsps->[$i] } )[9, 10];
412 :     return 0 if $q1 < $q1_0 || $q2 > $q2_0;
413 :     }
414 :    
415 :     return 1;
416 :     }
417 :    
418 :    
419 :     #-------------------------------------------------------------------------------
420 :     #
421 : overbeek 1.1 # $structured_blast = blastpgp( $db, $profile, $options )
422 :     #
423 :     # Required:
424 :     #
425 : golsen 1.3 # $db_file or \@db_seq
426 :     # $profile_file or \@profile_seq
427 :     #
428 :     # Options:
429 :     #
430 :     # e_value => $max_e_value # D = 0.01
431 :     # max_e_val => $max_e_value # D = 0.01
432 :     # nresult => $max_seq # D = 5000
433 :     # query => $q_file_name # most complete
434 :     # query => \@q_seq_entry # most complete
435 : overbeek 1.1 #
436 :     #-------------------------------------------------------------------------------
437 :    
438 :     sub blastpgp
439 :     {
440 :     my ( $db, $profile, $opts ) = @_;
441 :     $opts ||= {};
442 : golsen 1.3 my $tmp = SeedAware::location_of_tmp();
443 : overbeek 1.1
444 :     my ( $dbfile, $rm_db );
445 :     if ( defined $db && ref $db )
446 :     {
447 :     ref $db eq 'ARRAY'
448 :     && @$db
449 :     || print STDERR "blastpgp requires one or more database sequences.\n"
450 :     && return undef;
451 : golsen 1.3 $dbfile = "$tmp/blastpgp_db_$$";
452 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $dbfile, $db );
453 :     $rm_db = 1;
454 :     }
455 :     elsif ( defined $db && -f $db )
456 :     {
457 :     $dbfile = $db;
458 :     }
459 :     else
460 :     {
461 :     die "blastpgp requires database.";
462 :     }
463 : golsen 1.3 verify_db( $dbfile, 'P' ); # protein
464 : overbeek 1.1
465 :     my ( $proffile, $rm_profile );
466 :     if ( defined $profile && ref $profile )
467 :     {
468 :     ref $profile eq 'ARRAY'
469 :     && @$profile
470 :     || print STDERR "blastpgp requires one or more profile sequences.\n"
471 :     && return undef;
472 : golsen 1.3 $proffile = "$tmp/blastpgp_profile_$$";
473 : overbeek 1.1 &print_alignment_as_pseudoclustal( $profile, { file => $proffile, upper => 1 } );
474 :     $rm_profile = 1;
475 :     }
476 :     elsif ( defined $profile && -f $profile )
477 :     {
478 :     $proffile = $profile;
479 :     }
480 :     else
481 :     {
482 :     die "blastpgp requires profile.";
483 :     }
484 :    
485 :     my ( $qfile, $rm_query );
486 :     my $query = $opts->{ query };
487 :     if ( defined $query && ref $query )
488 :     {
489 : golsen 1.3 ref $query eq 'ARRAY' && @$query == 3
490 :     or print STDERR "blastpgp invalid query sequence.\n"
491 :     and return undef;
492 :    
493 :     $qfile = "$tmp/blastpgp_query_$$";
494 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
495 :     $rm_query = 1;
496 :     }
497 :     elsif ( defined $query && -f $query )
498 :     {
499 :     $qfile = $query;
500 :     }
501 :     elsif ( $profile ) # Build it from profile
502 :     {
503 :     $query = &representative_for_profile( $profile );
504 : golsen 1.3 $qfile = "$tmp/blastpgp_query_$$";
505 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
506 :     $rm_query = 1;
507 :     }
508 :     else
509 : golsen 1.3 {
510 : overbeek 1.1 die "blastpgp requires database.";
511 :     }
512 :    
513 : golsen 1.3 my $nkeep = $opts->{ nresult } || $opts->{ n_result } || 5000;
514 :     my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || 0.01;
515 :     my @cmd = ( SeedAware::executable_for( 'blastpgp' ),
516 :     '-B', $proffile,
517 :     '-j', 1,
518 :     '-d', $dbfile,
519 :     '-b', $nkeep,
520 :     '-v', $nkeep,
521 :     '-i', $qfile,
522 :     '-e', $e_val
523 :     );
524 :    
525 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts->{stderr} ? { stderr => $opts->{stderr} } : () )
526 :     or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
527 :     and return undef;
528 :     my $out = &gjoparseblast::structured_blast_output( $blastfh, 1 ); # selfmatches okay
529 :     close $blastfh;
530 : overbeek 1.1
531 :     if ( $rm_db )
532 :     {
533 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
534 :     unlink @files if @files;
535 :     }
536 :     unlink $proffile if $rm_profile;
537 :     unlink $qfile if $rm_query;
538 :    
539 :     return $out;
540 :     }
541 :    
542 :    
543 :     #-------------------------------------------------------------------------------
544 : golsen 1.3 # Get an input file handle, and boolean on whether to close or not:
545 :     #
546 :     # ( \*FH, $close ) = input_file_handle( $filename );
547 :     # ( \*FH, $close ) = input_file_handle( \*FH );
548 :     # ( \*FH, $close ) = input_file_handle( ); # D = STDIN
549 :     #
550 :     #-------------------------------------------------------------------------------
551 :    
552 :     sub input_file_handle
553 :     {
554 :     my ( $file, $umask ) = @_;
555 :    
556 :     my ( $fh, $close );
557 :     if ( defined $file )
558 :     {
559 :     if ( ref $file eq 'GLOB' )
560 :     {
561 :     $fh = $file;
562 :     $close = 0;
563 :     }
564 :     elsif ( -f $file )
565 :     {
566 :     open( $fh, "<$file") || die "input_file_handle could not open '$file'.\n";
567 :     $close = 1;
568 :     }
569 :     else
570 :     {
571 :     die "input_file_handle could not find file '$file'.\n";
572 :     }
573 :     }
574 :     else
575 :     {
576 :     $fh = \*STDIN;
577 :     $close = 0;
578 :     }
579 :    
580 :     return ( $fh, $close );
581 :     }
582 :    
583 :    
584 :     #-------------------------------------------------------------------------------
585 : overbeek 1.1 # Get an output file handle, and boolean on whether to close or not:
586 :     #
587 :     # ( \*FH, $close ) = output_file_handle( $filename );
588 :     # ( \*FH, $close ) = output_file_handle( \*FH );
589 :     # ( \*FH, $close ) = output_file_handle( ); # D = STDOUT
590 :     #
591 :     #-------------------------------------------------------------------------------
592 :    
593 :     sub output_file_handle
594 :     {
595 :     my ( $file, $umask ) = @_;
596 :    
597 :     my ( $fh, $close );
598 :     if ( defined $file )
599 :     {
600 :     if ( ref $file eq 'GLOB' )
601 :     {
602 :     $fh = $file;
603 :     $close = 0;
604 :     }
605 :     else
606 :     {
607 :     open( $fh, ">$file") || die "output_file_handle could not open '$file'.\n";
608 :     chmod 0664, $file; # Seems to work on open file!
609 :     $close = 1;
610 :     }
611 :     }
612 :     else
613 :     {
614 :     $fh = \*STDOUT;
615 :     $close = 0;
616 :     }
617 :    
618 :     return ( $fh, $close );
619 :     }
620 :    
621 :    
622 :     #-------------------------------------------------------------------------------
623 :     #
624 :     # print_blast_as_records( $file, \@queries )
625 :     # print_blast_as_records( \*FH, \@queries )
626 :     # print_blast_as_records( \@queries ) # STDOUT
627 :     #
628 :     # \@queries = read_blast_from_records( $file )
629 :     # \@queries = read_blast_from_records( \*FH )
630 :     # \@queries = read_blast_from_records( ) # STDIN
631 :     #
632 :     # Three output record types:
633 :     #
634 :     # Query \t qid \t qdef \t dlen
635 :     # > \t sid \t sdef \t slen
636 :     # HSP \t ...
637 :     #
638 :     #-------------------------------------------------------------------------------
639 :     sub print_blast_as_records
640 :     {
641 :     my $file = ( $_[0] && ! ( ref $_[0] eq 'ARRAY' ) ? shift : undef );
642 :     my ( $fh, $close ) = output_file_handle( $file );
643 :    
644 :     my $queries = shift;
645 :     $queries && ref $queries eq 'ARRAY'
646 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
647 :     and return;
648 :    
649 :     foreach my $qdata ( @$queries )
650 :     {
651 :     $qdata && ref $qdata eq 'ARRAY' && @$qdata == 4
652 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
653 :     and return;
654 :    
655 :     my $subjcts = pop @$qdata;
656 :     $subjcts && ref $subjcts eq 'ARRAY'
657 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
658 :     and return;
659 :    
660 :     print $fh join( "\t", 'Query', @$qdata ), "\n";
661 :     foreach my $sdata ( @$subjcts )
662 :     {
663 :     $sdata && ref $sdata eq 'ARRAY' && @$sdata == 4
664 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
665 :     and return;
666 :    
667 :     my $hsps = pop @$sdata;
668 :     $hsps && ref $hsps eq 'ARRAY' && @$hsps
669 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
670 :     and return;
671 :    
672 :     print $fh join( "\t", '>', @$sdata ), "\n";
673 :    
674 :     foreach my $hsp ( @$hsps )
675 :     {
676 :     $hsp && ref $hsp eq 'ARRAY' && @$hsp > 4
677 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
678 :     and return;
679 :     print $fh join( "\t", 'HSP', @$hsp ), "\n";
680 :     }
681 :     }
682 :     }
683 :    
684 :     close $fh if $close;
685 :     }
686 :    
687 :    
688 :     sub read_blast_from_records
689 :     {
690 :     my ( $file ) = @_;
691 : golsen 1.3
692 :     my ( $fh, $close ) = input_file_handle( $file );
693 :     $fh or print STDERR "read_blast_from_records could not open input file.\n"
694 : overbeek 1.1 and return wantarray ? () : [];
695 : golsen 1.3
696 : overbeek 1.1 my @queries = ();
697 :     my @qdata;
698 :     my @subjects = ();
699 :     my @sdata;
700 :     my @hsps = ();
701 :    
702 :     local $_;
703 :     while ( defined( $_ = <$fh> ) )
704 :     {
705 :     chomp;
706 :     my ( $type, @datum ) = split /\t/;
707 :     if ( $type eq 'Query' )
708 :     {
709 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
710 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
711 :     @qdata = @datum;
712 :     @sdata = ();
713 :     @subjects = ();
714 :     @hsps = ();
715 :     }
716 :     elsif ( $type eq '>' )
717 :     {
718 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
719 :     @sdata = @datum;
720 :     @hsps = ();
721 :     }
722 :     elsif ( $type eq 'HSP' )
723 :     {
724 :     push @hsps, \@datum;
725 :     }
726 :     }
727 :    
728 :     close $fh if $close;
729 :    
730 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
731 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
732 :    
733 :     wantarray ? @queries : \@queries;
734 :     }
735 :    
736 :    
737 : golsen 1.3 sub verify_db
738 :     {
739 :     my ( $db, $type ) = @_;
740 :    
741 :     my @args;
742 :     if ($type =~ /^p/i && ((! -s "$db.psq") || (-M "$db.psq" > -M $db)))
743 :     {
744 :     @args = ( "-p", "T", "-i", $db );
745 :     }
746 :     elsif ((! -s "$db.nsq") || (-M "$db.nsq" > -M $db))
747 :     {
748 :     @args = ( "-p", "F", "-i", $db );
749 :     }
750 :     else
751 :     {
752 :     return 1;
753 :     }
754 :    
755 :     #
756 :     # Find formatdb; if we're operating in a SEED environment
757 :     # use it from there.
758 :     #
759 :    
760 :     my $prog = SeedAware::executable_for( 'formatdb' );
761 :     if ( ! $prog )
762 :     {
763 :     warn "gjoalignandtree::verify_db: formatdb program not found.\n";
764 :     return 0;
765 :     }
766 :    
767 :     my $rc = system( $prog, @args );
768 :    
769 :     if ( $rc != 0 )
770 :     {
771 :     my $cmd = join( ' ', $prog, @args );
772 :     warn "gjoalignandtree::verify_db: formatdb failed with rc = $rc: $cmd\n";
773 :     return 0;
774 :     }
775 :    
776 :     return 1;
777 :     }
778 :    
779 :    
780 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3