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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.15 - (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 : golsen 1.8 # Profile blast against protein database:
27 :     #
28 : golsen 1.3 # $structured_blast = blastpgp( $dbfile, $profilefile, \%options )
29 :     # $structured_blast = blastpgp( \@dbseq, $profilefile, \%options )
30 :     # $structured_blast = blastpgp( $dbfile, \@profileseq, \%options )
31 :     # $structured_blast = blastpgp( \@dbseq, \@profileseq, \%options )
32 :     #
33 : golsen 1.8 # Profile blast against DNA database (PSI-tblastn):
34 :     #
35 :     # $structured_blast = blastpgpn( $dbfile, $profilefile, \%options )
36 :     # $structured_blast = blastpgpn( \@dbseq, $profilefile, \%options )
37 :     # $structured_blast = blastpgpn( $dbfile, \@profileseq, \%options )
38 :     # $structured_blast = blastpgpn( \@dbseq, \@profileseq, \%options )
39 :     #
40 : golsen 1.3 # print_blast_as_records( $file, \@queries )
41 :     # print_blast_as_records( \*FH, \@queries )
42 :     # print_blast_as_records( \@queries ) # STDOUT
43 :     #
44 :     # \@queries = read_blast_from_records( $file )
45 :     # \@queries = read_blast_from_records( \*FH )
46 :     # \@queries = read_blast_from_records( ) # STDIN
47 :     #
48 :     # $exists = verify_db( $db, $type );
49 :     #
50 :     #-------------------------------------------------------------------------------
51 : overbeek 1.1
52 :     use strict;
53 : golsen 1.3 use SeedAware;
54 : overbeek 1.1 use gjoseqlib;
55 : golsen 1.3 use gjoparseblast;
56 : overbeek 1.1
57 :     #-------------------------------------------------------------------------------
58 :     # Order a set of sequences:
59 :     #
60 :     # @seqs = order_sequences_by_length( \@seqs, \%opts )
61 :     # \@seqs = order_sequences_by_length( \@seqs, \%opts )
62 :     #
63 :     # Options:
64 :     #
65 : golsen 1.3 # order => $key # increasing (D), decreasing, median_outward,
66 :     # # closest_to_median, closest_to_length
67 :     # length => $prefered_length
68 : overbeek 1.1 #
69 :     #-------------------------------------------------------------------------------
70 :    
71 :     sub order_sequences_by_length
72 :     {
73 :     my ( $seqs, $opts ) = @_;
74 :     $seqs && ref $seqs eq 'ARRAY' && @$seqs
75 :     or print STDERR "order_sequences_by_length called with invalid sequence list.\n"
76 :     and return wantarray ? () : [];
77 :    
78 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
79 :     my $order = $opts->{ order } || 'increasing';
80 :    
81 :     my @seqs = ();
82 :     if ( $order =~ /^dec/i )
83 :     {
84 :     $opts->{ order } = 'decreaseing';
85 :     @seqs = sort { length( $b->[2] ) <=> length( $a->[2] ) } @$seqs;
86 :     }
87 :     elsif ( $order =~ /^med/i )
88 :     {
89 :     $opts->{ order } = 'median_outward';
90 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
91 :     while ( $_ = shift @seqs0 )
92 :     {
93 :     push @seqs, $_;
94 :     push @seqs, pop @seqs0 if @seqs0;
95 :     }
96 :     @seqs = reverse @seqs;
97 :     }
98 :     elsif ( $order =~ /^close/i )
99 :     {
100 :     my $pref_len;
101 :     if ( defined $opts->{ length } ) # caller supplies preferred length?
102 :     {
103 :     $opts->{ order } = 'closest_to_length';
104 :     $pref_len = $opts->{ length } + 0;
105 :     }
106 :     else # preferred length is median
107 :     {
108 :     $opts->{ order } = 'closest_to_median';
109 :     my @seqs0 = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
110 :     my $i = int( @seqs0 / 2 );
111 :     $pref_len = ( @seqs0 % 2 ) ? length( $seqs0[$i]->[2] )
112 :     : ( length( $seqs0[$i-1]->[2] ) + length( $seqs0[$i]->[2] ) ) / 2;
113 :     }
114 :    
115 :     @seqs = map { $_->[0] }
116 :     sort { $a->[1] <=> $b->[1] } # closest to preferred first
117 :     map { [ $_, abs( length( $_->[2] ) - $pref_len ) ] } # dist from pref?
118 :     @$seqs;
119 :     }
120 :     else
121 :     {
122 :     $opts->{ order } = 'increasing';
123 :     @seqs = sort { length( $a->[2] ) <=> length( $b->[2] ) } @$seqs;
124 :     }
125 :    
126 :     wantarray ? @seqs : \@seqs;
127 :     }
128 :    
129 :    
130 :     #-------------------------------------------------------------------------------
131 :     # Write representative seqeunce groups to a file
132 :     #
133 :     # write_representative_groups( \%rep, $file )
134 :     #
135 :     #-------------------------------------------------------------------------------
136 :    
137 :     sub write_representative_groups
138 :     {
139 :     my ( $repH, $file ) = @_;
140 :     $repH && (ref $repH eq 'HASH') && keys %$repH
141 :     or print STDERR "write_representative_groups called with invalid rep hash.\n"
142 :     and return;
143 :    
144 :     my ( $fh, $close ) = &output_file_handle( $file );
145 :    
146 :     # Order keys from largest to smallest set, then alphabetical
147 :    
148 :     my @keys = sort { @{ $repH->{$b} } <=> @{ $repH->{$a} } || lc $a cmp lc $b } keys %$repH;
149 :    
150 :     foreach ( @keys ) { print $fh join( "\t", ( $_, @{ $repH->{ $_ } } ) ), "\n" }
151 :    
152 :     close $fh if $close;
153 :     }
154 :    
155 :    
156 :     #-------------------------------------------------------------------------------
157 :     #
158 :     # @align = trim_align_to_median_ends( \@align, \%opts )
159 :     # \@align = trim_align_to_median_ends( \@align, \%opts )
160 :     #
161 :     # Options:
162 :     #
163 : fangfang 1.5 # begin => bool # Trim start (specifically)
164 :     # end => bool # Trim end (specifically)
165 :     # fract_cov => fract # Fraction of sequences to be covered (D: 0.75)
166 : overbeek 1.1 #
167 :     #-------------------------------------------------------------------------------
168 :     sub trim_align_to_median_ends
169 :     {
170 :     my ( $align, $opts ) = @_;
171 : golsen 1.8
172 : overbeek 1.1 $align && ref $align eq 'ARRAY' && @$align
173 :     or print STDERR "trim_align_to_median_ends called with invalid sequence list.\n"
174 :     and return wantarray ? () : [];
175 :    
176 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
177 :     my $tr_beg = $opts->{ begin } || $opts->{ beg } || $opts->{ start } ? 1 : 0;
178 :     my $tr_end = $opts->{ end } ? 1 : 0;
179 :     $tr_beg = $tr_end = 1 if ! ( $tr_beg || $tr_end );
180 : fangfang 1.5 my $frac = $opts->{ fract_cov } || 0.75;
181 : overbeek 1.1
182 : fangfang 1.5 my( @ngap1, @ngap2);
183 : overbeek 1.1 foreach my $seq ( @$align )
184 :     {
185 :     my( $b, $e ) = $seq->[2] =~ /^(-*).*[^-](-*)$/;
186 : fangfang 1.5 push @ngap1, length( $b || '' );
187 :     push @ngap2, length( $e || '' );
188 : overbeek 1.1 }
189 :    
190 : fangfang 1.5 @ngap1 = sort { $a <=> $b } @ngap1;
191 :     @ngap2 = sort { $a <=> $b } @ngap2;
192 : overbeek 1.1
193 : fangfang 1.5 my $ngap1 = $tr_beg ? $ngap1[ int( @ngap1 * $frac ) ] : 0;
194 :     my $ngap2 = $tr_end ? $ngap2[ int( @ngap2 * $frac ) ] : 0;
195 : overbeek 1.1
196 :     my $ori_len = length( $align->[0]->[2] );
197 : fangfang 1.5 my $new_len = $ori_len - ( $ngap1 + $ngap2 );
198 :     my @align2 = map { [ @$_[0,1], substr( $_->[2], $ngap1, $new_len ) ] }
199 : overbeek 1.1 @$align;
200 :    
201 :     wantarray ? @align2 : \@align2;
202 :     }
203 :    
204 :    
205 :     #-------------------------------------------------------------------------------
206 :     #
207 : golsen 1.3 # print_alignment_as_pseudoclustal( $align, \%opts )
208 :     #
209 : overbeek 1.1 # Options:
210 :     #
211 :     # file => $filename # supply a file name to open and write
212 :     # file => \*FH # supply an open file handle (D = STDOUT)
213 :     # line => $linelen # residues per line (D = 60)
214 :     # lower => $bool # all lower case sequence
215 :     # upper => $bool # all upper case sequence
216 :     #
217 :     #-------------------------------------------------------------------------------
218 :    
219 :     sub print_alignment_as_pseudoclustal
220 :     {
221 :     my ( $align, $opts ) = @_;
222 :     $align && ref $align eq 'ARRAY' && @$align
223 :     or print STDERR "print_alignment_as_pseudoclustal called with invalid sequence list.\n"
224 :     and return wantarray ? () : [];
225 :    
226 :     $opts = {} if ! ( $opts && ref $opts eq 'HASH' );
227 :     my $line_len = $opts->{ line } || 60;
228 :     my $case = $opts->{ upper } ? 1 : $opts->{ lower } ? -1 : 0;
229 :    
230 :     my ( $fh, $close ) = &output_file_handle( $opts->{ file } );
231 :    
232 :     my $namelen = 0;
233 :     foreach ( @$align ) { $namelen = length $_->[0] if $namelen < length $_->[0] }
234 :     my $fmt = "%-${namelen}s %s\n";
235 :    
236 :     my $id;
237 :     my @lines = map { $id = $_->[0]; [ map { sprintf $fmt, $id, $_ }
238 : golsen 1.3 map { $case < 0 ? lc $_ : $case > 0 ? uc $_ : $_ } # map sequence only
239 : overbeek 1.1 $_->[2] =~ m/.{1,$line_len}/g
240 :     ] }
241 :    
242 :     @$align;
243 :    
244 :     my $ngroup = @{ $lines[0] };
245 :     for ( my $i = 0; $i < $ngroup; $i++ )
246 :     {
247 :     foreach ( @lines ) { print $fh $_->[$i] if $_->[$i] }
248 :     print $fh "\n";
249 :     }
250 :    
251 :     close $fh if $close;
252 :     }
253 :    
254 :    
255 :     #-------------------------------------------------------------------------------
256 : golsen 1.3 #
257 :     # @seqs = read_pseudoclustal( ) # D = STDIN
258 :     # \@seqs = read_pseudoclustal( ) # D = STDIN
259 :     # @seqs = read_pseudoclustal( $file_name )
260 :     # \@seqs = read_pseudoclustal( $file_name )
261 :     # @seqs = read_pseudoclustal( \*FH )
262 :     # \@seqs = read_pseudoclustal( \*FH )
263 :     #
264 :     #-------------------------------------------------------------------------------
265 :    
266 :     sub read_pseudoclustal
267 :     {
268 :     my ( $file ) = @_;
269 :     my ( $fh, $close ) = input_file_handle( $file );
270 :     my %seq;
271 :     my @ids;
272 :     while ( <$fh> )
273 :     {
274 :     chomp;
275 :     my ( $id, $data ) = /^(\S+)\s+(\S.*)$/;
276 :     if ( defined $id && defined $data )
277 :     {
278 :     push @ids, $id if ! $seq{ $id };
279 :     $data =~ s/\s+//g;
280 :     push @{ $seq{ $id } }, $data;
281 :     }
282 :     }
283 :     close $fh if $close;
284 :    
285 :     my @seq = map { [ $_, '', join( '', @{ $seq{ $_ } } ) ] } @ids;
286 :     wantarray ? @seq : \@seq;
287 :     }
288 :    
289 :    
290 :     #-------------------------------------------------------------------------------
291 : overbeek 1.1 # The profile 'query' sequence:
292 :     # 1. No terminal gaps, if possible
293 :     # 1.1 Otherwise, no gaps at start
294 :     # 2. Longest sequence passing above
295 :     #
296 :     # $prof_rep = representative_for_profile( $align )
297 :     #-------------------------------------------------------------------------------
298 :     sub representative_for_profile
299 :     {
300 :     my ( $align ) = @_;
301 :     $align && ref $align eq 'ARRAY' && @$align
302 :     or die "representative_for_profile called with invalid sequence list.\n";
303 :    
304 :     my @cand;
305 :     @cand = grep { $_->[2] =~ /^[^-].*[^-]$/ } @$align; # No end gaps
306 :     @cand = grep { $_->[2] =~ /^[^-]/ } @$align if ! @cand; # No start gaps
307 : fangfang 1.4 @cand = @$align if ! @cand;
308 : overbeek 1.1
309 : golsen 1.3 my ( $rep ) = map { $_->[0] } # sequence entry
310 :     sort { $b->[1] <=> $a->[1] } # max nongaps
311 :     map { [ $_, $_->[2] =~ tr/-//c ] } # count nongaps
312 : overbeek 1.1 @cand;
313 :    
314 :     $rep = [ @$rep ]; # Make a copy
315 :     $rep->[2] =~ s/-+//g; # Compress sequence
316 :    
317 :     return $rep;
318 :     }
319 :    
320 :    
321 :     #-------------------------------------------------------------------------------
322 :     #
323 : golsen 1.3 # \@seq = extract_with_psiblast( $db, $profile, \%opts )
324 :     # ( \@seq, \%report ) = extract_with_psiblast( $db, $profile, \%opts )
325 :     #
326 :     # $db can be $db_file_name or \@db_seqs
327 :     # $profile can be $profile_file_name or \@profile_seqs
328 :     #
329 :     # If supplied as file, profile is pseudoclustal
330 :     #
331 :     # Report records:
332 :     #
333 : fangfang 1.13 # [ $sid, $score, $e_value, $slen, $status,
334 :     # $frac_id, $frac_pos,
335 : golsen 1.3 # $q_uncov_n_term, $q_uncov_c_term,
336 :     # $s_uncov_n_term, $s_uncov_c_term ]
337 :     #
338 :     # Options:
339 :     #
340 :     # e_value => $max_e_value # maximum blastpgp E-value (D = 0.01)
341 : golsen 1.8 # max_e_value => $max_e_value # maximum blastpgp E-value (D = 0.01)
342 : golsen 1.3 # max_q_uncov => $aa_per_end # maximum unmatched query (D = 20)
343 :     # max_q_uncov_c => $aa_per_end # maximum unmatched query, c-term (D = 20)
344 :     # max_q_uncov_n => $aa_per_end # maximum unmatched query, n-term (D = 20)
345 :     # min_ident => $frac_ident # minimum fraction identity (D = 0.15)
346 :     # min_positive => $frac_positive # minimum fraction positive scoring (D = 0.20)
347 : fangfang 1.11 # min_frac_cov => $frac_cov # minimum fraction coverage of query and subject sequence (D = 0.20)
348 : fangfang 1.10 # min_q_cov => $frac_cov # minimum fraction coverage of query sequence (D = 0.50)
349 : fangfang 1.11 # min_s_cov => $frac_cov # minimum fraction coverage of subject sequence (D = 0.01)
350 : golsen 1.8 # n_result => $max_results # maxiumn restuls returned (D = 1000)
351 :     # n_thread => $n_thread # number of blastpgp threads (D = 2)
352 : golsen 1.3 # query => $q_file_name # query sequence file (D = most complete)
353 :     # query => \@q_seq_entry # query sequence (D = most complete)
354 : fangfang 1.15 # query_used => \@q_seq_entry # output
355 : golsen 1.3 # stderr => $file # blastpgp stderr (D = /dev/stderr)
356 :     #
357 :     # If supplied, query must be identical to a profile sequence, but no gaps.
358 :     #
359 :     #-------------------------------------------------------------------------------
360 :    
361 :     sub extract_with_psiblast
362 :     {
363 : fangfang 1.6 my ( $seqs, $profile, $opts ) = @_;
364 : golsen 1.3
365 :     $opts && ref $opts eq 'HASH' or $opts = {};
366 :    
367 : golsen 1.8 $opts->{ e_value } ||= $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
368 : fangfang 1.12 $opts->{ n_result } ||= $opts->{ nresult } || 1000;
369 : fangfang 1.15 $opts->{ n_thread } ||= $opts->{ nthread } || 2;
370 : fangfang 1.11 my $max_q_uncov_c = $opts->{ max_q_uncov_c } || $opts->{ max_q_uncov } || 20;
371 :     my $max_q_uncov_n = $opts->{ max_q_uncov_n } || $opts->{ max_q_uncov } || 20;
372 :     my $min_q_cov = $opts->{ min_q_cov } || $opts->{ min_frac_cov } || 0.50;
373 :     my $min_s_cov = $opts->{ min_s_cov } || $opts->{ min_frac_cov } || 0.01;
374 :     my $min_ident = $opts->{ min_ident } || 0.15;
375 :     my $min_pos = $opts->{ min_positive } || 0.20;
376 : golsen 1.3
377 : fangfang 1.6 my $blast = blastpgp( $seqs, $profile, $opts );
378 : golsen 1.3 $blast && @$blast or return ();
379 :    
380 :     my ( $qid, $qdef, $qlen, $qhits ) = @{ $blast->[0] };
381 :    
382 :     my @trimmed;
383 :     my %report;
384 :    
385 :     foreach my $sdata ( @$qhits )
386 :     {
387 :     my ( $sid, $sdef, $slen, $hsps ) = @$sdata;
388 :    
389 :     if ( one_real_hsp( $hsps ) )
390 :     {
391 :     # [ scr, exp, p_n, pval, nmat, nid, nsim, ngap, dir, q1, q2, qseq, s1, s2, sseq ]
392 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
393 :    
394 :     my $hsp0 = $hsps->[0];
395 : fangfang 1.13 my ( $scr, $exp, $nmat, $nid, $npos, $q1, $q2, $qseq, $s1, $s2, $sseq ) = ( @$hsp0 )[ 0, 1, 4, 5, 6, 9, 10, 11, 12, 13, 14 ];
396 : golsen 1.3
397 :     my $status;
398 : fangfang 1.10 if ( $q1-1 > $max_q_uncov_n ) { $status = 'missing start' }
399 :     elsif ( $qlen-$q2 > $max_q_uncov_c ) { $status = 'missing end' }
400 :     elsif ( $nid / $nmat < $min_ident ) { $status = 'low identity' }
401 :     elsif ( $npos / $nmat < $min_pos ) { $status = 'low positives' }
402 :     elsif ( ($q2-$q1+1)/$qlen < $min_q_cov ) { $status = 'low coverage' }
403 :     elsif ( ($s2-$s1+1)/$slen < $min_s_cov ) { $status = 'long subject' }
404 : golsen 1.3 else
405 :     {
406 :     $sseq =~ s/-+//g;
407 :     $sdef .= " ($s1-$s2/$slen)" if ( ( $s1 > 1 ) || ( $s2 < $slen ) );
408 :     push @trimmed, [ $sid, $sdef, $sseq ];
409 :     $status = 'included';
410 :     }
411 :    
412 : fangfang 1.13 my $frac_id = sprintf("%.3f", $nid/$nmat);
413 :     my $frac_pos = sprintf("%.3f", $npos/$nmat);
414 :    
415 :     $report{ $sid } = [ $sid, $scr, $exp, $slen, $status, $frac_id, $frac_pos, $q1-1, $qlen-$q2, $s1-1, $slen-$s2 ];
416 : golsen 1.3 }
417 :     }
418 :    
419 :     wantarray ? ( \@trimmed, \%report ) : \@trimmed;
420 :     }
421 :    
422 :    
423 :     #-------------------------------------------------------------------------------
424 :     #
425 :     # Allow fragmentary matches inside of the highest-scoring hsp:
426 :     #
427 :     #-------------------------------------------------------------------------------
428 :    
429 :     sub one_real_hsp
430 :     {
431 :     my ( $hsps ) = @_;
432 :     return 0 if ! ( $hsps && ( ref( $hsps ) eq 'ARRAY' ) && @$hsps );
433 :     return 1 if @$hsps == 1;
434 :    
435 :     my ( $q1_0, $q2_0 ) = ( @{ $hsps->[0] } )[9, 10];
436 :     for ( my $i = 1; $i < @$hsps; $i++ )
437 :     {
438 :     my ($q1, $q2) = ( @{ $hsps->[$i] } )[9, 10];
439 :     return 0 if $q1 < $q1_0 || $q2 > $q2_0;
440 :     }
441 :    
442 :     return 1;
443 :     }
444 :    
445 :    
446 :     #-------------------------------------------------------------------------------
447 :     #
448 : golsen 1.8 # $structured_blast = blastpgp( $dbfile, $profilefile, \%options )
449 :     # $structured_blast = blastpgp( \@dbseq, $profilefile, \%options )
450 :     # $structured_blast = blastpgp( $dbfile, \@profileseq, \%options )
451 :     # $structured_blast = blastpgp( \@dbseq, \@profileseq, \%options )
452 : overbeek 1.1 #
453 :     # Required:
454 :     #
455 : golsen 1.3 # $db_file or \@db_seq
456 :     # $profile_file or \@profile_seq
457 :     #
458 :     # Options:
459 :     #
460 : golsen 1.8 # e_value => $max_e_value # maximum E-value of matches (D = 0.01)
461 :     # max_e_value => $max_e_value # maximum E-value of matches (D = 0.01)
462 :     # n_result => $max_results # maximum matches returned (D = 1000)
463 :     # n_thread => $n_thread # number of blastpgp threads (D = 2)
464 :     # stderr => $file # place to send blast stderr (D = /dev/null)
465 :     # query => $q_file_name # most complete
466 :     # query => \@q_seq_entry # most complete
467 : fangfang 1.14 # query_used => \@q_seq_entry # output
468 : overbeek 1.1 #
469 :     #-------------------------------------------------------------------------------
470 :    
471 :     sub blastpgp
472 :     {
473 :     my ( $db, $profile, $opts ) = @_;
474 :     $opts ||= {};
475 : golsen 1.8 my $tmp = SeedAware::location_of_tmp( $opts );
476 : overbeek 1.1
477 :     my ( $dbfile, $rm_db );
478 :     if ( defined $db && ref $db )
479 :     {
480 :     ref $db eq 'ARRAY'
481 :     && @$db
482 :     || print STDERR "blastpgp requires one or more database sequences.\n"
483 :     && return undef;
484 : golsen 1.8 $dbfile = SeedAware::new_file_name( "$tmp/blastpgp_db" );
485 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $dbfile, $db );
486 :     $rm_db = 1;
487 :     }
488 :     elsif ( defined $db && -f $db )
489 :     {
490 :     $dbfile = $db;
491 :     }
492 :     else
493 :     {
494 :     die "blastpgp requires database.";
495 :     }
496 : golsen 1.3 verify_db( $dbfile, 'P' ); # protein
497 : overbeek 1.1
498 :     my ( $proffile, $rm_profile );
499 :     if ( defined $profile && ref $profile )
500 :     {
501 :     ref $profile eq 'ARRAY'
502 :     && @$profile
503 :     || print STDERR "blastpgp requires one or more profile sequences.\n"
504 :     && return undef;
505 : golsen 1.8 $proffile = SeedAware::new_file_name( "$tmp/blastpgp_profile" );
506 : overbeek 1.1 &print_alignment_as_pseudoclustal( $profile, { file => $proffile, upper => 1 } );
507 :     $rm_profile = 1;
508 :     }
509 :     elsif ( defined $profile && -f $profile )
510 :     {
511 :     $proffile = $profile;
512 :     }
513 :     else
514 :     {
515 :     die "blastpgp requires profile.";
516 :     }
517 :    
518 :     my ( $qfile, $rm_query );
519 :     my $query = $opts->{ query };
520 :     if ( defined $query && ref $query )
521 :     {
522 : golsen 1.3 ref $query eq 'ARRAY' && @$query == 3
523 :     or print STDERR "blastpgp invalid query sequence.\n"
524 :     and return undef;
525 :    
526 : golsen 1.8 $qfile = SeedAware::new_file_name( "$tmp/blastpgp_query" );
527 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
528 :     $rm_query = 1;
529 :     }
530 :     elsif ( defined $query && -f $query )
531 :     {
532 :     $qfile = $query;
533 :     }
534 :     elsif ( $profile ) # Build it from profile
535 :     {
536 :     $query = &representative_for_profile( $profile );
537 : golsen 1.8 $qfile = SeedAware::new_file_name( "$tmp/blastpgp_query" );
538 : overbeek 1.1 gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
539 :     $rm_query = 1;
540 :     }
541 :     else
542 : golsen 1.3 {
543 : overbeek 1.1 die "blastpgp requires database.";
544 :     }
545 :    
546 : fangfang 1.15 ($query) = gjoseqlib::read_fasta($qfile);
547 :     $opts->{ query_used } = $query;
548 : fangfang 1.13
549 : golsen 1.8 my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
550 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 2;
551 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
552 : fangfang 1.13
553 : golsen 1.9 my $blastpgp = SeedAware::executable_for( 'blastpgp' )
554 :     or print STDERR "Could not find executable for program 'blastpgp'.\n"
555 :     and return undef;
556 :    
557 :     my @cmd = ( $blastpgp,
558 : golsen 1.8 '-j', 1, # function is profile alignment
559 :     '-B', $proffile, # location of profile
560 : golsen 1.3 '-d', $dbfile,
561 : golsen 1.8 '-i', $qfile,
562 :     '-e', $e_val,
563 : golsen 1.3 '-b', $nkeep,
564 :     '-v', $nkeep,
565 : golsen 1.8 '-t', 1 # issues warning if not set
566 : golsen 1.3 );
567 : golsen 1.8 push @cmd, ( '-a', $n_cpu ) if $n_cpu;
568 : golsen 1.3
569 : golsen 1.8 my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
570 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
571 : golsen 1.3 or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
572 :     and return undef;
573 :     my $out = &gjoparseblast::structured_blast_output( $blastfh, 1 ); # selfmatches okay
574 :     close $blastfh;
575 : overbeek 1.1
576 :     if ( $rm_db )
577 :     {
578 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
579 :     unlink @files if @files;
580 :     }
581 :     unlink $proffile if $rm_profile;
582 :     unlink $qfile if $rm_query;
583 :    
584 :     return $out;
585 :     }
586 :    
587 :    
588 :     #-------------------------------------------------------------------------------
589 : golsen 1.8 # Do psiblast against tranlated genomic DNA
590 :     #
591 :     # $structured_blast = blastpgpn( $dbfile, $profilefile, \%options )
592 :     # $structured_blast = blastpgpn( \@dbseq, $profilefile, \%options )
593 :     # $structured_blast = blastpgpn( $dbfile, \@profileseq, \%options )
594 :     # $structured_blast = blastpgpn( \@dbseq, \@profileseq, \%options )
595 :     #
596 :     # Required:
597 :     #
598 :     # $db_file or \@db_seq
599 :     # $profile_file or \@profile_seq
600 :     #
601 :     # Options:
602 :     #
603 :     # aa_db_file => $trans_file # put translation db here
604 :     # e_value => $max_e_value # D = 0.01
605 :     # max_e_val => $max_e_value # D = 0.01
606 :     # n_cpu => $n_cpu # synonym of n_thread
607 :     # n_result => $max_seq # D = 1000
608 :     # n_thread => $n_cpu # D = 1
609 :     # query => $q_file_name # most complete
610 :     # query => \@q_seq_entry # most complete
611 : fangfang 1.14 # query_used => \@q_seq_entry # output
612 : golsen 1.8 # stderr => $file # place to send blast stderr (D = /dev/null)
613 :     # tmp => $temp_dir # location for temp files
614 :     #
615 :     # depricated alternatives:
616 :     #
617 :     # prot_file => $trans_file # put translation db here
618 :     #-------------------------------------------------------------------------------
619 :    
620 :     sub blastpgpn
621 :     {
622 :     my ( $ndb, $profile, $opts ) = @_;
623 :     $opts ||= {};
624 :     my $tmp = SeedAware::location_of_tmp( $opts );
625 :    
626 :     $opts->{ aa_db_file } ||= $opts->{ prot_file } if defined $opts->{ prot_file };
627 :     my $dbfile = $opts->{ aa_db_file } || SeedAware::new_file_name( "$tmp/blastpgpn_db" );
628 :     my $rm_db = ! $opts->{ aa_db_file };
629 :     if ( defined $dbfile && -f $dbfile && -s $dbfile )
630 :     {
631 :     # The tranaslated sequence database exists
632 :     }
633 :     elsif ( defined $ndb )
634 :     {
635 :     if ( ref $ndb eq 'ARRAY' && @$ndb )
636 :     {
637 :     ref $ndb eq 'ARRAY'
638 :     && @$ndb
639 :     || print STDERR "Bad sequence reference passed to blastpgpn.\n"
640 :     && return undef;
641 :     my @pdb = map { six_translations( $_ ) } @$ndb;
642 :     gjoseqlib::print_alignment_as_fasta( $dbfile, \@pdb );
643 :     }
644 :     elsif ( -f $ndb && -s $ndb )
645 :     {
646 :     open( PDB, ">$dbfile" )
647 :     or print STDERR "Could not open protein database file '$dbfile'.\n"
648 :     and return undef;
649 :     my $entry;
650 :     while ( $entry = gjoseqlib::next_fasta_entry( $ndb ) )
651 :     {
652 :     gjoseqlib::write_alignment_as_fasta( \*PDB, six_translations( $entry ) );
653 :     }
654 :     close PDB;
655 :     }
656 :     else
657 :     {
658 :     die "blastpgpn requires a sequence database.";
659 :     }
660 :     }
661 :     verify_db( $dbfile, 'P' ); # protein
662 :    
663 :     my ( $proffile, $rm_profile );
664 :     if ( defined $profile && ref $profile )
665 :     {
666 :     ref $profile eq 'ARRAY'
667 :     && @$profile
668 :     || print STDERR "blastpgpn requires one or more profile sequences.\n"
669 :     && return undef;
670 :     $proffile = SeedAware::new_file_name( "$tmp/blastpgpn_profile" );
671 :     &print_alignment_as_pseudoclustal( $profile, { file => $proffile, upper => 1 } );
672 :     $rm_profile = 1;
673 :     }
674 :     elsif ( defined $profile && -f $profile )
675 :     {
676 :     $proffile = $profile;
677 :     }
678 :     else
679 :     {
680 :     die "blastpgpn requires profile.";
681 :     }
682 :    
683 :     my ( $qfile, $rm_query );
684 :     my $query = $opts->{ query };
685 :     if ( defined $query && ref $query )
686 :     {
687 :     ref $query eq 'ARRAY' && @$query == 3
688 :     or print STDERR "blastpgpn invalid query sequence.\n"
689 :     and return undef;
690 :    
691 :     $qfile = SeedAware::new_file_name( "$tmp/blastpgpn_query" );
692 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
693 :     $rm_query = 1;
694 :     }
695 :     elsif ( defined $query && -f $query )
696 :     {
697 :     $qfile = $query;
698 :     }
699 :     elsif ( $profile ) # Build it from profile
700 :     {
701 :     $query = &representative_for_profile( $profile );
702 :     $qfile = SeedAware::new_file_name( "$tmp/blastpgpn_query" );
703 :     gjoseqlib::print_alignment_as_fasta( $qfile, [$query] );
704 :     $rm_query = 1;
705 :     }
706 :     else
707 :     {
708 :     die "blastpgpn requires database.";
709 :     }
710 :    
711 : fangfang 1.15 ($query) = gjoseqlib::read_fasta($qfile);
712 :     $opts->{ query_used } = $query;
713 : fangfang 1.13
714 : golsen 1.8 my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
715 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 2;
716 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
717 :    
718 : golsen 1.9 my $blastpgp = SeedAware::executable_for( 'blastpgp' )
719 :     or print STDERR "Could not find executable for program 'blastpgp'.\n"
720 :     and return undef;
721 :    
722 :     my @cmd = ( $blastpgp,
723 : golsen 1.8 '-j', 1, # function is profile alignment
724 :     '-B', $proffile, # location of protein profile
725 :     '-d', $dbfile, # protein database
726 :     '-i', $qfile, # one protein from the profile
727 :     '-e', $e_val,
728 :     '-b', $nkeep,
729 :     '-v', $nkeep,
730 :     '-t', 1, # issues warning if not set
731 :     );
732 :     push @cmd, ( '-a', $n_cpu ) if $n_cpu;
733 :    
734 :     my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
735 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
736 :     or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
737 :     and return undef;
738 :     my $out = &gjoparseblast::structured_blast_output( $blastfh );
739 :     close $blastfh;
740 :    
741 :     if ( $rm_db )
742 :     {
743 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
744 :     unlink @files if @files;
745 :     }
746 :     unlink $proffile if $rm_profile;
747 :     unlink $qfile if $rm_query;
748 :    
749 :     # Fix the blastp output to look like tblastn output
750 :    
751 :     foreach my $qdata ( @$out )
752 :     {
753 :     my %sdata; # There are now multiple sequences per subject DNA
754 :     foreach my $sdata ( @{ $qdata->[3] } )
755 :     {
756 :     my ( $sid, $sdef, undef, $hsps ) = @$sdata;
757 :     my $fr;
758 :     ( $sid, $fr ) = $sid =~ m/^(.*)\.([-+]\d)$/;
759 :     my ( $b, $e, $slen ) = $sdef =~ m/(\d+)-(\d+)\/(\d+)$/;
760 :     $sdef =~ s/ ?\S+$//;
761 :     my $sd2 = $sdata{ $sid } ||= [ $sid, $sdef, $slen, [] ];
762 :     foreach ( @$hsps ) { adjust_hsp( $_, $fr, $b ); push @{$sd2->[3]}, $_ }
763 :     }
764 :    
765 :     # Order the hsps for each subject
766 :    
767 :     foreach my $sid ( keys %sdata )
768 :     {
769 :     my $hsps = $sdata{ $sid }->[3];
770 :     @$hsps = sort { $b->[0] <=> $a->[0] } @$hsps;
771 :     }
772 :    
773 :     # Order the subjects for the query
774 :    
775 :     @{$qdata->[3]} = map { $_->[0] } # remove score
776 :     sort { $b->[1] <=> $a->[1] } # sort by score
777 :     map { [ $_, $_->[3]->[0]->[0] ] } # tag with top score
778 :     map { $sdata{ $_ } } # subject data
779 :     keys %sdata; # subject ids
780 :     }
781 :    
782 :     return $out;
783 :     }
784 :    
785 :    
786 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
787 :     # When search is versus six frame translation, there is a need to adjust
788 :     # the frame and location information in the hsp back to the DNA coordinates.
789 :     #
790 :     # adjust_hsp( $hsp, $frame, $begin )
791 :     #
792 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
793 :     # scr Eval nseg Eval naln nid npos ngap fr q1 q2 qseq s1 s2 sseq
794 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
795 :     sub adjust_hsp
796 :     {
797 :     my ( $hsp, $fr, $b ) = @_;
798 :     $hsp->[8] = $fr;
799 :     if ( $fr > 0 )
800 :     {
801 :     $hsp->[12] = $b + 3 * ( $hsp->[12] - 1 );
802 :     $hsp->[13] = $b + 3 * ( $hsp->[13] - 1 ) + 2;
803 :     }
804 :     else
805 :     {
806 :     $hsp->[12] = $b - 3 * ( $hsp->[12] - 1 );
807 :     $hsp->[13] = $b - 3 * ( $hsp->[13] - 1 ) - 2;
808 :     }
809 :     }
810 :    
811 :    
812 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
813 :     # Do a six frame translation for use by blastpgpn. This modifications of
814 :     # the identifiers and definitions are essential to the interpretation of
815 :     # the blast results. The program 'translate_fasta_6' produces the same
816 :     # output format, and is much faster.
817 :     #
818 :     # @translations = six_translations( $nucleotide_entry )
819 :     #
820 :     # The ids are modified by adding ".frame" (+1, +2, +3, -1, -2, -3).
821 :     # The definition is mofidified by adding " begin-end/of_length".
822 :     # NCBI frame numbers reverse strand translation frames from the end of the
823 :     # sequence (i.e., the beginning of the complement of the strand).
824 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
825 :    
826 :     sub six_translations
827 :     {
828 :     my ( $id, $def, $seq ) = map { defined($_) ? $_ : '' } @{ $_[0] };
829 :     my $l = length( $seq );
830 :    
831 :     return () if $l < 15;
832 :    
833 :     # fr beg end
834 :     my @intervals = ( [ '+1', 1, $l - ( $l % 3 ) ],
835 :     [ '+2', 2, $l - ( ($l-1) % 3 ) ],
836 :     [ '+3', 3, $l - ( ($l-2) % 3 ) ],
837 :     [ '-1', $l, 1 + ( $l % 3 ) ],
838 :     [ '-2', $l-1, 1 + ( ($l-1) % 3 ) ],
839 :     [ '-3', $l-2, 1 + ( ($l-2) % 3 ) ]
840 :     );
841 :     my ( $fr, $b, $e );
842 :    
843 :     map { ( $fr, $b, $e ) = @$_;
844 :     [ "$id.$fr", "$def $b-$e/$l", gjoseqlib::translate_seq( gjoseqlib::DNA_subseq( \$seq, $b, $e ) ) ]
845 :     } @intervals;
846 :     }
847 :    
848 :    
849 :     #-------------------------------------------------------------------------------
850 : golsen 1.3 # Get an input file handle, and boolean on whether to close or not:
851 :     #
852 :     # ( \*FH, $close ) = input_file_handle( $filename );
853 :     # ( \*FH, $close ) = input_file_handle( \*FH );
854 :     # ( \*FH, $close ) = input_file_handle( ); # D = STDIN
855 :     #
856 :     #-------------------------------------------------------------------------------
857 :    
858 :     sub input_file_handle
859 :     {
860 :     my ( $file, $umask ) = @_;
861 :    
862 :     my ( $fh, $close );
863 :     if ( defined $file )
864 :     {
865 :     if ( ref $file eq 'GLOB' )
866 :     {
867 :     $fh = $file;
868 :     $close = 0;
869 :     }
870 :     elsif ( -f $file )
871 :     {
872 :     open( $fh, "<$file") || die "input_file_handle could not open '$file'.\n";
873 :     $close = 1;
874 :     }
875 :     else
876 :     {
877 :     die "input_file_handle could not find file '$file'.\n";
878 :     }
879 :     }
880 :     else
881 :     {
882 :     $fh = \*STDIN;
883 :     $close = 0;
884 :     }
885 :    
886 :     return ( $fh, $close );
887 :     }
888 :    
889 :    
890 :     #-------------------------------------------------------------------------------
891 : overbeek 1.1 # Get an output file handle, and boolean on whether to close or not:
892 :     #
893 :     # ( \*FH, $close ) = output_file_handle( $filename );
894 :     # ( \*FH, $close ) = output_file_handle( \*FH );
895 :     # ( \*FH, $close ) = output_file_handle( ); # D = STDOUT
896 :     #
897 :     #-------------------------------------------------------------------------------
898 :    
899 :     sub output_file_handle
900 :     {
901 :     my ( $file, $umask ) = @_;
902 :    
903 :     my ( $fh, $close );
904 :     if ( defined $file )
905 :     {
906 :     if ( ref $file eq 'GLOB' )
907 :     {
908 :     $fh = $file;
909 :     $close = 0;
910 :     }
911 :     else
912 :     {
913 :     open( $fh, ">$file") || die "output_file_handle could not open '$file'.\n";
914 :     chmod 0664, $file; # Seems to work on open file!
915 :     $close = 1;
916 :     }
917 :     }
918 :     else
919 :     {
920 :     $fh = \*STDOUT;
921 :     $close = 0;
922 :     }
923 :    
924 :     return ( $fh, $close );
925 :     }
926 :    
927 :    
928 :     #-------------------------------------------------------------------------------
929 :     #
930 :     # print_blast_as_records( $file, \@queries )
931 :     # print_blast_as_records( \*FH, \@queries )
932 :     # print_blast_as_records( \@queries ) # STDOUT
933 :     #
934 :     # \@queries = read_blast_from_records( $file )
935 :     # \@queries = read_blast_from_records( \*FH )
936 :     # \@queries = read_blast_from_records( ) # STDIN
937 :     #
938 :     # Three output record types:
939 :     #
940 :     # Query \t qid \t qdef \t dlen
941 :     # > \t sid \t sdef \t slen
942 :     # HSP \t ...
943 :     #
944 :     #-------------------------------------------------------------------------------
945 :     sub print_blast_as_records
946 :     {
947 :     my $file = ( $_[0] && ! ( ref $_[0] eq 'ARRAY' ) ? shift : undef );
948 :     my ( $fh, $close ) = output_file_handle( $file );
949 :    
950 :     my $queries = shift;
951 :     $queries && ref $queries eq 'ARRAY'
952 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
953 :     and return;
954 :    
955 :     foreach my $qdata ( @$queries )
956 :     {
957 :     $qdata && ref $qdata eq 'ARRAY' && @$qdata == 4
958 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
959 :     and return;
960 :    
961 :     my $subjcts = pop @$qdata;
962 :     $subjcts && ref $subjcts eq 'ARRAY'
963 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
964 :     and return;
965 :    
966 :     print $fh join( "\t", 'Query', @$qdata ), "\n";
967 :     foreach my $sdata ( @$subjcts )
968 :     {
969 :     $sdata && ref $sdata eq 'ARRAY' && @$sdata == 4
970 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
971 :     and return;
972 :    
973 :     my $hsps = pop @$sdata;
974 :     $hsps && ref $hsps eq 'ARRAY' && @$hsps
975 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
976 :     and return;
977 :    
978 :     print $fh join( "\t", '>', @$sdata ), "\n";
979 :    
980 :     foreach my $hsp ( @$hsps )
981 :     {
982 :     $hsp && ref $hsp eq 'ARRAY' && @$hsp > 4
983 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
984 :     and return;
985 :     print $fh join( "\t", 'HSP', @$hsp ), "\n";
986 :     }
987 :     }
988 :     }
989 :    
990 :     close $fh if $close;
991 :     }
992 :    
993 :    
994 :     sub read_blast_from_records
995 :     {
996 :     my ( $file ) = @_;
997 : golsen 1.3
998 :     my ( $fh, $close ) = input_file_handle( $file );
999 :     $fh or print STDERR "read_blast_from_records could not open input file.\n"
1000 : overbeek 1.1 and return wantarray ? () : [];
1001 : golsen 1.3
1002 : overbeek 1.1 my @queries = ();
1003 :     my @qdata;
1004 :     my @subjects = ();
1005 :     my @sdata;
1006 :     my @hsps = ();
1007 :    
1008 :     local $_;
1009 :     while ( defined( $_ = <$fh> ) )
1010 :     {
1011 :     chomp;
1012 :     my ( $type, @datum ) = split /\t/;
1013 :     if ( $type eq 'Query' )
1014 :     {
1015 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1016 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
1017 :     @qdata = @datum;
1018 :     @sdata = ();
1019 :     @subjects = ();
1020 :     @hsps = ();
1021 :     }
1022 :     elsif ( $type eq '>' )
1023 :     {
1024 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1025 :     @sdata = @datum;
1026 :     @hsps = ();
1027 :     }
1028 :     elsif ( $type eq 'HSP' )
1029 :     {
1030 :     push @hsps, \@datum;
1031 :     }
1032 :     }
1033 :    
1034 :     close $fh if $close;
1035 :    
1036 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1037 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
1038 :    
1039 :     wantarray ? @queries : \@queries;
1040 :     }
1041 :    
1042 :    
1043 : golsen 1.3 sub verify_db
1044 :     {
1045 :     my ( $db, $type ) = @_;
1046 :    
1047 :     my @args;
1048 : golsen 1.8 if ( $type =~ m/^p/i )
1049 : golsen 1.3 {
1050 : golsen 1.8 @args = ( "-p", "T", "-i", $db ) if ((! -s "$db.psq") || (-M "$db.psq" > -M $db));
1051 : golsen 1.3 }
1052 :     else
1053 :     {
1054 : golsen 1.8 @args = ( "-p", "F", "-i", $db ) if ((! -s "$db.nsq") || (-M "$db.nsq" > -M $db));
1055 : golsen 1.3 }
1056 : golsen 1.8 @args or return ( -s $db ? 1 : 0 );
1057 :    
1058 : golsen 1.3 #
1059 : golsen 1.8 # Find formatdb appropriate for the excecution environemnt.
1060 : golsen 1.3 #
1061 :    
1062 :     my $prog = SeedAware::executable_for( 'formatdb' );
1063 :     if ( ! $prog )
1064 :     {
1065 :     warn "gjoalignandtree::verify_db: formatdb program not found.\n";
1066 :     return 0;
1067 :     }
1068 :    
1069 :     my $rc = system( $prog, @args );
1070 :    
1071 :     if ( $rc != 0 )
1072 :     {
1073 :     my $cmd = join( ' ', $prog, @args );
1074 :     warn "gjoalignandtree::verify_db: formatdb failed with rc = $rc: $cmd\n";
1075 :     return 0;
1076 :     }
1077 :    
1078 :     return 1;
1079 :     }
1080 :    
1081 :    
1082 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3