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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3