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

Annotation of /FigKernelPackages/gjoalignandtree.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.13 - (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 :     # stderr => $file # blastpgp stderr (D = /dev/stderr)
355 :     #
356 :     # If supplied, query must be identical to a profile sequence, but no gaps.
357 : fangfang 1.13 # If not supplied, $opts{ query } will be set to the most compete sequence.
358 : golsen 1.3 #
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 : golsen 1.8 $opts->{ n_thread } ||= 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 : overbeek 1.1 #
468 : fangfang 1.13 # If not supplied, $opts{ query } will be set to the most compete sequence.
469 : overbeek 1.1 #-------------------------------------------------------------------------------
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.13 $opts->{ query } ||= $query;
547 :    
548 : golsen 1.8 my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
549 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 2;
550 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
551 : fangfang 1.13
552 : golsen 1.9 my $blastpgp = SeedAware::executable_for( 'blastpgp' )
553 :     or print STDERR "Could not find executable for program 'blastpgp'.\n"
554 :     and return undef;
555 :    
556 :     my @cmd = ( $blastpgp,
557 : golsen 1.8 '-j', 1, # function is profile alignment
558 :     '-B', $proffile, # location of profile
559 : golsen 1.3 '-d', $dbfile,
560 : golsen 1.8 '-i', $qfile,
561 :     '-e', $e_val,
562 : golsen 1.3 '-b', $nkeep,
563 :     '-v', $nkeep,
564 : golsen 1.8 '-t', 1 # issues warning if not set
565 : golsen 1.3 );
566 : golsen 1.8 push @cmd, ( '-a', $n_cpu ) if $n_cpu;
567 : golsen 1.3
568 : golsen 1.8 my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
569 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
570 : golsen 1.3 or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
571 :     and return undef;
572 :     my $out = &gjoparseblast::structured_blast_output( $blastfh, 1 ); # selfmatches okay
573 :     close $blastfh;
574 : overbeek 1.1
575 :     if ( $rm_db )
576 :     {
577 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
578 :     unlink @files if @files;
579 :     }
580 :     unlink $proffile if $rm_profile;
581 :     unlink $qfile if $rm_query;
582 :    
583 :     return $out;
584 :     }
585 :    
586 :    
587 :     #-------------------------------------------------------------------------------
588 : golsen 1.8 # Do psiblast against tranlated genomic DNA
589 :     #
590 :     # $structured_blast = blastpgpn( $dbfile, $profilefile, \%options )
591 :     # $structured_blast = blastpgpn( \@dbseq, $profilefile, \%options )
592 :     # $structured_blast = blastpgpn( $dbfile, \@profileseq, \%options )
593 :     # $structured_blast = blastpgpn( \@dbseq, \@profileseq, \%options )
594 :     #
595 :     # Required:
596 :     #
597 :     # $db_file or \@db_seq
598 :     # $profile_file or \@profile_seq
599 :     #
600 :     # Options:
601 :     #
602 :     # aa_db_file => $trans_file # put translation db here
603 :     # e_value => $max_e_value # D = 0.01
604 :     # max_e_val => $max_e_value # D = 0.01
605 :     # n_cpu => $n_cpu # synonym of n_thread
606 :     # n_result => $max_seq # D = 1000
607 :     # n_thread => $n_cpu # D = 1
608 :     # query => $q_file_name # most complete
609 :     # query => \@q_seq_entry # most complete
610 :     # stderr => $file # place to send blast stderr (D = /dev/null)
611 :     # tmp => $temp_dir # location for temp files
612 :     #
613 : fangfang 1.13 # If not supplied, $opts{ query } will be set to the most compete sequence.
614 :     #
615 : golsen 1.8 # 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.13 $opts->{ query } ||= $query;
712 :    
713 : golsen 1.8 my $e_val = $opts->{ e_value } || $opts->{ max_e_val } || $opts->{ max_e_value } || 0.01;
714 :     my $n_cpu = $opts->{ n_thread } || $opts->{ nthread } || $opts->{ n_cpu } || $opts->{ ncpu } || 2;
715 :     my $nkeep = $opts->{ n_result } || $opts->{ nresult } || 1000;
716 :    
717 : golsen 1.9 my $blastpgp = SeedAware::executable_for( 'blastpgp' )
718 :     or print STDERR "Could not find executable for program 'blastpgp'.\n"
719 :     and return undef;
720 :    
721 :     my @cmd = ( $blastpgp,
722 : golsen 1.8 '-j', 1, # function is profile alignment
723 :     '-B', $proffile, # location of protein profile
724 :     '-d', $dbfile, # protein database
725 :     '-i', $qfile, # one protein from the profile
726 :     '-e', $e_val,
727 :     '-b', $nkeep,
728 :     '-v', $nkeep,
729 :     '-t', 1, # issues warning if not set
730 :     );
731 :     push @cmd, ( '-a', $n_cpu ) if $n_cpu;
732 :    
733 :     my $opts2 = { stderr => ( $opts->{stderr} || '/dev/null' ) };
734 :     my $blastfh = SeedAware::read_from_pipe_with_redirect( @cmd, $opts2 )
735 :     or print STDERR "Failed to open: '" . join( ' ', @cmd ), "'.\n"
736 :     and return undef;
737 :     my $out = &gjoparseblast::structured_blast_output( $blastfh );
738 :     close $blastfh;
739 :    
740 :     if ( $rm_db )
741 :     {
742 :     my @files = grep { -f $_ } map { ( $_, "$_.psq", "$_.pin", "$_.phr" ) } $dbfile;
743 :     unlink @files if @files;
744 :     }
745 :     unlink $proffile if $rm_profile;
746 :     unlink $qfile if $rm_query;
747 :    
748 :     # Fix the blastp output to look like tblastn output
749 :    
750 :     foreach my $qdata ( @$out )
751 :     {
752 :     my %sdata; # There are now multiple sequences per subject DNA
753 :     foreach my $sdata ( @{ $qdata->[3] } )
754 :     {
755 :     my ( $sid, $sdef, undef, $hsps ) = @$sdata;
756 :     my $fr;
757 :     ( $sid, $fr ) = $sid =~ m/^(.*)\.([-+]\d)$/;
758 :     my ( $b, $e, $slen ) = $sdef =~ m/(\d+)-(\d+)\/(\d+)$/;
759 :     $sdef =~ s/ ?\S+$//;
760 :     my $sd2 = $sdata{ $sid } ||= [ $sid, $sdef, $slen, [] ];
761 :     foreach ( @$hsps ) { adjust_hsp( $_, $fr, $b ); push @{$sd2->[3]}, $_ }
762 :     }
763 :    
764 :     # Order the hsps for each subject
765 :    
766 :     foreach my $sid ( keys %sdata )
767 :     {
768 :     my $hsps = $sdata{ $sid }->[3];
769 :     @$hsps = sort { $b->[0] <=> $a->[0] } @$hsps;
770 :     }
771 :    
772 :     # Order the subjects for the query
773 :    
774 :     @{$qdata->[3]} = map { $_->[0] } # remove score
775 :     sort { $b->[1] <=> $a->[1] } # sort by score
776 :     map { [ $_, $_->[3]->[0]->[0] ] } # tag with top score
777 :     map { $sdata{ $_ } } # subject data
778 :     keys %sdata; # subject ids
779 :     }
780 :    
781 :     return $out;
782 :     }
783 :    
784 :    
785 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
786 :     # When search is versus six frame translation, there is a need to adjust
787 :     # the frame and location information in the hsp back to the DNA coordinates.
788 :     #
789 :     # adjust_hsp( $hsp, $frame, $begin )
790 :     #
791 :     # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
792 :     # scr Eval nseg Eval naln nid npos ngap fr q1 q2 qseq s1 s2 sseq
793 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
794 :     sub adjust_hsp
795 :     {
796 :     my ( $hsp, $fr, $b ) = @_;
797 :     $hsp->[8] = $fr;
798 :     if ( $fr > 0 )
799 :     {
800 :     $hsp->[12] = $b + 3 * ( $hsp->[12] - 1 );
801 :     $hsp->[13] = $b + 3 * ( $hsp->[13] - 1 ) + 2;
802 :     }
803 :     else
804 :     {
805 :     $hsp->[12] = $b - 3 * ( $hsp->[12] - 1 );
806 :     $hsp->[13] = $b - 3 * ( $hsp->[13] - 1 ) - 2;
807 :     }
808 :     }
809 :    
810 :    
811 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
812 :     # Do a six frame translation for use by blastpgpn. This modifications of
813 :     # the identifiers and definitions are essential to the interpretation of
814 :     # the blast results. The program 'translate_fasta_6' produces the same
815 :     # output format, and is much faster.
816 :     #
817 :     # @translations = six_translations( $nucleotide_entry )
818 :     #
819 :     # The ids are modified by adding ".frame" (+1, +2, +3, -1, -2, -3).
820 :     # The definition is mofidified by adding " begin-end/of_length".
821 :     # NCBI frame numbers reverse strand translation frames from the end of the
822 :     # sequence (i.e., the beginning of the complement of the strand).
823 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
824 :    
825 :     sub six_translations
826 :     {
827 :     my ( $id, $def, $seq ) = map { defined($_) ? $_ : '' } @{ $_[0] };
828 :     my $l = length( $seq );
829 :    
830 :     return () if $l < 15;
831 :    
832 :     # fr beg end
833 :     my @intervals = ( [ '+1', 1, $l - ( $l % 3 ) ],
834 :     [ '+2', 2, $l - ( ($l-1) % 3 ) ],
835 :     [ '+3', 3, $l - ( ($l-2) % 3 ) ],
836 :     [ '-1', $l, 1 + ( $l % 3 ) ],
837 :     [ '-2', $l-1, 1 + ( ($l-1) % 3 ) ],
838 :     [ '-3', $l-2, 1 + ( ($l-2) % 3 ) ]
839 :     );
840 :     my ( $fr, $b, $e );
841 :    
842 :     map { ( $fr, $b, $e ) = @$_;
843 :     [ "$id.$fr", "$def $b-$e/$l", gjoseqlib::translate_seq( gjoseqlib::DNA_subseq( \$seq, $b, $e ) ) ]
844 :     } @intervals;
845 :     }
846 :    
847 :    
848 :     #-------------------------------------------------------------------------------
849 : golsen 1.3 # Get an input file handle, and boolean on whether to close or not:
850 :     #
851 :     # ( \*FH, $close ) = input_file_handle( $filename );
852 :     # ( \*FH, $close ) = input_file_handle( \*FH );
853 :     # ( \*FH, $close ) = input_file_handle( ); # D = STDIN
854 :     #
855 :     #-------------------------------------------------------------------------------
856 :    
857 :     sub input_file_handle
858 :     {
859 :     my ( $file, $umask ) = @_;
860 :    
861 :     my ( $fh, $close );
862 :     if ( defined $file )
863 :     {
864 :     if ( ref $file eq 'GLOB' )
865 :     {
866 :     $fh = $file;
867 :     $close = 0;
868 :     }
869 :     elsif ( -f $file )
870 :     {
871 :     open( $fh, "<$file") || die "input_file_handle could not open '$file'.\n";
872 :     $close = 1;
873 :     }
874 :     else
875 :     {
876 :     die "input_file_handle could not find file '$file'.\n";
877 :     }
878 :     }
879 :     else
880 :     {
881 :     $fh = \*STDIN;
882 :     $close = 0;
883 :     }
884 :    
885 :     return ( $fh, $close );
886 :     }
887 :    
888 :    
889 :     #-------------------------------------------------------------------------------
890 : overbeek 1.1 # Get an output file handle, and boolean on whether to close or not:
891 :     #
892 :     # ( \*FH, $close ) = output_file_handle( $filename );
893 :     # ( \*FH, $close ) = output_file_handle( \*FH );
894 :     # ( \*FH, $close ) = output_file_handle( ); # D = STDOUT
895 :     #
896 :     #-------------------------------------------------------------------------------
897 :    
898 :     sub output_file_handle
899 :     {
900 :     my ( $file, $umask ) = @_;
901 :    
902 :     my ( $fh, $close );
903 :     if ( defined $file )
904 :     {
905 :     if ( ref $file eq 'GLOB' )
906 :     {
907 :     $fh = $file;
908 :     $close = 0;
909 :     }
910 :     else
911 :     {
912 :     open( $fh, ">$file") || die "output_file_handle could not open '$file'.\n";
913 :     chmod 0664, $file; # Seems to work on open file!
914 :     $close = 1;
915 :     }
916 :     }
917 :     else
918 :     {
919 :     $fh = \*STDOUT;
920 :     $close = 0;
921 :     }
922 :    
923 :     return ( $fh, $close );
924 :     }
925 :    
926 :    
927 :     #-------------------------------------------------------------------------------
928 :     #
929 :     # print_blast_as_records( $file, \@queries )
930 :     # print_blast_as_records( \*FH, \@queries )
931 :     # print_blast_as_records( \@queries ) # STDOUT
932 :     #
933 :     # \@queries = read_blast_from_records( $file )
934 :     # \@queries = read_blast_from_records( \*FH )
935 :     # \@queries = read_blast_from_records( ) # STDIN
936 :     #
937 :     # Three output record types:
938 :     #
939 :     # Query \t qid \t qdef \t dlen
940 :     # > \t sid \t sdef \t slen
941 :     # HSP \t ...
942 :     #
943 :     #-------------------------------------------------------------------------------
944 :     sub print_blast_as_records
945 :     {
946 :     my $file = ( $_[0] && ! ( ref $_[0] eq 'ARRAY' ) ? shift : undef );
947 :     my ( $fh, $close ) = output_file_handle( $file );
948 :    
949 :     my $queries = shift;
950 :     $queries && ref $queries eq 'ARRAY'
951 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
952 :     and return;
953 :    
954 :     foreach my $qdata ( @$queries )
955 :     {
956 :     $qdata && ref $qdata eq 'ARRAY' && @$qdata == 4
957 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
958 :     and return;
959 :    
960 :     my $subjcts = pop @$qdata;
961 :     $subjcts && ref $subjcts eq 'ARRAY'
962 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
963 :     and return;
964 :    
965 :     print $fh join( "\t", 'Query', @$qdata ), "\n";
966 :     foreach my $sdata ( @$subjcts )
967 :     {
968 :     $sdata && ref $sdata eq 'ARRAY' && @$sdata == 4
969 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
970 :     and return;
971 :    
972 :     my $hsps = pop @$sdata;
973 :     $hsps && ref $hsps eq 'ARRAY' && @$hsps
974 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
975 :     and return;
976 :    
977 :     print $fh join( "\t", '>', @$sdata ), "\n";
978 :    
979 :     foreach my $hsp ( @$hsps )
980 :     {
981 :     $hsp && ref $hsp eq 'ARRAY' && @$hsp > 4
982 :     or print STDERR "Bad blast data supplied to print_blast_as_records.\n"
983 :     and return;
984 :     print $fh join( "\t", 'HSP', @$hsp ), "\n";
985 :     }
986 :     }
987 :     }
988 :    
989 :     close $fh if $close;
990 :     }
991 :    
992 :    
993 :     sub read_blast_from_records
994 :     {
995 :     my ( $file ) = @_;
996 : golsen 1.3
997 :     my ( $fh, $close ) = input_file_handle( $file );
998 :     $fh or print STDERR "read_blast_from_records could not open input file.\n"
999 : overbeek 1.1 and return wantarray ? () : [];
1000 : golsen 1.3
1001 : overbeek 1.1 my @queries = ();
1002 :     my @qdata;
1003 :     my @subjects = ();
1004 :     my @sdata;
1005 :     my @hsps = ();
1006 :    
1007 :     local $_;
1008 :     while ( defined( $_ = <$fh> ) )
1009 :     {
1010 :     chomp;
1011 :     my ( $type, @datum ) = split /\t/;
1012 :     if ( $type eq 'Query' )
1013 :     {
1014 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1015 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
1016 :     @qdata = @datum;
1017 :     @sdata = ();
1018 :     @subjects = ();
1019 :     @hsps = ();
1020 :     }
1021 :     elsif ( $type eq '>' )
1022 :     {
1023 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1024 :     @sdata = @datum;
1025 :     @hsps = ();
1026 :     }
1027 :     elsif ( $type eq 'HSP' )
1028 :     {
1029 :     push @hsps, \@datum;
1030 :     }
1031 :     }
1032 :    
1033 :     close $fh if $close;
1034 :    
1035 :     push @subjects, [ @sdata, [ @hsps ] ] if @hsps;
1036 :     push @queries, [ @qdata, [ @subjects ] ] if @qdata;
1037 :    
1038 :     wantarray ? @queries : \@queries;
1039 :     }
1040 :    
1041 :    
1042 : golsen 1.3 sub verify_db
1043 :     {
1044 :     my ( $db, $type ) = @_;
1045 :    
1046 :     my @args;
1047 : golsen 1.8 if ( $type =~ m/^p/i )
1048 : golsen 1.3 {
1049 : golsen 1.8 @args = ( "-p", "T", "-i", $db ) if ((! -s "$db.psq") || (-M "$db.psq" > -M $db));
1050 : golsen 1.3 }
1051 :     else
1052 :     {
1053 : golsen 1.8 @args = ( "-p", "F", "-i", $db ) if ((! -s "$db.nsq") || (-M "$db.nsq" > -M $db));
1054 : golsen 1.3 }
1055 : golsen 1.8 @args or return ( -s $db ? 1 : 0 );
1056 :    
1057 : golsen 1.3 #
1058 : golsen 1.8 # Find formatdb appropriate for the excecution environemnt.
1059 : golsen 1.3 #
1060 :    
1061 :     my $prog = SeedAware::executable_for( 'formatdb' );
1062 :     if ( ! $prog )
1063 :     {
1064 :     warn "gjoalignandtree::verify_db: formatdb program not found.\n";
1065 :     return 0;
1066 :     }
1067 :    
1068 :     my $rc = system( $prog, @args );
1069 :    
1070 :     if ( $rc != 0 )
1071 :     {
1072 :     my $cmd = join( ' ', $prog, @args );
1073 :     warn "gjoalignandtree::verify_db: formatdb failed with rc = $rc: $cmd\n";
1074 :     return 0;
1075 :     }
1076 :    
1077 :     return 1;
1078 :     }
1079 :    
1080 :    
1081 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3