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

Annotation of /FigKernelPackages/gjophylip.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : fangfang 1.6 # This is a SAS component
2 :     #
3 :     # Copyright (c) 2003-2010 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     #
18 :    
19 : golsen 1.1 package gjophylip;
20 :    
21 :     #===============================================================================
22 :     # A perl interface to programs in the PHYLIP program package:
23 :     #===============================================================================
24 :     # proml:
25 :     #
26 :     # @tree_likelihood_pairs = proml( \@alignment, \%options )
27 :     # @tree_likelihood_pairs = proml( \@alignment, %options )
28 : golsen 1.3 # @tree_likelihood_pairs = proml( \%options )
29 :     # @tree_likelihood_pairs = proml( %options )
30 : golsen 1.1 #
31 : golsen 1.3 # estimate_protein_site_rates -- Use proml to estimate site-specific rates of change
32 :     #
33 :     # ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, \%proml_opts )
34 :     #
35 : golsen 1.4 # optimize_protein_coef_of_var:
36 :     #
37 :     # ( $coef, $tree, $score ) = optimize_protein_coef_of_var( \@alignment, \%options )
38 :     #
39 : golsen 1.3 # dnadist:
40 :     #
41 :     # $dists = dnadist( \@alignment, \%options )
42 :     # $dists = dnadist( %options )
43 :     #
44 :     # dnadist and neighbor:
45 :     #
46 :     # $tree = dnadist_neighbor( \@alignment, \%options )
47 :     # $tree = dnadist_neighbor( \@alignment, %options )
48 :     # $tree = dnadist_neighbor( \%options )
49 :     # $tree = dnadist_neighbor( %options )
50 : golsen 1.1 #
51 :     # protdist:
52 :     #
53 :     # $distance_matrix = protdist( \@alignment, \%options )
54 : golsen 1.3 # $distance_matrix = protdist( %options )
55 : golsen 1.1 #
56 :     # protdist and neighbor:
57 :     #
58 :     # $tree = protdist_neighbor( \@alignment, \%options )
59 :     # $tree = protdist_neighbor( \@alignment, %options )
60 : golsen 1.3 # $tree = protdist_neighbor( \%options )
61 :     # $tree = protdist_neighbor( %options )
62 : golsen 1.1 #
63 : golsen 1.3 # neighbor:
64 : golsen 1.1 #
65 : golsen 1.3 # $tree = neighbor( \@distances, \%options )
66 :     # $tree = neighbor( \@distances, %options )
67 :     # $tree = neighbor( \%options )
68 :     # $tree = neighbor( %options )
69 : golsen 1.1 #
70 : golsen 1.3 # consense:
71 : golsen 1.2 #
72 :     # $consensus_with_boot_as_branch_len = consensus_tree( @trees, \%options );
73 :     # $consensus_with_boot_as_branch_len = consensus_tree( \@trees, \%options );
74 :     # $consensus_with_boot_as_branch_len = consensus_tree( @trees );
75 :     # $consensus_with_boot_as_branch_len = consensus_tree( \@trees );
76 :     #
77 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, @trees, \%options );
78 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, \@trees, \%options );
79 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, @trees );
80 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, \@trees );
81 :     #
82 : golsen 1.3 #
83 :     # Common options for PHYLIP programs:
84 :     #
85 :     # alignment => \@alignment supply the alignment as an option
86 :     # alpha => float alpha param of gamma distrib. (0.5 - inf)
87 :     # categories => $categories category rates and site-specific rate categs
88 :     # coef_of_var => float 1/sqrt(alpha) of gamma distrib. (D = 0)
89 :     # gamma_bins => int rates in approx. gamma distrib. (D = 5)
90 :     # global => bool global rearrangements
91 :     # invar_frac => 0 - 1 fraction of site that are invariant (D = 0)
92 :     # jumble_seed => odd_int jumble random seed
93 :     # model => model protein evolution model: JTT (D) | PMB | PAM
94 :     # n_jumble => int number of order of addition jumbles
95 :     # persistance => float persistance length of rate category (D = 0)
96 :     # rate_hmm => $rate_hmm # not implimented
97 :     # rearrange => \@trees rearrange user trees
98 :     # slow => bool more accurate but slower search (D = 0)
99 :     # user_lengths => bool use user-supplied branch lengths (D = 0)
100 :     # user_trees => \@trees use user-supplied trees
101 :     # weights => $site_weights site-specific weights
102 :     #
103 :     #
104 :     # Data types:
105 :     #
106 :     # @alignment = ( [ id, def, seq ], [ id, def, seq ] ... )
107 :     #
108 :     # $categories = [ [ cat1_rate, cat2_rate, ... ], $category_per_site_string ]
109 :     #
110 :     # $rate_hmm = [ [ rate_1, prob_1 ], [ rate_2, prob2 ], ... ]
111 :     #
112 :     # $weights = $weight_per_site_string
113 :     #
114 :     #
115 :     # Common options for programming inferface:
116 :     #
117 :     # keep_duplicates => bool do not remove duplicate sequences [NOT IMPLIMENTED]
118 :     # tmp => directory directory for tmp_dir
119 :     # tmp_dir => directory directory for temporary files
120 :     # tree_format => format format of output tree: overbeek | gjo | fig
121 :     #
122 :     # Options that permit defining path or nonstandard name for program:
123 :     #
124 :     # dnadist => dnadist just of dnadist_neighbor
125 :     # program => program program name, possibly with full path
126 :     # protdist => protdist just of protdist_neighbor
127 :     # neighbor => neighbor dnadist_neighbor and protdist_neighbor
128 :     #
129 :     #
130 :     # Auxilliary functions:
131 :     #
132 :     # [ \@rates, $site_categs ] = read_categories( $categories_file )
133 :     #
134 :     # $site_weights = read_weights( $weights_file )
135 :     #
136 :     # @distance_matrix = read_distances() # D = outfile
137 :     # \@distance_matrix = read_distances()
138 :     # @distance_matrix = read_distances( $file )
139 :     # \@distance_matrix = read_distances( $file )
140 :     #
141 :     # ( \%id_index, \%bs_frac ) = read_consense_partitians() # D = outfile
142 :     # ( \%id_index, \%bs_frac ) = read_consense_partitians( $file )
143 :     #
144 : golsen 1.1 #===============================================================================
145 :    
146 :     use strict;
147 : golsen 1.3 use SeedAware;
148 : golsen 1.1 use gjonewicklib qw( gjonewick_to_overbeek
149 :     newick_is_unrooted
150 :     newick_relabel_nodes
151 :     newick_rescale_branches
152 :     newick_tree_length
153 :     overbeek_to_gjonewick
154 :     parse_newick_tree_str
155 :     strNewickTree
156 :     uproot_newick
157 :     );
158 : golsen 1.3 use Data::Dumper;
159 : golsen 1.1
160 :    
161 :     #-------------------------------------------------------------------------------
162 :     # proml -- A perl interface to the proml program
163 :     #
164 :     # @tree_likelihood_pairs = proml( \@alignment, \%options )
165 :     # @tree_likelihood_pairs = proml( \@alignment, %options )
166 :     # @tree_likelihood_pairs = proml( \%options )
167 :     # @tree_likelihood_pairs = proml( %options )
168 :     #
169 :     # @alignment = ( [ id, seq ], [ id, seq ] ... )
170 :     # or
171 :     # ( [ id, def, seq ], [ id, def, seq ] ... )
172 :     #
173 :     # options:
174 :     #
175 :     # proml:
176 :     # alignment => \@alignment supply the alignment as an option
177 :     # alpha => float alpha parameter of gamma distribution (0.5 - inf)
178 :     # categories => [ [ rate1, ... ], site_categories ]
179 :     # coef_of_var => float 1/sqrt(alpha) for gamma distribution (D = 0)
180 :     # gamma_bins => int number of rate categories used to approximate gamma (D=5)
181 :     # global => bool global rearrangements
182 :     # invar_frac => 0 - 1 fraction of site that are invariant
183 :     # jumble_seed => odd int jumble random seed
184 :     # model => model evolution model JTT (D) | PMB | PAM
185 :     # n_jumble => int number of jumbles
186 :     # persistance => float persistance length of rate category
187 : golsen 1.3 # rate_hmm => rate_prob_pairs general hidden markov model with 2 to 9 catefories
188 : golsen 1.1 # rearrange => [ trees ] rearrange user trees
189 :     # slow => bool more accurate but slower search (D = 0)
190 :     # user_lengths => bool use supplied branch lengths
191 :     # user_trees => [ trees ] user trees
192 :     # weights => site_weights
193 :     #
194 :     # other:
195 : golsen 1.3 # definitions => bool include definition in sequence label
196 : golsen 1.1 # keep_duplicates => bool do not remove duplicate sequences [NOT IMPLIMENTED]
197 : golsen 1.3 # no_ids => bool leave the sequence id out of the sequence label
198 : golsen 1.1 # program => program allows fully defined path
199 :     # tmp => directory directory for tmp_dir
200 :     # tmp_dir => directory directory for temporary files
201 :     # tree_format => overbeek | gjo | fig format of output tree
202 :     #
203 :     # tmp_dir is created and deleted unless its name is supplied, and it already
204 :     # exists.
205 :     #
206 :     #-------------------------------------------------------------------------------
207 :    
208 :     sub proml
209 :     {
210 :     my $align;
211 : golsen 1.3 $align = shift @_ if ( ref( $_[0] ) eq 'ARRAY' );
212 :    
213 :     # Collect the options in a hash, and canonicalize their keys:
214 : golsen 1.1
215 :     my %options;
216 :     if ( $_[0] )
217 :     {
218 :     my %opts = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
219 :     %options = map { canonical_key( $_ ) => $opts{ $_ } } keys %opts;
220 :     }
221 :    
222 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
223 : golsen 1.3 # Clean up a local copy of the alignment. Assign local ids
224 :     # Alignment can also be an option:
225 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226 :    
227 :     $align ||= $options{ alignment } || $options{ align };
228 :    
229 : golsen 1.3 my ( $align2, $id, $local_id ) = &process_protein_alignment( $align );
230 :    
231 :     if ( ! $align2 )
232 :     {
233 :     print STDERR " gjophylip::proml requires an alignment.\n";
234 :     return ();
235 :     }
236 : golsen 1.1
237 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238 :     # Process proml options:
239 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240 :    
241 :     my $categories = $options{ categorie }; # The 's' is lost in canonicalizing keys
242 :     if ( $categories )
243 :     {
244 : golsen 1.3 ( $categories = &process_categories( $categories ) )
245 :     or print STDERR " gjophylip::proml bad categories option value.\n"
246 :     and return ();
247 : golsen 1.1 }
248 :    
249 :     my $coef_of_var = $options{ coefofvar }
250 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
251 :     || 0;
252 :     if ( $coef_of_var < 0 )
253 :     {
254 :     print STDERR "gjophylip::proml coef_of_var option value must be >= 0\n";
255 :     return ();
256 :     }
257 :    
258 : golsen 1.3 my $no_ids = $options{ noid } || 0;
259 :     my $definitions = $options{ definition } || $no_ids; # No id requires definition
260 :     if ( $definitions )
261 :     {
262 :     my %new_id = map { $_->[0] => ( ! length $_->[1] ? $_->[0] : $no_ids ? $_->[1] : "$_->[0] $_->[1]" ) }
263 :     map { [ $_->[0], @$_ == 3 && defined $_->[1] ? $_->[1] : '' ] }
264 :     @$align;
265 :     foreach ( keys %$id ) { $id->{ $_ } = $new_id{ $id->{ $_ } } }
266 :     }
267 :    
268 : golsen 1.1 my $gamma_bins = int( $options{ gammabin } || ( $coef_of_var ? 5 : 2 ) );
269 :     if ( ( $gamma_bins < 2 ) || ( $gamma_bins > 9 ) )
270 :     {
271 :     print STDERR "gjophylip::proml gamma_bins option value must be > 1 and <= 9\n";
272 :     return ();
273 :     }
274 :    
275 :     my $global = $options{ global } || 0;
276 :    
277 :     my $invar_frac = $options{ invarfrac } || 0;
278 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
279 :     {
280 :     print STDERR "gjophylip::proml invar_frac option value must be >= 0 and < 1\n";
281 :     return ();
282 :     }
283 :    
284 :     my $n_jumble = int( $options{ njumble } || ( $options{ jumbleseed } ? 1 : 0) );
285 :     if ( $n_jumble < 0 )
286 :     {
287 :     print STDERR "gjophylip::proml n_jumble option value must be >= 0\n";
288 :     return ();
289 :     }
290 :    
291 :     my $jumble_seed = int( $options{ jumbleseed } || 4 * int( 499999999 * rand() ) + 1 );
292 :     if ( ( $jumble_seed <= 0) || ( $jumble_seed % 2 != 1 ) )
293 :     {
294 :     print STDERR "gjophylip::proml jumble_seed option value must be an odd number > 0\n";
295 :     return ();
296 :     }
297 :    
298 : golsen 1.3 my $model = &process_protein_model( \%options );
299 : golsen 1.1
300 : golsen 1.3 my $persistance = $options{ persistance } || 0;
301 : golsen 1.1 if ( $persistance && ( $persistance <= 1 ) )
302 :     {
303 :     print STDERR "gjophylip::proml persistance option value must be > 1\n";
304 :     return ();
305 :     }
306 :    
307 : golsen 1.3 my $rate_hmm = $options{ ratehmm };
308 :     if ( $rate_hmm )
309 :     {
310 :     ( $rate_hmm = &process_rate_hmm( $rate_hmm ) )
311 :     or print STDERR " gjophylip::proml bad rate_hmm value\n"
312 :     and return ();
313 :     }
314 : golsen 1.1
315 : golsen 1.3 my $rearrange = $options{ rearrange };
316 :    
317 :     my $slow = $options{ slow };
318 : golsen 1.1
319 :     my $user_lengths = $options{ userlength };
320 :    
321 : golsen 1.3 my $user_trees = $options{ usertree } || $rearrange;
322 : golsen 1.1 if ( $user_trees )
323 :     {
324 :     if ( ( ref( $user_trees ) ne 'ARRAY' ) || ( ! @$user_trees ) )
325 :     {
326 :     $user_trees = undef; # No trees
327 :     }
328 :     elsif ( ref( $user_trees->[0] ) ne 'ARRAY' ) # First element not tree
329 :     {
330 :     print STDERR "gjophylip::proml user_trees or rearrange option value must be reference to list of trees\n";
331 :     return ();
332 :     }
333 :     }
334 :    
335 :     my $weights = $options{ weight };
336 :    
337 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
338 :     # Options that are not proml options per se:
339 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
340 :    
341 : golsen 1.3 my $program = SeedAware::executable_for( $options{ proml } || $options{ program } || 'proml' );
342 : golsen 1.1
343 : golsen 1.3 my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( \%options );
344 : golsen 1.1 $tmp_dir or return ();
345 :    
346 : golsen 1.3 my $tree_format = process_tree_format( \%options );
347 : golsen 1.1
348 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349 :     # Prepare data:
350 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351 :    
352 :     my @user_trees = ();
353 :     if ( $user_trees && @$user_trees )
354 :     {
355 : golsen 1.3 if ( gjonewicklib::is_overbeek_tree( @$user_trees[0] ) )
356 : golsen 1.1 {
357 :     @user_trees = map { gjonewicklib::overbeek_to_gjonewick( $_ ) }
358 :     @$user_trees;
359 :     }
360 :     else
361 :     {
362 :     @user_trees = map { gjonewicklib::copy_newick_tree( $_ ) }
363 :     @$user_trees;
364 :     }
365 :    
366 :     # Relabel and make sure trees are unrooted:
367 :    
368 :     @user_trees = map { gjonewicklib::newick_is_unrooted( $_ )
369 :     ? $_
370 :     : gjonewicklib::uproot_newick( $_ )
371 :     }
372 : golsen 1.3 map { gjonewicklib::newick_relabel_nodes( $_, $local_id ); $_ }
373 : golsen 1.1 @user_trees;
374 :     }
375 :    
376 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
377 :     # Write the files and run the program:
378 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
379 :    
380 :     my $cwd = $ENV{ cwd } || `pwd`;
381 :     chomp $cwd;
382 :     chdir $tmp_dir;
383 :    
384 : golsen 1.3 my ( $trees, $likelihoods ) = run_proml( align => $align2,
385 :     categories => $categories,
386 :     coef_of_var => $coef_of_var,
387 :     gamma_bins => $gamma_bins,
388 :     global => $global,
389 :     invar_frac => $invar_frac,
390 :     jumble_seed => $jumble_seed,
391 :     model => $model,
392 :     n_jumble => $n_jumble,
393 :     program => $program,
394 :     persistance => $persistance,
395 :     rate_hmm => $rate_hmm,
396 :     rearrange => $rearrange,
397 :     slow => $slow,
398 :     user_lengths => $user_lengths,
399 :     user_trees => \@user_trees,
400 :     weights => $weights,
401 : golsen 1.1 );
402 :    
403 :     # We are done, go back to the original directory:
404 :    
405 :     chdir $cwd;
406 :    
407 :     # Delete the temporary directory unless it already existed:
408 :    
409 : golsen 1.3 system( '/bin/rm', '-r', $tmp_dir ) if ! $save_tmp;
410 : golsen 1.1
411 :     # Fix the tree labels:
412 :    
413 : golsen 1.3 my @trees = map { gjonewicklib::newick_relabel_nodes( $_, $id ) }
414 : golsen 1.1 @$trees;
415 :    
416 :     # Fix the tree format, if necessary:
417 :    
418 :     if ( $tree_format =~ m/overbeek/i )
419 :     {
420 :     @trees = map { gjonewicklib::gjonewick_to_overbeek( $_ ) } @trees;
421 :     }
422 :    
423 :     # Merge into a list of tree-likelihood pairs:
424 :    
425 :     return map { [ $_, shift @$likelihoods ] } @trees;
426 :     }
427 :    
428 :    
429 :     #-------------------------------------------------------------------------------
430 :     # A local routine to run proml
431 :     #-------------------------------------------------------------------------------
432 :     sub run_proml
433 :     {
434 : golsen 1.3 my %data = @_;
435 : golsen 1.1
436 :     unlink 'outfile' if -f 'outfile'; # Just checking
437 :     unlink 'outtree' if -f 'outtree'; # ditto
438 :    
439 : golsen 1.3 &write_seq_infile( @{ $data{align} } )
440 : golsen 1.1 or print STDERR "gjophylip::run_proml: Could not write infile\n"
441 :     and return ();
442 :    
443 :     open( PROML, ">proml_cmd" )
444 : golsen 1.3 or print STDERR "gjophylip::run_proml: Could not open command file\n"
445 : golsen 1.1 and return ();
446 :    
447 : golsen 1.3 if ( $data{categories} )
448 : golsen 1.1 {
449 : golsen 1.3 &write_categories( $data{categories}->[1] )
450 : golsen 1.1 or print STDERR "gjophylip::run_proml: Could not write categories\n"
451 :     and return ();
452 :     print PROML "C\n",
453 : golsen 1.3 scalar @{$data{categories}->[0]}, "\n",
454 :     join( ' ', @{ $data{categories}->[0] } ), "\n";
455 : golsen 1.1 }
456 :    
457 : golsen 1.3 if ( $data{invar_frac} || $data{coef_of_var} || $data{rate_hmm} )
458 : golsen 1.1 {
459 :     print PROML "R\n";
460 : golsen 1.3 print PROML "R\n" if $data{invar_frac} || $data{rate_hmm};
461 :     print PROML "R\n" if $data{rate_hmm};
462 :     print PROML "A\n", "$data{persistance}\n" if $data{persistance};
463 : golsen 1.1
464 :     }
465 :    
466 : golsen 1.3 print PROML "G\n" if $data{global};
467 : golsen 1.1
468 : golsen 1.3 print PROML "J\n", "$data{jumble_seed}\n", "$data{n_jumble}\n" if $data{n_jumble};
469 : golsen 1.1
470 : golsen 1.3 print PROML "P\n" if $data{model} =~ m/PMB/i;
471 :     print PROML "P\nP\n" if $data{model} =~ m/PAM/i;
472 : golsen 1.1
473 : golsen 1.3 if ( @{$data{user_trees}} )
474 : golsen 1.1 {
475 : golsen 1.3 &write_intree( @{$data{user_trees}} )
476 : golsen 1.1 or print STDERR "gjophylip::run_proml: Could not write intree\n"
477 :     and return ();
478 :     print PROML "U\n";
479 : golsen 1.3 print PROML "V\n" if $data{rearrange} || $data{global};
480 :     print PROML "L\n" if $data{user_lengths} && ! $data{rearrange} && ! $data{global};
481 : golsen 1.1 }
482 : golsen 1.3 elsif ( $data{slow} ) # Slow and user trees are mutually exclusive
483 : golsen 1.1 {
484 :     print PROML "S\n";
485 :     }
486 :    
487 : golsen 1.3 if ( $data{weights} )
488 : golsen 1.1 {
489 : golsen 1.3 &write_weights( $data{weights} )
490 : golsen 1.1 or print STDERR "gjophylip::run_proml: Could not write weights\n"
491 :     and return ();
492 :     print PROML "W\n";
493 :     }
494 :    
495 :     # We have identified all the options; indicate this with 'Y':
496 :    
497 :     print PROML "Y\n";
498 :    
499 :     # Becuase of the options interface, some values still must be supplied:
500 :    
501 : golsen 1.3 if ( $data{invar_frac} || $data{coef_of_var} )
502 : golsen 1.1 {
503 : golsen 1.3 if ( $data{invar_frac} )
504 : golsen 1.1 {
505 : golsen 1.3 if ( $data{coef_of_var} ) { $data{gamma_bins}++ if $data{gamma_bins} < 9 }
506 :     else { $data{gamma_bins} = 2 }
507 : golsen 1.1 }
508 : golsen 1.3 print PROML "$data{coef_of_var}\n";
509 :     print PROML "$data{gamma_bins}\n";
510 :     print PROML "$data{invar_frac}\n" if $data{invar_frac};
511 :     }
512 :     elsif ( $data{rate_hmm} )
513 :     {
514 :     my @hmm = @{ $data{rate_hmm} };
515 :     print PROML scalar @hmm, "\n";
516 :     primt PROML join( ' ', map { sprintf( "%.6f", $_->[0] ) } @hmm ), "\n";
517 :     primt PROML join( ' ', map { sprintf( "%.6f", $_->[1] ) } @hmm ), "\n";
518 : golsen 1.1 }
519 :    
520 : golsen 1.3 if ( @{$data{user_trees}} )
521 : golsen 1.1 {
522 :     print PROML "13\n"; # Random number seed of unknown use
523 :     }
524 :    
525 :     close PROML;
526 :    
527 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
528 :     # Run the program
529 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
530 :    
531 : golsen 1.3 my $program = $data{ program } || $data{ proml } || 'proml';
532 :     my $redirects = { stdin => 'proml_cmd',
533 :     stdout => 'proml_out', # Gets deleted with tmp dir
534 :     stderr => 'proml_err' # Gets deleted with tmp dir
535 :     };
536 :     SeedAware::system_with_redirect( $program, $redirects )
537 :     and print STDERR "gjophylip::run_proml: Could not run $program.\n"
538 :     and return undef;
539 : golsen 1.1
540 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
541 :     # Gather the results
542 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
543 :    
544 :     my @likelihoods = &read_likelihoods();
545 :    
546 :     my @trees = gjonewicklib::read_newick_trees( 'outtree' );
547 :     @trees or print STDERR "gjophylip::run_proml: Could not read proml outtree file\n"
548 :     and return ();
549 :    
550 :     return ( \@trees, \@likelihoods );
551 :     }
552 :    
553 :    
554 :     #-------------------------------------------------------------------------------
555 :     # estimate_protein_site_rates -- Use proml to estimate site-specific rates of change
556 :     #
557 :     # ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, %proml_opts )
558 :     # ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, \%proml_opts )
559 :     #
560 :     # $categories = [ [ $rate1, ... ], $site_categories ];
561 :     #
562 :     # $alignment = [ [ id, def, seq ], ... ]
563 :     # or
564 :     # [ [ id, seq ], ... ]
565 :     #
566 :     # proml_opts is list of key value pairs, or reference to a hash
567 :     #-------------------------------------------------------------------------------
568 :    
569 :     sub estimate_protein_site_rates
570 :     {
571 :     my ( $align, $tree, @proml_opts ) = @_;
572 :    
573 : golsen 1.3 # Clean-up the alignment, and make local ids that are compatible with
574 :     # the PHYLIP programs:
575 :    
576 :     my ( $id, $local_id );
577 :     ( $align, $id, $local_id ) = &process_protein_alignment( $align );
578 :    
579 :     if ( ! $align )
580 :     {
581 :     print STDERR " gjophylip::estimate_protein_site_rates requires an alignment.\n";
582 :     return ();
583 :     }
584 :    
585 :     my @align = @$align;
586 : golsen 1.1
587 :     # Make the tree a gjonewick tree, uproot it, and change to the local ids.
588 :    
589 : golsen 1.3 if ( gjonewicklib::is_overbeek_tree( $tree ) )
590 : golsen 1.1 {
591 :     $tree = gjonewicklib::overbeek_to_gjonewick( $tree );
592 :     }
593 :     else
594 :     {
595 :     $tree = gjonewicklib::copy_newick_tree( $tree );
596 :     }
597 :    
598 :     $tree = gjonewicklib::uproot_newick( $tree ) if ! gjonewicklib::newick_is_unrooted( $tree );
599 :    
600 : golsen 1.3 gjonewicklib::newick_relabel_nodes( $tree, $local_id );
601 : golsen 1.1
602 :     # The minimum rate will be 1/2 change per total tree branch length.
603 :     # This needs to be checked for proml. The intent is that he optimal
604 :     # rate for a site with one amino acid change is twice this value.
605 :    
606 :     my $kmin = 1 / ( gjonewicklib::newick_tree_length( $tree ) || 1 );
607 :    
608 :     # Generate "rate variation" by rescaling the supplied tree. We could use a
609 :     # finer grain estimator, then categorize the inferred values. This might
610 :     # work slightly better (this is what DNArates currently does).
611 :    
612 :     my $f = exp( log( 2 ) / 1 ); # Interval of 2
613 :     my @rates = map { $kmin * $f**$_ } ( 0 .. 16 ); # kmin .. 65000 * kmin in 17 bins
614 :     my @cat_vals = ( 1 .. 17 );
615 :     my @trees;
616 :     my $rate;
617 :     foreach $rate ( @rates )
618 :     {
619 :     my $tr = gjonewicklib::copy_newick_tree( $tree );
620 :     gjonewicklib::newick_rescale_branches( $tr, $rate ); # Rescales in place
621 :     push @trees, $tr;
622 :     }
623 :    
624 :     # Adjust (a copy of) the proml opts:
625 :    
626 :     my %proml_opts = ( ref( $proml_opts[0] ) eq 'HASH' ) ? %{ $proml_opts[0] } : @proml_opts;
627 :    
628 : golsen 1.3 $proml_opts{ user_lengths } = 1;
629 : golsen 1.1 $proml_opts{ user_trees } = \@trees;
630 : golsen 1.3 $proml_opts{ tree_format } = 'gjonewick';
631 : golsen 1.1
632 :     delete $proml_opts{ alpha } if exists $proml_opts{ alpha };
633 :     delete $proml_opts{ categories } if exists $proml_opts{ categories };
634 :     delete $proml_opts{ coef_of_var } if exists $proml_opts{ coef_of_var };
635 :     delete $proml_opts{ gamma_bins } if exists $proml_opts{ gamma_bins };
636 :     delete $proml_opts{ invar_frac } if exists $proml_opts{ invar_frac };
637 :     delete $proml_opts{ jumble_seed } if exists $proml_opts{ jumble_seed };
638 :     delete $proml_opts{ n_jumble } if exists $proml_opts{ n_jumble };
639 : golsen 1.3 delete $proml_opts{ rate_hmm } if exists $proml_opts{ rate_hmm };
640 : golsen 1.1 delete $proml_opts{ rearrange } if exists $proml_opts{ rearrange };
641 :    
642 :     # Work throught the sites, finding their optimal rates/categories:
643 :    
644 :     my @categories;
645 :     my @weights;
646 :     my $imax = length( $align[0]->[-1] );
647 :     for ( my $i = 0; $i < $imax; $i++ )
648 :     {
649 :     my $inform = 0;
650 :     my @align2 = map { my $c = substr( $_->[-1], $i, 1 );
651 :     $inform++ if ( $c =~ m/[ACDEFGHIKLMNPQRSTVWY]/i );
652 :     [ $_->[0], $c ]
653 :     }
654 : golsen 1.3 @$align;
655 : golsen 1.1
656 :     # Only analyze the rate if there are 4 or more informative sequences:
657 :    
658 :     if ( $inform >= 4 )
659 :     {
660 :     my @results = proml( \@align2, \%proml_opts );
661 :    
662 :     my ( $best ) = sort { $b->[1] <=> $a->[1] }
663 :     map { [ $_, @{ shift @results }[1] ] } # get the likelihoods
664 :     @cat_vals;
665 :    
666 :     # printf STDERR "%6d %2d => %12.4f\n", $i+1, @$best; ## DEBUG ##
667 :     push @categories, $best->[0];
668 :     push @weights, 1;
669 :     }
670 :     else
671 :     {
672 :     push @categories, 9;
673 :     push @weights, 0;
674 :     }
675 :     }
676 :    
677 :     # Find the minimum category value to appear:
678 :    
679 :     my ( $mincat ) = sort { $a <=> $b } @categories;
680 :     my $adjust = $mincat - 1;
681 :    
682 :     @categories = map { min( $_ - $adjust, 9 ) } @categories;
683 :     @rates = @rates[ $adjust .. ( $adjust+8 ) ];
684 :    
685 :     # Return category and weight data:
686 :    
687 :     ( [ \@rates, join( '', @categories ) ], join( '', @weights ) )
688 :     }
689 :    
690 :    
691 :     #===============================================================================
692 : golsen 1.4 # optimize_protein_coef_of_var -- Find optimal tree and coeficient of variation.
693 :     #
694 :     # ( $coef, $tree, $score ) = optimize_protein_coef_of_var( \@alignment, \%options )
695 :     #
696 : golsen 1.5 # Options are same as proml, plus:
697 :     #
698 :     # max => max_coef # Largest value of search range (D = 2)
699 :     # step => step_size # Step size for intial grid search (D = 0.32)
700 :     # verbose => boolean # Print progress to STDERR (D = 0)
701 :     #
702 : golsen 1.4 #===============================================================================
703 :     sub optimize_protein_coef_of_var
704 :     {
705 :     my ( $align, $opts ) = @_;
706 :     $align && ref $align eq 'ARRAY' && @$align > 3 && $align->[0] && ref $align->[0] eq 'ARRAY'
707 :     or return ();
708 :    
709 : golsen 1.5 my $max_coef = $opts->{ max } || $opts->{ maximum } || 2;
710 :     my $step = $opts->{ step } || $opts->{ step0 } || 0.32;
711 :     my $verbose = $opts->{ verbose } || 0;
712 :    
713 :     my @coefs;
714 :     for ( my $c = 0; $c <= $max_coef; $c += $step ) { push @coefs, $c }
715 : golsen 1.4 my %tree_score;
716 :     while ( 1 )
717 :     {
718 :     # Evaluate any values without a score (and tree)
719 :     foreach my $coef ( grep { ! $tree_score{ $_ } } @coefs )
720 :     {
721 :     $opts->{ coef_of_var } = $coef;
722 :     ( $tree_score{ $coef } ) = gjophylip::proml( $align, $opts );
723 : golsen 1.5 print STDERR "Coef = $coef, LnLikelihood = $tree_score{$coef}->[1]\n" if $verbose;
724 : golsen 1.4 }
725 :    
726 :     # Sort highest score to lowest score
727 :     @coefs = sort { $tree_score{$b}->[1] <=> $tree_score{$a}->[1] } @coefs;
728 :     last if $step <= 0.011; # Done if reached desired resolution
729 :    
730 :     splice @coefs, 3; # Keep best 3
731 :     $step *= 0.5; # Go to half the step size
732 :     @coefs = sort { $a <=> $b } @coefs; # Put in numerical order
733 :     push @coefs, $coefs[0]+$step, $coefs[2]-$step; # Add intermediate coefs
734 :     }
735 :    
736 :     ( $coefs[0], @{ $tree_score{$coefs[0]} } ); # coef_of_var, tree, lnlikelihood
737 :     }
738 :    
739 :    
740 :     #===============================================================================
741 : golsen 1.3 # dnadist -- A perl interface to the dnadist program
742 : golsen 1.1 #
743 : golsen 1.3 # $dists = dnadist( \@alignment, \%options )
744 :     # $dists = dnadist( %options )
745 : golsen 1.1 #
746 :     # @alignment = ( [ id, seq ], [ id, seq ] ... )
747 :     # or
748 :     # ( [ id, def, seq ], [ id, def, seq ] ... )
749 :     #
750 :     # Options:
751 :     #
752 : golsen 1.3 # dnadist:
753 : golsen 1.1 # alignment => \@alignment supply the alignment as an option
754 :     # alpha => float alpha parameter of gamma distribution (0.5 - inf)
755 :     # categories => [ [ rates ], site_categories ]
756 :     # coef_of_var => float 1/sqrt(alpha) for gamma distribution (D = 0)
757 :     # invar_frac => 0 - 1 fraction of site that are invariant
758 :     # model => model evolution model JTT (D) | PMB | PAM
759 : golsen 1.3 # rate_hmm => rate_prob_pairs does not exist in current dnadist
760 : golsen 1.1 # weights => site_weights
761 :     #
762 :     # other:
763 :     # keep_duplicates => bool do not remove duplicate sequences [NOT IMPLIMENTED]
764 : golsen 1.3 # dnadist => dnadist allows fully defined path
765 : golsen 1.1 # tmp => directory directory for tmp_dir
766 :     # tmp_dir => directory directory for temporary files
767 :     #
768 :     # tmp_dir is created and deleted unless its name is supplied, and it already
769 :     # exists.
770 :     #===============================================================================
771 :    
772 : golsen 1.3 sub dnadist
773 : golsen 1.1 {
774 :     my $align;
775 : golsen 1.3 $align = shift @_ if ( ref( $_[0] ) eq 'ARRAY' );
776 : golsen 1.1
777 :     my %options;
778 :     if ( $_[0] )
779 :     {
780 :     my %opts = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
781 :     %options = map { canonical_key( $_ ) => $opts{ $_ } } keys %opts;
782 :     }
783 :    
784 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
785 :     # Clean up a local copy of the alignment. Assign local ids.
786 :     # Alignment can also be an option:
787 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
788 : golsen 1.1
789 :     $align ||= $options{ alignment } || $options{ align };
790 :    
791 : golsen 1.3 my ( $id, $local_id );
792 :     ( $align, $id, $local_id ) = &process_dna_alignment( $align );
793 : golsen 1.1
794 : golsen 1.3 if ( ! $align )
795 :     {
796 :     print STDERR " gjophylip::dnadist requires an alignment.\n";
797 :     return ();
798 :     }
799 : golsen 1.1
800 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
801 :     # Process dnadist options:
802 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
803 : golsen 1.1
804 :     my $categories = $options{ categorie }; # The 's' is lost in canonicalizing keys
805 :     if ( $categories )
806 :     {
807 : golsen 1.3 ( $categories = &process_categories( $categories ) )
808 :     or print STDERR " gjophylip::dnadist bad categories option value.\n"
809 :     and return ();
810 : golsen 1.1 }
811 :    
812 :     my $coef_of_var = $options{ coefofvar }
813 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
814 :     || 0;
815 :     if ( $coef_of_var < 0 )
816 :     {
817 : golsen 1.3 print STDERR "gjophylip::dnadist coef_of_var option value must be >= 0\n";
818 : golsen 1.1 return ();
819 :     }
820 :    
821 : golsen 1.3 my $frequencies = $options{ frequencie };
822 :     if ( $frequencies && ( ref( $frequencies ) ne 'ARRAY' || ( @$frequencies != 4 ) ) )
823 : golsen 1.1 {
824 : golsen 1.3 print STDERR "gjophylip::dnadist frequences must be [ f(A), f(C), f(G), f(T) ]\n";
825 : golsen 1.1 return ();
826 :     }
827 :    
828 : golsen 1.3 my $invar_frac = $options{ invarfrac } || 0;
829 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
830 : golsen 1.1 {
831 : golsen 1.3 print STDERR "gjophylip::dnadist invar_frac option value must be >= 0 and < 1\n";
832 : golsen 1.1 return ();
833 :     }
834 :    
835 : golsen 1.3 my $model = &process_dna_model( \%options );
836 : golsen 1.1
837 :     my $persistance = $options{ persistance } || 0;
838 :     if ( $persistance && ( $persistance <= 1 ) )
839 :     {
840 : golsen 1.3 print STDERR "gjophylip::dnadist persistance option value must be > 1\n";
841 :     return ();
842 :     }
843 :    
844 :     # Does not exist in current dnadist program
845 :     #
846 :     # my $rate_hmm = $options{ ratehmm };
847 :     # if ( $rate_hmm )
848 :     # {
849 :     # ( $rate_hmm = &process_rate_hmm( $rate_hmm ) )
850 :     # or print STDERR " gjophylip::proml bad rate_hmm value\n"
851 :     # and return ();
852 :     # }
853 :    
854 :     my $tt_param = $options{ transitiontransversion } || $options{ ttparam };
855 :     if ( $tt_param && ( $tt_param <= 0 ) )
856 :     {
857 :     print STDERR "gjophylip::dnadist tt_param option value must be > 0\n";
858 : golsen 1.1 return ();
859 :     }
860 :    
861 : golsen 1.2 my $weights = $options{ weight };
862 : golsen 1.1
863 :    
864 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
865 :     # Options that are not dnadist options per se:
866 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
867 : golsen 1.1
868 : golsen 1.3 my $dnadist = SeedAware::executable_for( $options{ dnadist } || 'dnadist' );
869 : golsen 1.1
870 : golsen 1.3 my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( \%options );
871 : golsen 1.1 $tmp_dir or return ();
872 :    
873 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
874 : golsen 1.1 # Write the files and run the program:
875 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
876 : golsen 1.1
877 :     my $cwd = $ENV{ cwd } || `pwd`;
878 :     chomp $cwd;
879 :     chdir $tmp_dir;
880 :    
881 : golsen 1.3 my $distances = run_dnadist( align => $align,
882 :     categories => $categories,
883 :     coef_of_var => $coef_of_var,
884 :     program => $dnadist,
885 :     frequencies => $frequencies,
886 :     invar_frac => $invar_frac,
887 :     model => $model,
888 :     persistance => $persistance,
889 :     # rate_hmm => $rate_hmm,
890 :     tt_param => $tt_param,
891 :     weights => $weights,
892 :     );
893 :    
894 :     # We are done, go back to the original directory:
895 :    
896 :     chdir $cwd;
897 :    
898 :     # Delete the temporary directory unless it already existed:
899 :    
900 :     system( '/bin/rm', '-r', $tmp_dir ) if ! $save_tmp;
901 :    
902 :     return $distances;
903 :     }
904 :    
905 :    
906 :     #-------------------------------------------------------------------------------
907 :     # A local routine to run dnadist
908 :     #-------------------------------------------------------------------------------
909 :     sub run_dnadist
910 :     {
911 :     my %data = @_;
912 :    
913 :     unlink 'outfile' if -f 'outfile'; # Just checking
914 :    
915 :     &write_seq_infile( @{$data{align}} )
916 :     or print STDERR "gjophylip::run_dnadist: Could write infile\n"
917 :     and return undef;
918 :    
919 :     open( DNAD, ">dnadist_cmd" )
920 :     or print STDERR "gjophylip::run_dnadist: Could not open command file for $data{dnadist}\n"
921 :     and return undef;
922 :    
923 :     # Some sanity checks (dnadist does not like to get something that it
924 :     # does not expect):
925 :    
926 :     if ( $data{model} ne 'F84' )
927 :     {
928 :     $data{frequencies} = undef;
929 :     }
930 : golsen 1.1
931 : golsen 1.3 if ( $data{model} =~ m/LogDet/i
932 :     || $data{model} =~ m/identity/i
933 :     )
934 : golsen 1.1 {
935 : golsen 1.3 $data{categories} = undef;
936 :     $data{coef_of_var} = undef;
937 :     $data{invar_frac} = undef;
938 : golsen 1.1 }
939 :    
940 :    
941 : golsen 1.3 # Start writing options for dnadist (order matters):
942 :    
943 :     # Cycle through to correct model (F84 is default):
944 :    
945 :     print DNAD "D\n" if $data{model} =~ m/Kimura/i;
946 :     print DNAD "D\nD\n" if $data{model} =~ m/JC/i;
947 :     print DNAD "D\nD\nD\n" if $data{model} =~ m/LogDet/i;
948 :     print DNAD "D\nD\nD\nD\n" if $data{model} =~ m/identity/i;
949 :    
950 :     if ( $data{categories} )
951 :     {
952 :     &write_categories( $data{categories}->[1] )
953 :     or print STDERR "gjophylip::run_dnadist: Could not write categories\n"
954 :     and return undef;
955 :     print DNAD "C\n",
956 :     scalar @{$data{categories}->[0]}, "\n",
957 :     join( ' ', map { sprintf( "%.6f", $_ ) } @{ $data{categories}->[0] } ), "\n";
958 :     }
959 :     elsif ( $data{invar_frac} || $data{coef_of_var} )
960 :     {
961 :     print DNAD "G\n";
962 :     print DNAD "G\n" if $data{invar_frac};
963 :     # print DNAD "A\n", "$data{persistance}\n" if $data{persistance};
964 :     }
965 : golsen 1.1
966 : golsen 1.3 if ( $data{tt_param} && ( $data{model} =~ m/Kimura/i
967 :     || $data{model} =~ m/Kimura/i
968 :     )
969 :     )
970 :     {
971 :     print DNAD "T\n";
972 :     printf DNAD "%.6f\n", $data{tt_param};
973 :     }
974 : golsen 1.1
975 : golsen 1.3 if ( $data{weights} )
976 : golsen 1.1 {
977 : golsen 1.3 &write_weights( $data{weights} )
978 :     or print STDERR "gjophylip::run_dnadist: Could not write weights\n"
979 :     and return undef;
980 :     print DNAD "W\n";
981 :     }
982 :    
983 :     # We have identified all the options; indicate this with 'Y':
984 :    
985 :     print DNAD "Y\n";
986 : golsen 1.1
987 : golsen 1.3 # Becuase of the options interface, some values still must be supplied:
988 : golsen 1.1
989 : golsen 1.3 if ( $data{invar_frac} || $data{coef_of_var} )
990 :     {
991 :     print DNAD "$data{coef_of_var}\n";
992 :     print DNAD "$data{invar_frac}\n" if $data{invar_frac};
993 : golsen 1.1 }
994 :    
995 : golsen 1.3 close DNAD;
996 :    
997 :     #---------------------------------------------------------------------------
998 :     # Run the program
999 :     #---------------------------------------------------------------------------
1000 :    
1001 :     my $program = $data{ program } || $data{ dnadist } || 'dnadist';
1002 :     my $redirects = { stdin => 'dnadist_cmd',
1003 :     stdout => 'dnadist_out', # Gets deleted with tmp dir
1004 :     stderr => 'dnadist_err' # Gets deleted with tmp dir
1005 :     };
1006 :     SeedAware::system_with_redirect( $program, $redirects )
1007 :     and print STDERR "gjophylip::run_dnadist: Could not run $program.\n"
1008 :     and return undef;
1009 :    
1010 :     #---------------------------------------------------------------------------
1011 :     # Gather the results
1012 :     #---------------------------------------------------------------------------
1013 :    
1014 :     my $distances = read_distances();
1015 :     $distances or print STDERR "gjophylip::run_dnadist: Could not read 'outfile'\n"
1016 :     and return undef;
1017 :    
1018 :     return $distances;
1019 : golsen 1.1 }
1020 :    
1021 :    
1022 :     #===============================================================================
1023 : golsen 1.3 # dnadist_neighbor -- A perl interface to the dnadist and neighbor programs
1024 : golsen 1.1 #
1025 : golsen 1.3 # $tree = dnadist_neighbor( \@alignment, \%options )
1026 :     # $tree = dnadist_neighbor( \@alignment, %options )
1027 :     # $tree = dnadist_neighbor( \%options )
1028 :     # $tree = dnadist_neighbor( %options )
1029 : golsen 1.1 #
1030 :     # @alignment = ( [ id, seq ], [ id, seq ] ... )
1031 :     # or
1032 :     # ( [ id, def, seq ], [ id, def, seq ] ... )
1033 :     #
1034 :     # Options:
1035 :     #
1036 : golsen 1.3 # dnadist:
1037 : golsen 1.1 # alignment => \@alignment supply the alignment as an option
1038 :     # alpha => float alpha parameter of gamma distribution (0.5 - inf)
1039 :     # categories => [ [ rates ], site_categories ]
1040 :     # coef_of_var => float 1/sqrt(alpha) for gamma distribution (D = 0)
1041 :     # invar_frac => 0 - 1 fraction of site that are invariant
1042 :     # model => model evolution model JTT (D) | PMB | PAM
1043 : golsen 1.3 # rate_hmm => rate_prob_pairs does not exist in current dnadist
1044 : golsen 1.1 # weights => site_weights
1045 :     #
1046 : golsen 1.3 # neighbor:
1047 :     # jumble_seed => odd_int jumble random seed
1048 :     #
1049 : golsen 1.1 # other:
1050 : golsen 1.3 # definitions => bool include definition in sequence label
1051 :     # dnadist => dnadist allows fully defined path
1052 :     # keep_duplicates => bool do not remove duplicate sequences [NOT IMPLIMENTED]
1053 :     # neighbor => neighbor allows fully defined path
1054 :     # no_ids => bool leave the sequence id out of the sequence label
1055 : golsen 1.1 # tmp => directory directory for tmp_dir
1056 :     # tmp_dir => directory directory for temporary files
1057 : golsen 1.3 # tree_format => overbeek | gjo | fig format of output tree
1058 : golsen 1.1 #
1059 :     # tmp_dir is created and deleted unless its name is supplied, and it already
1060 :     # exists.
1061 :     #===============================================================================
1062 : golsen 1.3
1063 :     sub dnadist_neighbor
1064 : golsen 1.1 {
1065 :     my $align;
1066 : golsen 1.3 $align = shift @_ if ( ref( $_[0] ) eq 'ARRAY' );
1067 : golsen 1.1
1068 :     my %options;
1069 :     if ( $_[0] )
1070 :     {
1071 :     my %opts = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
1072 :     %options = map { canonical_key( $_ ) => $opts{ $_ } } keys %opts;
1073 :     }
1074 :    
1075 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1076 :     # Clean up a local copy of the alignment. Assign local ids.
1077 :     # Alignment can also be an option:
1078 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1079 : golsen 1.1
1080 :     $align ||= $options{ alignment } || $options{ align };
1081 :    
1082 : golsen 1.3 my ( $align2, $id, $local_id ) = &process_dna_alignment( $align );
1083 : golsen 1.1
1084 : golsen 1.3 if ( ! $align2 )
1085 :     {
1086 :     print STDERR " gjophylip::dnadist_neighbor requires an alignment.\n";
1087 :     return ();
1088 :     }
1089 : golsen 1.1
1090 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1091 :     # Process dnadist_neighbor options:
1092 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1093 : golsen 1.1
1094 :     my $categories = $options{ categorie }; # The 's' is lost in canonicalizing keys
1095 :     if ( $categories )
1096 :     {
1097 : golsen 1.3 ( $categories = &process_categories( $categories ) )
1098 :     or print STDERR " gjophylip::dnadist_neighbor bad categories option value.\n"
1099 :     and return ();
1100 :     }
1101 :    
1102 :     my $coef_of_var = $options{ coefofvar }
1103 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
1104 :     || 0;
1105 :     if ( $coef_of_var < 0 )
1106 :     {
1107 :     print STDERR "gjophylip::dnadist_neighbor coef_of_var option value must be >= 0\n";
1108 :     return ();
1109 :     }
1110 :    
1111 :     my $no_ids = $options{ noid } || 0;
1112 :     my $definitions = $options{ definition } || $no_ids; # No id requires definition
1113 :     if ( $definitions )
1114 :     {
1115 :     my %new_id = map { $_->[0] => ( ! length $_->[1] ? $_->[0] : $no_ids ? $_->[1] : "$_->[0] $_->[1]" ) }
1116 :     map { [ $_->[0], @$_ == 3 && defined $_->[1] ? $_->[1] : '' ] }
1117 :     @$align;
1118 :     foreach ( keys %$id ) { $id->{ $_ } = $new_id{ $id->{ $_ } } }
1119 :     }
1120 :    
1121 :     my $frequencies = $options{ frequencie };
1122 :     if ( $frequencies && ( ref( $frequencies ) ne 'ARRAY' || ( @$frequencies != 4 ) ) )
1123 :     {
1124 :     print STDERR "gjophylip::dnadist_neighbor frequences must be [ f(A), f(C), f(G), f(T) ]\n";
1125 :     return ();
1126 :     }
1127 :    
1128 :     my $invar_frac = $options{ invarfrac } || 0;
1129 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
1130 :     {
1131 :     print STDERR "gjophylip::dnadist_neighbor invar_frac option value must be >= 0 and < 1\n";
1132 :     return ();
1133 :     }
1134 :    
1135 :     my $jumble_seed = int( $options{ jumbleseed } ) || 0;
1136 :     if ( $jumble_seed && ( ( $jumble_seed < 0) || ( $jumble_seed % 2 != 1 ) ) )
1137 :     {
1138 :     print STDERR "gjophylip::dnadist_neighbor jumble_seed option value must be an odd number > 0\n";
1139 :     return ();
1140 :     }
1141 :    
1142 :     my $model = &process_dna_model( \%options );
1143 :    
1144 :     my $persistance = $options{ persistance } || 0;
1145 :     if ( $persistance && ( $persistance <= 1 ) )
1146 :     {
1147 :     print STDERR "gjophylip::dnadist_neighbor persistance option value must be > 1\n";
1148 :     return ();
1149 :     }
1150 :    
1151 :     # Does not exist in current dnadist program
1152 :     #
1153 :     # my $rate_hmm = $options{ ratehmm };
1154 :     # if ( $rate_hmm )
1155 :     # {
1156 :     # ( $rate_hmm = &process_rate_hmm( $rate_hmm ) )
1157 :     # or print STDERR " gjophylip::proml bad rate_hmm value\n"
1158 :     # and return ();
1159 :     # }
1160 :    
1161 :     my $tt_param = $options{ transitiontransversion } || $options{ ttparam };
1162 :     if ( $tt_param && ( $tt_param <= 0 ) )
1163 :     {
1164 :     print STDERR "gjophylip::dnadist_neighbor tt_param option value must be > 0\n";
1165 :     return ();
1166 :     }
1167 :    
1168 :     my $weights = $options{ weight };
1169 :    
1170 :    
1171 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1172 :     # Options that are not dnadist_neighbor options per se:
1173 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1174 :    
1175 :     my $dnadist = SeedAware::executable_for( $options{ dnadist } || 'dnadist' );
1176 :    
1177 :     my $neighbor = $options{ neighbor } || 'neighbor';
1178 :    
1179 :     my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( \%options );
1180 :     $tmp_dir or return ();
1181 :    
1182 :     my $tree_format = process_tree_format( \%options );
1183 :    
1184 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1185 :     # Write the files and run the program:
1186 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1187 :    
1188 :     my $cwd = $ENV{ cwd } || `pwd`;
1189 :     chomp $cwd;
1190 :     chdir $tmp_dir;
1191 :    
1192 :     my $distances = run_dnadist( align => $align2,
1193 :     categories => $categories,
1194 :     coef_of_var => $coef_of_var,
1195 :     program => $dnadist,
1196 :     frequencies => $frequencies,
1197 :     invar_frac => $invar_frac,
1198 :     model => $model,
1199 :     persistance => $persistance,
1200 :     # rate_hmm => $rate_hmm,
1201 :     tt_param => $tt_param,
1202 :     weights => $weights,
1203 :     );
1204 :    
1205 :     my $tree = undef;
1206 :     if ( $distances )
1207 :     {
1208 :     $tree = run_neighbor( distances => $distances,
1209 :     jumble_seed => $jumble_seed,
1210 :     program => $neighbor
1211 :     );
1212 :     }
1213 :    
1214 :     # We are done, go back to the original directory:
1215 :    
1216 :     chdir $cwd;
1217 :    
1218 :     # Delete the temporary directory unless it already existed:
1219 :    
1220 :     system( '/bin/rm', '-r', $tmp_dir ) if ! $save_tmp;
1221 :    
1222 :     if ( $tree )
1223 :     {
1224 :     # Returned trees have our labels:
1225 :    
1226 :     gjonewicklib::newick_relabel_nodes( $tree, $id );
1227 :    
1228 :     if ( $tree_format =~ m/overbeek/i )
1229 : golsen 1.1 {
1230 : golsen 1.3 $tree = gjonewicklib::gjonewick_to_overbeek( $tree );
1231 : golsen 1.1 }
1232 : golsen 1.3 }
1233 :    
1234 :     return $tree;
1235 :     }
1236 : golsen 1.1
1237 :    
1238 : golsen 1.3 #===============================================================================
1239 :     # protdist -- A perl interface for the PHYLIP protdist program:
1240 :     #
1241 :     # \@distance_matrix = protdist( \@alignment, \%options )
1242 :     # \@distance_matrix = protdist( \@alignment, %options )
1243 :     # \@distance_matrix = protdist( \%options )
1244 :     # \@distance_matrix = protdist( %options )
1245 :     #
1246 :     # @alignment = ( [ id, seq ], [ id, seq ] ... )
1247 :     # or
1248 :     # ( [ id, def, seq ], [ id, def, seq ] ... )
1249 :     #
1250 :     # @distance_matrix = ( [ id1, dist11, dist12, dist13, ... ],
1251 :     # [ id2, dist21, dist22, dist14, ... ],
1252 :     # ...
1253 :     # )
1254 :     #
1255 :     # Options:
1256 :     #
1257 :     # protdist
1258 :     # alignment => \@alignment supply the alignment as an option
1259 :     # alpha => float alpha parameter of gamma distribution (0.5 - inf)
1260 :     # categories => [ [ rates ], site_categories ]
1261 :     # coef_of_var => float 1/sqrt(alpha) for gamma distribution (D = 0)
1262 :     # invar_frac => 0 - 1 fraction of site that are invariant
1263 :     # model => model evolution model JTT (D) | PMB | PAM
1264 :     # persistance => float persistance length of rate category
1265 :     # rate_hmm => [ [ rates ], [ probabilies ] ]
1266 :     # weights => site_weights
1267 :     #
1268 :     # other:
1269 :     # protdist => protdist allows fully defined path
1270 :     # tmp => directory directory for tmp_dir
1271 :     # tmp_dir => directory directory for temporary files
1272 :     #
1273 :     # tmp_dir is created and deleted unless its name is supplied, and it already
1274 :     # exists.
1275 :     #===============================================================================
1276 :     sub protdist
1277 :     {
1278 :     my $align;
1279 :     $align = shift @_ if ( ref( $_[0] ) eq 'ARRAY' );
1280 :    
1281 :     my %options;
1282 :     if ( $_[0] )
1283 :     {
1284 :     my %opts = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
1285 :     %options = map { canonical_key( $_ ) => $opts{ $_ } } keys %opts;
1286 :     }
1287 :    
1288 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1289 :     # Clean up a local copy of the alignment. Assign local ids.
1290 :     # Alignment can also be an option:
1291 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1292 :    
1293 :     $align ||= $options{ alignment } || $options{ align };
1294 :    
1295 :     my ( $id, $local_id );
1296 :     ( $align, $id, $local_id ) = &process_protein_alignment( $align );
1297 :    
1298 :     if ( ! $align )
1299 :     {
1300 :     print STDERR " gjophylip::protdist requires an alignment.\n";
1301 :     return ();
1302 :     }
1303 :    
1304 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1305 :     # Process protdist options:
1306 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1307 :    
1308 :     my $categories = $options{ categorie }; # The 's' is lost in canonicalizing keys
1309 :     if ( $categories )
1310 :     {
1311 :     ( $categories = &process_categories( $categories ) )
1312 :     or print STDERR " gjophylip::protdist bad categories option value.\n"
1313 :     and return ();
1314 : golsen 1.1 }
1315 :    
1316 :     my $coef_of_var = $options{ coefofvar }
1317 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
1318 :     || 0;
1319 :     if ( $coef_of_var < 0 )
1320 :     {
1321 :     print STDERR "gjophylip::protdist coef_of_var option value must be >= 0\n";
1322 :     return undef;
1323 :     }
1324 :    
1325 :     my $invar_frac = $options{ invarfrac } || 0;
1326 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
1327 :     {
1328 :     print STDERR "gjophylip::protdist invar_frac option value must be >= 0 and < 1\n";
1329 :     return undef;
1330 :     }
1331 :    
1332 : golsen 1.3 my $model = &process_protein_model( \%options );
1333 : golsen 1.1
1334 : golsen 1.3 my $persistance = $options{ persistance } || 0;
1335 : golsen 1.1 if ( $persistance && ( $persistance <= 1 ) )
1336 :     {
1337 :     print STDERR "gjophylip::protdist persistance option value must be > 1\n";
1338 :     return undef;
1339 :     }
1340 :    
1341 : golsen 1.3 # Does not exist in current protdist
1342 :     #
1343 :     # my $rate_hmm = $options{ ratehmm };
1344 :     # if ( $rate_hmm )
1345 :     # {
1346 :     # ( $rate_hmm = &process_rate_hmm( $rate_hmm ) )
1347 :     # or print STDERR " gjophylip::proml bad rate_hmm value\n"
1348 :     # and return ();
1349 :     # }
1350 :    
1351 :     my $weights = $options{ weight };
1352 : golsen 1.1
1353 :    
1354 :     #---------------------------------------------------------------------------
1355 :     # Options that are not protdist options per se:
1356 :     #---------------------------------------------------------------------------
1357 :    
1358 : golsen 1.3 my $protdist = SeedAware::executable_for( $options{ protdist } || 'protdist' );
1359 : golsen 1.1
1360 : golsen 1.3 my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( \%options );
1361 : golsen 1.1 $tmp_dir or return ();
1362 :    
1363 :     #---------------------------------------------------------------------------
1364 :     # Write the files and run the program:
1365 :     #---------------------------------------------------------------------------
1366 :    
1367 :     my $cwd = $ENV{ cwd } || `pwd`;
1368 :     chomp $cwd;
1369 :     chdir $tmp_dir;
1370 :    
1371 : golsen 1.3 my $distances = run_protdist( align => $align,
1372 :     categories => $categories,
1373 :     coef_of_var => $coef_of_var,
1374 :     invar_frac => $invar_frac,
1375 :     model => $model,
1376 :     persistance => $persistance,
1377 :     program => $protdist,
1378 :     # rate_hmm => $rate_hmm,
1379 :     weights => $weights,
1380 : golsen 1.1 );
1381 :    
1382 :     chdir $cwd;
1383 :    
1384 :     # Delete the temporary directory unless it already existed:
1385 :    
1386 : golsen 1.3 system( '/bin/rm', '-r', $tmp_dir ) if ! $save_tmp;
1387 : golsen 1.1
1388 :     # Fix the taxon labels:
1389 :    
1390 : golsen 1.3 foreach ( @$distances ) { $_->[0] = $id->{ $_->[0] } }
1391 : golsen 1.1
1392 :     return $distances;
1393 :     }
1394 :    
1395 :    
1396 :     #-------------------------------------------------------------------------------
1397 :     # A local routine to run protdist
1398 :     #-------------------------------------------------------------------------------
1399 : golsen 1.3 sub run_protdist
1400 :     {
1401 :     my %data = @_;
1402 :    
1403 :     unlink 'outfile' if -f 'outfile'; # Just checking
1404 :    
1405 :     &write_seq_infile( @{$data{align}} )
1406 :     or print STDERR "gjophylip::run_protdist: Could write infile\n"
1407 :     and return undef;
1408 :    
1409 :     open( PROTD, ">protdist_cmd" )
1410 :     or print STDERR "gjophylip::run_protdist: Could not open command file for $data{protdist}\n"
1411 :     and return undef;
1412 :    
1413 :     # Start writing options for protdist:
1414 :    
1415 :     if ( $data{categories} )
1416 :     {
1417 :     &write_categories( $data{categories}->[1] )
1418 :     or print STDERR "gjophylip::run_protdist: Could not write categories\n"
1419 :     and return undef;
1420 :     print PROTD "C\n",
1421 :     scalar @{$data{categories}->[0]}, "\n",
1422 :     join( ' ', map { sprintf( "%.6f", $_ ) } @{ $data{categories}->[0] } ), "\n";
1423 :     }
1424 :    
1425 :     if ( $data{invar_frac} || $data{coef_of_var} )
1426 :     {
1427 :     print PROTD "G\n";
1428 :     print PROTD "G\n" if $data{invar_frac};
1429 :     print PROTD "A\n", "$data{persistance}\n" if $data{persistance};
1430 :     }
1431 :    
1432 :     print PROTD "P\n" if $data{model} =~ m/PMB/i;
1433 :     print PROTD "P\nP\n" if $data{model} =~ m/PAM/i;
1434 :     print PROTD "P\nP\nP\n" if $data{model} =~ m/Kimura/i;
1435 :    
1436 :     if ( $data{weights} )
1437 :     {
1438 :     &write_weights( $data{weights} )
1439 :     or print STDERR "gjophylip::run_protdist: Could not write weights\n"
1440 :     and return undef;
1441 :     print PROTD "W\n";
1442 :     }
1443 :    
1444 :     # We have identified all the options; indicate this with 'Y':
1445 :    
1446 :     print PROTD "Y\n";
1447 :    
1448 :     # Becuase of the options interface, some values still must be supplied:
1449 :    
1450 :     if ( $data{invar_frac} || $data{coef_of_var} )
1451 :     {
1452 :     print PROTD "$data{coef_of_var}\n";
1453 :     print PROTD "$data{invar_frac}\n" if $data{invar_frac};
1454 :     }
1455 :    
1456 :     close PROTD;
1457 :    
1458 :     #---------------------------------------------------------------------------
1459 :     # Run the program
1460 :     #---------------------------------------------------------------------------
1461 :    
1462 :     my $program = $data{ program } || $data{ protdist } || 'protdist';
1463 :     my $redirects = { stdin => 'protdist_cmd',
1464 :     stdout => 'protdist_out', # Gets deleted with tmp dir
1465 :     stderr => 'protdist_err' # Gets deleted with tmp dir
1466 :     };
1467 :     SeedAware::system_with_redirect( $program, $redirects )
1468 :     and print STDERR "gjophylip::run_protdist: Could not run $program.\n"
1469 :     and return undef;
1470 :    
1471 :     #---------------------------------------------------------------------------
1472 :     # Gather the results
1473 :     #---------------------------------------------------------------------------
1474 :    
1475 :     my $distances = read_distances();
1476 :     $distances or print STDERR "gjophylip::run_protdist: Could not read 'outfile'\n"
1477 :     and return undef;
1478 :    
1479 :     return $distances;
1480 :     }
1481 :    
1482 :    
1483 :     #===============================================================================
1484 :     # protdist_neighbor -- A perl interface to the protdist and neighbor programs
1485 :     #
1486 :     # $tree = protdist_neighbor( \@alignment, \%options )
1487 :     # $tree = protdist_neighbor( \@alignment, %options )
1488 :     # $tree = protdist_neighbor( \%options )
1489 :     # $tree = protdist_neighbor( %options )
1490 :     #
1491 :     # @alignment = ( [ id, seq ], [ id, seq ] ... )
1492 :     # or
1493 :     # ( [ id, def, seq ], [ id, def, seq ] ... )
1494 :     #
1495 :     # Options:
1496 :     #
1497 :     # protdist:
1498 :     # alignment => \@alignment supply the alignment as an option
1499 :     # alpha => float alpha parameter of gamma distribution (0.5 - inf)
1500 :     # categories => [ [ rates ], site_categories ]
1501 :     # coef_of_var => float 1/sqrt(alpha) for gamma distribution (D = 0)
1502 :     # invar_frac => 0 - 1 fraction of site that are invariant
1503 :     # model => model evolution model JTT (D) | PMB | PAM
1504 :     # rate_hmm => rate_prob_pairs does not exist in current protdist
1505 :     # weights => site_weights
1506 :     #
1507 :     # neighbor:
1508 :     # jumble_seed => odd_int jumble random seed
1509 :     #
1510 :     # other:
1511 :     # definitions => bool include definition in sequence label
1512 :     # keep_duplicates => bool do not remove duplicate sequences [NOT IMPLIMENTED]
1513 :     # no_ids => bool leave the sequence id out of the sequence label
1514 :     # protdist => protdist allows fully defined path
1515 :     # neighbor => neighbor allows fully defined path
1516 :     # tmp => directory directory for tmp_dir
1517 :     # tmp_dir => directory directory for temporary files
1518 :     # tree_format => overbeek | gjo | fig format of output tree
1519 :     #
1520 :     # tmp_dir is created and deleted unless its name is supplied, and it already
1521 :     # exists.
1522 :     #===============================================================================
1523 :    
1524 :     sub protdist_neighbor
1525 :     {
1526 :     my $align;
1527 :     $align = shift @_ if ( ref( $_[0] ) eq 'ARRAY' );
1528 :    
1529 :     my %options;
1530 :     if ( $_[0] )
1531 :     {
1532 :     my %opts = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
1533 :     %options = map { canonical_key( $_ ) => $opts{ $_ } } keys %opts;
1534 :     }
1535 :    
1536 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1537 :     # Clean up a local copy of the alignment. Assign local ids.
1538 :     # Alignment can also be an option:
1539 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1540 :    
1541 :     $align ||= $options{ alignment } || $options{ align };
1542 :    
1543 :     my ( $align2, $id, $local_id ) = &process_protein_alignment( $align );
1544 :    
1545 :     if ( ! $align2 )
1546 :     {
1547 :     print STDERR " gjophylip::protdist_neighbor requires an alignment.\n";
1548 :     return ();
1549 :     }
1550 :    
1551 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1552 :     # Process protdist_neighbor options:
1553 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1554 :    
1555 :     my $categories = $options{ categorie }; # The 's' is lost in canonicalizing keys
1556 :     if ( $categories )
1557 :     {
1558 :     ( $categories = &process_categories( $categories ) )
1559 :     or print STDERR " gjophylip::protdist_neighbor bad categories option value.\n"
1560 :     and return ();
1561 :     }
1562 :    
1563 :     my $coef_of_var = $options{ coefofvar }
1564 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
1565 :     || 0;
1566 :     if ( $coef_of_var < 0 )
1567 :     {
1568 :     print STDERR "gjophylip::protdist_neighbor coef_of_var option value must be >= 0\n";
1569 :     return ();
1570 :     }
1571 :    
1572 :     my $no_ids = $options{ noid } || 0;
1573 :     my $definitions = $options{ definition } || $no_ids; # No id requires definition
1574 :     if ( $definitions )
1575 :     {
1576 :     my %new_id = map { $_->[0] => ( ! length $_->[1] ? $_->[0] : $no_ids ? $_->[1] : "$_->[0] $_->[1]" ) }
1577 :     map { [ $_->[0], @$_ == 3 && defined $_->[1] ? $_->[1] : '' ] }
1578 :     @$align;
1579 :     foreach ( keys %$id ) { $id->{ $_ } = $new_id{ $id->{ $_ } } }
1580 :     }
1581 :    
1582 :     my $invar_frac = $options{ invarfrac } || 0;
1583 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
1584 :     {
1585 :     print STDERR "gjophylip::protdist_neighbor invar_frac option value must be >= 0 and < 1\n";
1586 :     return ();
1587 :     }
1588 :    
1589 :     my $jumble_seed = int( $options{ jumbleseed } || 0 );
1590 :     if ( $jumble_seed && ( ( $jumble_seed < 0) || ( $jumble_seed % 2 != 1 ) ) )
1591 :     {
1592 :     print STDERR "gjophylip::protdist_neighbor jumble_seed option value must be an odd number > 0\n";
1593 :     return ();
1594 :     }
1595 :    
1596 :     my $model = &process_protein_model( \%options );
1597 :    
1598 :     my $persistance = $options{ persistance } || 0;
1599 :     if ( $persistance && ( $persistance <= 1 ) )
1600 :     {
1601 :     print STDERR "gjophylip::protdist_neighbor persistance option value must be > 1\n";
1602 :     return ();
1603 :     }
1604 :    
1605 :     # Does not exist in current protdist program
1606 :     #
1607 :     # my $rate_hmm = $options{ ratehmm };
1608 :     # if ( $rate_hmm )
1609 :     # {
1610 :     # ( $rate_hmm = &process_rate_hmm( $rate_hmm ) )
1611 :     # or print STDERR " gjophylip::proml bad rate_hmm value\n"
1612 :     # and return ();
1613 :     # }
1614 :    
1615 :     my $weights = $options{ weight };
1616 :    
1617 :    
1618 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1619 :     # Options that are not protdist_neighbor options per se:
1620 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1621 :    
1622 :     my $protdist = SeedAware::executable_for( $options{ protdist } || 'protdist' );
1623 :    
1624 :     my $neighbor = SeedAware::executable_for( $options{ neighbor } || 'neighbor' );
1625 :    
1626 :     my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( \%options );
1627 :     $tmp_dir or return ();
1628 :    
1629 :     my $tree_format = process_tree_format( \%options );
1630 :    
1631 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1632 :     # Write the files and run the program:
1633 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1634 :    
1635 :     my $cwd = $ENV{ cwd } || `pwd`;
1636 :     chomp $cwd;
1637 :     chdir $tmp_dir;
1638 :    
1639 :     my $distances = run_protdist( align => $align2,
1640 :     categories => $categories,
1641 :     coef_of_var => $coef_of_var,
1642 :     invar_frac => $invar_frac,
1643 :     model => $model,
1644 :     persistance => $persistance,
1645 :     program => $protdist,
1646 :     # rate_hmm => $rate_hmm,
1647 :     weights => $weights,
1648 :     );
1649 :    
1650 :     my $tree = undef;
1651 :     if ( $distances )
1652 :     {
1653 :     $tree = run_neighbor( distances => $distances,
1654 :     jumble_seed => $jumble_seed,
1655 :     program => $neighbor
1656 :     );
1657 :     }
1658 :    
1659 :     # We are done, go back to the original directory:
1660 :    
1661 :     chdir $cwd;
1662 :    
1663 :     # Delete the temporary directory unless it already existed:
1664 :    
1665 :     system( '/bin/rm', '-r', $tmp_dir ) if ! $save_tmp;
1666 :    
1667 :     if ( $tree )
1668 :     {
1669 :     # Returned trees have our labels:
1670 :    
1671 :     gjonewicklib::newick_relabel_nodes( $tree, $id );
1672 :    
1673 :     if ( $tree_format =~ m/overbeek/i )
1674 :     {
1675 :     $tree = gjonewicklib::gjonewick_to_overbeek( $tree );
1676 :     }
1677 :     }
1678 :    
1679 :     return $tree;
1680 :     }
1681 :    
1682 :    
1683 :     #===============================================================================
1684 :     # neighbor -- A perl interface to the neighbor program
1685 :     #
1686 :     # $tree = neighbor( \@distances, \%options )
1687 :     # $tree = neighbor( \@distances, %options )
1688 :     # $tree = neighbor( \%options )
1689 :     # $tree = neighbor( %options )
1690 :     #
1691 :     # @distances = ( [ id1, dist11, dist12, ... ],
1692 :     # [ id2, dist21, dist22, ... ],
1693 :     # ...
1694 :     # )
1695 :     #
1696 :     # Options:
1697 :     #
1698 :     # neighbor:
1699 :     # jumble_seed => odd_int jumble random seed
1700 :     # upgma => bool use UPGMA
1701 :     #
1702 :     # other:
1703 :     # neighbor => neighbor allows fully defined path
1704 :     # program => neighbor allows fully defined path
1705 :     # tmp => directory directory for tmp_dir
1706 :     # tmp_dir => directory directory for temporary files
1707 :     # tree_format => overbeek | gjo | fig format of output tree
1708 :     #
1709 :     # tmp_dir is created and deleted unless its name is supplied, and it already
1710 :     # exists.
1711 :     #===============================================================================
1712 :     sub neighbor
1713 :     {
1714 :     my $distances;
1715 :     $distances = shift @_ if ( ref( $_[0] ) eq 'ARRAY' );
1716 :    
1717 :     my %options;
1718 :     if ( $_[0] )
1719 :     {
1720 :     my %opts = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
1721 :     %options = map { canonical_key( $_ ) => $opts{ $_ } } keys %opts;
1722 :     }
1723 : golsen 1.1
1724 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1725 :     # Clean up a local copy of the alignment. Assign local ids.
1726 :     # Alignment can also be an option:
1727 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1728 : golsen 1.1
1729 : golsen 1.3 $distances ||= $options{ distance } || $options{ dist };
1730 : golsen 1.1
1731 : golsen 1.3 my ( $id, $local_id );
1732 :     ( $distances, $id, $local_id ) = &process_distance_matrix( $distances );
1733 : golsen 1.1
1734 : golsen 1.3 if ( ! $distances )
1735 : golsen 1.1 {
1736 : golsen 1.3 print STDERR " gjophylip::neighbor requires a distance matrix.\n";
1737 :     return ();
1738 : golsen 1.1 }
1739 :    
1740 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1741 :     # Process neighbor options:
1742 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1743 :    
1744 :     my $jumble_seed = int( $options{ jumbleseed } || 0 );
1745 :     if ( $jumble_seed && ( ( $jumble_seed < 0) || ( $jumble_seed % 2 != 1 ) ) )
1746 : golsen 1.1 {
1747 : golsen 1.3 print STDERR "gjophylip::neighbor jumble_seed option value must be an odd number > 0\n";
1748 :     return ();
1749 : golsen 1.1 }
1750 :    
1751 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1752 :     # Options that are not neighbor options per se:
1753 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1754 :    
1755 :     my $program = SeedAware::executable_for( $options{ neighbor } || $options{ program } || 'neighbor' );
1756 :    
1757 :     my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( \%options );
1758 :     $tmp_dir or return ();
1759 :    
1760 :     my $tree_format = process_tree_format( \%options );
1761 : golsen 1.1
1762 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1763 :     # Write the files and run the program:
1764 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1765 : golsen 1.1
1766 : golsen 1.3 my $cwd = $ENV{ cwd } || `pwd`;
1767 :     chomp $cwd;
1768 :     chdir $tmp_dir;
1769 : golsen 1.1
1770 : golsen 1.3 my $tree = run_neighbor( distances => $distances,
1771 :     jumble_seed => $jumble_seed,
1772 :     program => $program,
1773 :     upgma => $options{ upgma },
1774 :     );
1775 : golsen 1.1
1776 : golsen 1.3 # We are done, go back to the original directory:
1777 : golsen 1.1
1778 : golsen 1.3 chdir $cwd;
1779 : golsen 1.1
1780 : golsen 1.3 # Delete the temporary directory unless it already existed:
1781 : golsen 1.1
1782 : golsen 1.3 system( '/bin/rm', '-r', $tmp_dir ) if ! $save_tmp;
1783 : golsen 1.1
1784 : golsen 1.3 if ( $tree )
1785 :     {
1786 :     # Returned trees have our labels:
1787 : golsen 1.1
1788 : golsen 1.3 gjonewicklib::newick_relabel_nodes( $tree, $id );
1789 : golsen 1.1
1790 : golsen 1.3 if ( $tree_format =~ m/overbeek/i )
1791 :     {
1792 :     $tree = gjonewicklib::gjonewick_to_overbeek( $tree );
1793 :     }
1794 :     }
1795 : golsen 1.1
1796 : golsen 1.3 return $tree;
1797 : golsen 1.1 }
1798 :    
1799 :    
1800 :     #-------------------------------------------------------------------------------
1801 :     # A local routine to run neighbor
1802 :     #-------------------------------------------------------------------------------
1803 :     sub run_neighbor
1804 :     {
1805 : golsen 1.3 my %data = @_;
1806 : golsen 1.1
1807 :     unlink 'outfile' if -f 'outfile'; # Just checking
1808 :     unlink 'outtree' if -f 'outtree'; # ditto
1809 :    
1810 : golsen 1.3 write_dist_infile( $data{distances} )
1811 : golsen 1.1 or print STDERR "gjophylip::run_neighbor: Could not write distance infile\n"
1812 :     and return undef;
1813 :    
1814 : golsen 1.3 open( NEIGH, ">neighbor_cmd" )
1815 : golsen 1.1 or print STDERR "gjophylip::run_neighbor: Could not open neighbor command file\n"
1816 :     and return undef;
1817 : golsen 1.3 print NEIGH "J\n", "$data{jumble_seed}\n" if $data{jumble_seed};
1818 :     print NEIGH "N\n" if $data{upgma};
1819 : golsen 1.1 print NEIGH "Y\n";
1820 :     close NEIGH;
1821 :    
1822 : golsen 1.3 #---------------------------------------------------------------------------
1823 :     # Run the program
1824 :     #---------------------------------------------------------------------------
1825 :    
1826 :     my $program = $data{ program } || $data{ neighbor } || 'neighbor';
1827 :     my $redirects = { stdin => 'neighbor_cmd',
1828 :     stdout => 'neighbpr_out', # Gets deleted with tmp dir
1829 :     stderr => 'neighbpr_err' # Gets deleted with tmp dir
1830 :     };
1831 :     SeedAware::system_with_redirect( $program, $redirects )
1832 :     and print STDERR "gjophylip::run_neighbor: Could not run $program.\n"
1833 :     and return undef;
1834 :    
1835 :     #---------------------------------------------------------------------------
1836 :     # Gather the results
1837 :     #---------------------------------------------------------------------------
1838 : golsen 1.1
1839 :     my ( $tree ) = gjonewicklib::read_newick_tree( 'outtree' );
1840 :    
1841 :     $tree or print STDERR "gjophylip::run_neighbor: Could not read neighbor outtree file\n";
1842 :    
1843 :     return $tree || undef;
1844 :     }
1845 :    
1846 :    
1847 :     #-------------------------------------------------------------------------------
1848 : golsen 1.2 # Consensus tree:
1849 :     #
1850 :     # $consensus_with_boot_as_branch_len = consensus_tree( @trees, \%options );
1851 :     # $consensus_with_boot_as_branch_len = consensus_tree( \@trees, \%options );
1852 :     # $consensus_with_boot_as_branch_len = consensus_tree( @trees );
1853 :     # $consensus_with_boot_as_branch_len = consensus_tree( \@trees );
1854 :     #
1855 :     # Option keys:
1856 :     #
1857 :     # consense => program # Define alternative name or path to executable
1858 :     # savetmp => boolean # Do not delete temporary directory or files.
1859 :     # tmpdir => directory # Directory for temporary files. If directory
1860 :     # # exists, temporary files are not deleted.
1861 :     # tmp => directory # Directory for tmpdir
1862 :     #
1863 :     #-------------------------------------------------------------------------------
1864 :     sub consensus_tree
1865 :     {
1866 :     my $options = ( ( @_ > 1 ) && ( ref( $_[-1] ) eq 'HASH' ) ) ? pop @_ : {};
1867 :    
1868 :     return undef if ( ! @_ ) || ( ref( $_[0] ) ne 'ARRAY' );
1869 :    
1870 :     my @trees = ( @_ > 1 ) ? @_ : @{ $_[0] };
1871 :     my $n_tree = @trees;
1872 :     my $tip = 'tip000001';
1873 :     my %tips = map { $_ => $tip++ } gjonewicklib::newick_tip_list( $trees[0] );
1874 :    
1875 :     # Make a copy of trees with clean labels
1876 :    
1877 :     my @trees2 = map { gjonewicklib::newick_relabel_tips( gjonewicklib::copy_newick_tree( $_ ), \%tips ) }
1878 :     @trees;
1879 :    
1880 : golsen 1.3 my $program = SeedAware::executable_for( $options->{ consense } || $options->{ program } || 'consense' );
1881 : golsen 1.2
1882 : golsen 1.3 my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( $options );
1883 : golsen 1.2 $tmp_dir or return ();
1884 :    
1885 :     #---------------------------------------------------------------------------
1886 :     # Write the files and run the program:
1887 :     #---------------------------------------------------------------------------
1888 :    
1889 :     my $cwd = $ENV{ cwd } || `pwd`;
1890 :     chomp $cwd;
1891 :     chdir $tmp_dir;
1892 :    
1893 :     my ( $tree )= run_consense( { trees => \@trees2, consense => $program } );
1894 :    
1895 :     chdir $cwd;
1896 :    
1897 :     # Delete the temporary directory unless it already existed:
1898 :    
1899 :     system "/bin/rm -r '$tmp_dir'" if ! $save_tmp;
1900 :    
1901 :     return undef if ! $tree;
1902 :    
1903 :     # run_consense() fixes the root and branch lengths. All we need to do is
1904 :     # fix the taxon labels and return the tree:
1905 :    
1906 :     gjonewicklib::newick_relabel_tips( $tree, { reverse %tips } );
1907 :     }
1908 :    
1909 :    
1910 :     #-------------------------------------------------------------------------------
1911 :     # Consensus tree:
1912 :     #
1913 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, @trees, \%options );
1914 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, \@trees, \%options );
1915 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, @trees );
1916 :     # $tree_with_labeled_nodes = bootstrap_label_nodes( $tree, \@trees );
1917 :     #
1918 :     # Option keys:
1919 :     #
1920 :     # consense => program # Define alternative name or path to executable
1921 : golsen 1.3 # midpoint => boolean # Midpoint root reference tree before labeling
1922 : golsen 1.2 # savetmp => boolean # Do not delete temporary directory or files.
1923 :     # tmpdir => directory # Directory for temporary files. If directory
1924 :     # # exists, temporary files are not deleted.
1925 :     # tmp => directory # Directory for tmpdir (D = SEED tmp or /tmp)
1926 :     #
1927 :     #-------------------------------------------------------------------------------
1928 :     sub bootstrap_label_nodes
1929 :     {
1930 :     my $reftree = shift;
1931 :     my $options = ( ( @_ > 1 ) && ( ref( $_[-1] ) eq 'HASH' ) ) ? pop @_ : {};
1932 :    
1933 :     return undef if ( ! @_ ) || ( ref( $_[0] ) ne 'ARRAY' );
1934 :    
1935 : golsen 1.3 my @trees = ( @_ > 1 ) ? @_ : @{ $_[0] };
1936 : golsen 1.2 my $n_tree = @trees;
1937 : golsen 1.3 my $tip = 'tip000001';
1938 :     my %tips = map { $_ => $tip++ } gjonewicklib::newick_tip_list( $trees[0] );
1939 : golsen 1.2
1940 :     # Make a copy of trees with clean labels
1941 :    
1942 :     my @trees2 = map { gjonewicklib::newick_relabel_tips( gjonewicklib::copy_newick_tree( $_ ), \%tips ) }
1943 :     @trees;
1944 :    
1945 : golsen 1.3 my $program = SeedAware::executable_for( $options->{ consense } || $options->{ program } || 'consense' );
1946 : golsen 1.2
1947 : golsen 1.3 my ( $tmp_dir, $save_tmp ) = SeedAware::temporary_directory( $options );
1948 : golsen 1.2 $tmp_dir or return ();
1949 :    
1950 :     #---------------------------------------------------------------------------
1951 :     # Write the files and run the program:
1952 :     #---------------------------------------------------------------------------
1953 :    
1954 :     my $cwd = $ENV{ cwd } || `pwd`;
1955 : golsen 1.3
1956 : golsen 1.2 chomp $cwd;
1957 :     chdir $tmp_dir;
1958 :     my $data = { trees => \@trees2, consense => $program };
1959 :     my ( undef, $id_index, $bs_frac ) = run_consense( $data );
1960 : golsen 1.3
1961 : golsen 1.2 chdir $cwd;
1962 :    
1963 :     # Delete the temporary directory unless it already existed:
1964 :    
1965 : golsen 1.3 system( '/bin/rm', '-r', $tmp_dir ) if ! $save_tmp;
1966 : golsen 1.2
1967 :     $id_index && $bs_frac or return undef;
1968 :    
1969 :     # We need to map taxon labels and then traverse the tree, putting the
1970 :     # bootstrap fractions as internal node labels:
1971 :    
1972 :     my %id_index_2 = map { $_ => $id_index->{ $tips{ $_ } } }
1973 :     keys %tips;
1974 :    
1975 : golsen 1.3 my $format = $n_tree > 100 ? '%.3f' : '%.2f';
1976 : golsen 1.2
1977 :     my ( $set1, $set2, $frac );
1978 :     my %bs_frac_2 = map { ( $set1 = $_ ) =~ tr/.*/01/;
1979 :     ( $set2 = $set1 ) =~ tr/01/10/;
1980 :     $frac = sprintf( $format, $bs_frac->{ $_ } );
1981 :     ( $set1 => $frac, $set2 => $frac )
1982 :     }
1983 :     keys %$bs_frac;
1984 :    
1985 : golsen 1.3 $reftree = gjonewicklib::copy_newick_tree( $reftree ); # Copy reference tree
1986 :     $reftree = gjonewicklib::reroot_newick_to_midpoint_w( $reftree ) if $options->{ midpoint };
1987 :    
1988 : golsen 1.2 label_nodes_with_bs( $reftree, \%id_index_2, \%bs_frac_2 );
1989 :     gjonewicklib::set_newick_lbl( $reftree, undef ); # Root node label
1990 :    
1991 :     $reftree;
1992 :     }
1993 :    
1994 :    
1995 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1996 :     # Label tree nodes with their bootstrap proportion. Work from the tips to
1997 :     # root. Return the "tip_set" representing the set of tips below the node.
1998 :     #
1999 :     # $tip_set = label_nodes_with_bs( $reftree, \%id_index, $bs_frac );
2000 :     #
2001 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2002 :     sub label_nodes_with_bs
2003 :     {
2004 :     my ( $node, $id_index, $bs_frac ) = @_;
2005 :    
2006 :     my $tip_set = '0' x keys %$id_index;
2007 :    
2008 :     my @dl = gjonewicklib::newick_desc_list( $node );
2009 :     if ( @dl )
2010 :     {
2011 :     foreach ( @dl )
2012 :     {
2013 :     $tip_set |= label_nodes_with_bs( $_, $id_index, $bs_frac );
2014 :     }
2015 :     gjonewicklib::set_newick_lbl( $node, $bs_frac->{ $tip_set } || '0.00' );
2016 :     }
2017 :     else
2018 :     {
2019 :     my $index = $id_index->{ gjonewicklib::newick_lbl( $node ) };
2020 :     substr( $tip_set, $index, 1 ) = "1";
2021 :     }
2022 :    
2023 :     $tip_set;
2024 :     }
2025 :    
2026 :    
2027 :     #-------------------------------------------------------------------------------
2028 :     # Internal routine to run consense:
2029 :     #
2030 :     # ( $tree, \%id_index, \%boot_frac ) = run_consense( \%data )
2031 :     #
2032 :     # Data:
2033 :     #
2034 :     # consense => $program_name_or_path # D = consense
2035 :     # trees => \@trees
2036 :     #
2037 :     #-------------------------------------------------------------------------------
2038 :     sub run_consense
2039 :     {
2040 :     my $data = shift;
2041 :     $data and ref( $data ) eq 'HASH'
2042 :     or print STDERR "gjophylip::run_consense: Called without data hash.\n"
2043 :     and return undef;
2044 :    
2045 :     $data->{trees} && ref( $data->{trees} ) eq 'ARRAY'
2046 :     or print STDERR "gjophylip::run_consense: Array of trees not provided.\n"
2047 :     and return undef;
2048 :    
2049 :     my $n_tree = @{ $data->{trees} }
2050 :     or print STDERR "gjophylip::run_consense: No trees.\n"
2051 :     and return undef;
2052 :    
2053 :     write_intree_no_count( @{ $data->{trees} } );
2054 :    
2055 :     unlink 'outfile' if -f 'outfile'; # Just checking
2056 :     unlink 'outtree' if -f 'outtree'; # ditto
2057 :    
2058 :    
2059 :     # No command options, so just echo Y to exectute
2060 :    
2061 : golsen 1.3 my $program = $data->{ program } || $data->{ consense } || 'consense';
2062 :     my $redirects = { stdout => 'consense_out', # Gets deleted with tmp dir
2063 :     stderr => 'consense_err' # Gets deleted with tmp dir
2064 :     };
2065 :     my $consenseFH = SeedAware::write_to_pipe_with_redirect( $program, $redirects )
2066 :     or return ();
2067 :     print $consenseFH "Y\n";
2068 :     close $consenseFH;
2069 : golsen 1.2
2070 :     my $tree = gjonewicklib::read_newick_tree( 'outtree' );
2071 :    
2072 : golsen 1.3 $tree or print STDERR "gjophylip::run_consense: Could not read consense outtree file.\n"
2073 : golsen 1.2 and return ();
2074 :    
2075 :     $tree = gjonewicklib::uproot_newick( $tree );
2076 :     $tree = gjonewicklib::newick_modify_branches( $tree, \&rescale_bootstrap, [ 1/$n_tree ] );
2077 :    
2078 :     my ( $id_index, $bs_frac ) = read_consense_partitians();
2079 :    
2080 : golsen 1.3 $id_index && $bs_frac
2081 :     or print STDERR "gjophylip::run_consense: Could not read consense outfile.\n"
2082 :     and return ();
2083 : golsen 1.2
2084 :     return ( $tree, $id_index, $bs_frac );
2085 :     }
2086 :    
2087 :    
2088 :     sub rescale_bootstrap
2089 :     {
2090 :     my( $x, $scale ) = @_;
2091 :     return undef if ! defined $x;
2092 :     $x *= $scale;
2093 : golsen 1.3 $x = 1 if $x > 1; # This deals with an artifact in some PHYLIP consensus trees
2094 : golsen 1.2 return $x;
2095 :     }
2096 :    
2097 :    
2098 :     #-------------------------------------------------------------------------------
2099 : golsen 1.1 # Utility functions:
2100 :     #-------------------------------------------------------------------------------
2101 :     # Allow variations of option keys including uppercase, underscores and
2102 :     # terminal 's':
2103 :     #
2104 :     # $key = canonical_key( $key );
2105 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2106 :     sub canonical_key
2107 :     {
2108 :     my $key = lc shift;
2109 :     $key =~ s/_//g;
2110 :     $key =~ s/s$//; # This is dangerous if an s is part of a word!
2111 :     return $key;
2112 :     }
2113 :    
2114 :    
2115 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2116 : golsen 1.3 # min
2117 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2118 :     sub min { $_[0] < $_[1] ? $_[0] : $_[1] }
2119 :    
2120 :    
2121 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2122 :     # pack_by_mask
2123 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2124 : golsen 1.3 sub pack_by_mask
2125 : golsen 1.1 {
2126 : golsen 1.3 my ( $seq, $mask ) = @_;
2127 :     $seq &= $mask; # Mask the string
2128 :     $seq =~ tr/\000//d; # Compress out X'00'
2129 :     $seq;
2130 :     }
2131 :    
2132 : golsen 1.1
2133 : golsen 1.3 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2134 :     # process_dna_alignment
2135 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2136 :     sub process_dna_alignment
2137 :     {
2138 :     my ( $align ) = @_;
2139 : golsen 1.1
2140 : golsen 1.3 if ( ! $align || ( ref( $align ) ne 'ARRAY' ) || ( @$align < 2 ) )
2141 : golsen 1.1 {
2142 : golsen 1.3 print STDERR "gjophylip function called without a required alignment.\n";
2143 :     return ();
2144 : golsen 1.1 }
2145 : golsen 1.3
2146 :     # Make a clean copy of the alignment. Id is always first, seq is always last.
2147 :    
2148 :     my ( $seq, $id, %id, %local_id );
2149 :     my $local_id = 'seq0000000';
2150 :     my @align = map { $id = $_->[0];
2151 :     $local_id++;
2152 :     $id{ $local_id } = $id;
2153 :     $local_id{ $id } = $local_id;
2154 :     $seq = $_->[-1];
2155 :     $seq =~ s/[^ABCDGHKMNRSTUVWY]/N/gi; # Bad letters go to N
2156 :     $seq =~ s/[^A-Z]/-/gi; # Anything else becomes -
2157 :     [ $local_id, '', $seq ]
2158 :     }
2159 :     @$align;
2160 :    
2161 :     # Should probably add a length check
2162 :    
2163 :     return ( \@align, \%id, \%local_id );
2164 :     }
2165 :    
2166 :    
2167 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2168 :     # process_protein_alignment
2169 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2170 :     sub process_protein_alignment
2171 :     {
2172 :     my ( $align ) = @_;
2173 :    
2174 :     if ( ! $align || ( ref( $align ) ne 'ARRAY' ) || ( @$align < 2 ) )
2175 : golsen 1.1 {
2176 : golsen 1.3 print STDERR "gjophylip function called without a required alignment.\n";
2177 :     return ();
2178 : golsen 1.1 }
2179 :    
2180 : golsen 1.3 # Make a clean copy of the alignment. Id is always first, seq is always last.
2181 :    
2182 :     my ( $seq, $id, %id, %local_id );
2183 :     my $local_id = 'seq0000000';
2184 :     my @align = map { $id = $_->[0];
2185 :     $local_id++;
2186 :     $id{ $local_id } = $id;
2187 :     $local_id{ $id } = $local_id;
2188 :     $seq = $_->[-1];
2189 :     $seq =~ s/[BJOUZ]/X/gi; # Bad letters go to X
2190 :     $seq =~ s/[^A-Z]/-/gi; # Anything else becomes -
2191 :     [ $local_id, '', $seq ]
2192 :     }
2193 :     @$align;
2194 :    
2195 :     # Should probably add a length check
2196 : golsen 1.1
2197 : golsen 1.3 return ( \@align, \%id, \%local_id );
2198 : golsen 1.1 }
2199 :    
2200 :    
2201 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2202 : golsen 1.3 # process_categories
2203 : golsen 1.2 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2204 : golsen 1.3 sub process_categories
2205 : golsen 1.2 {
2206 : golsen 1.3 my ( $categories ) = @_;
2207 :    
2208 :     # Verify the list:
2209 : golsen 1.2
2210 : golsen 1.3 if ( ref( $categories ) ne 'ARRAY'
2211 :     || ( @$categories != 2 )
2212 :     || ( ref( $categories->[0] ) ne 'ARRAY' )
2213 :     || ( @{ $categories->[0] } < 2 )
2214 :     || ( @{ $categories->[0] } > 9 )
2215 :     )
2216 :     {
2217 :     print STDERR "gjophylip categories option value must be [ [ rate1, ... ], site_cats ]\n";
2218 :     print STDERR " with 2 - 9 rate categories\n";
2219 :     return undef;
2220 :     }
2221 :    
2222 :     # Rate values cannot have many decimal places or proml can't read them:
2223 : golsen 1.2
2224 : golsen 1.3 @{ $categories->[0] } = map { sprintf( "%.6f", $_ ) } @{ $categories->[0] };
2225 :    
2226 :     return $categories;
2227 : golsen 1.2 }
2228 :    
2229 :    
2230 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2231 : golsen 1.3 # process_distance_matrix
2232 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2233 : golsen 1.3 sub process_distance_matrix
2234 :     {
2235 :     my ( $distances ) = @_;
2236 :    
2237 :     if ( ! $distances || ( ref( $distances ) ne 'ARRAY' ) || ( @$distances < 2 ) )
2238 :     {
2239 :     print STDERR "gjophylip function called without required distances.\n";
2240 :     return ();
2241 :     }
2242 :    
2243 :     # Make a clean copy of the distances with a local id.
2244 :    
2245 :     my ( $id, %id, %local_id );
2246 :     my $local_id = 'tax0000000';
2247 :     my @distances = map { $id = $_->[0];
2248 :     $local_id++;
2249 :     $id{ $local_id } = $id;
2250 :     $local_id{ $id } = $local_id;
2251 :     [ $local_id, @$_[ 1 .. $#{@$_} ] ]
2252 :     }
2253 :     @$distances;
2254 :    
2255 :     # Should probably add a length check
2256 :    
2257 :     return ( \@distances, \%id, \%local_id );
2258 :     }
2259 : golsen 1.1
2260 :    
2261 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2262 : golsen 1.3 # process_dna_model
2263 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2264 : golsen 1.3 sub process_dna_model
2265 : golsen 1.1 {
2266 : golsen 1.3 my ( $options ) = @_;
2267 :    
2268 :     $options && ref $options eq 'HASH' && $options->{ model }
2269 :     or return 'F84';
2270 :    
2271 :     # Allow a ridiculous list of synonyms:
2272 :    
2273 :     my $model = ( $options->{ model } =~ m/F84/i ) ? 'F84'
2274 :     : ( $options->{ model } =~ m/Churchill/i ) ? 'F84'
2275 :     : ( $options->{ model } =~ m/Felsenstein/i ) ? 'F84'
2276 :     : ( $options->{ model } =~ m/Hasegawa/i ) ? 'F84'
2277 :     : ( $options->{ model } =~ m/Kishino/i ) ? 'F84'
2278 :    
2279 :     : ( $options->{ model } =~ m/Kimura/i ) ? 'Kimura'
2280 :    
2281 :     : ( $options->{ model } =~ m/JC/i ) ? 'JC'
2282 :     : ( $options->{ model } =~ m/Jukes/i ) ? 'JC'
2283 :     : ( $options->{ model } =~ m/Cantor/i ) ? 'JC'
2284 :    
2285 :     : ( $options->{ model } =~ m/LogDet/i ) ? 'LogDet'
2286 :     : ( $options->{ model } =~ m/Barry/i ) ? 'LogDet'
2287 :     : ( $options->{ model } =~ m/Hartigan/i ) ? 'LogDet'
2288 :     : ( $options->{ model } =~ m/Lake/i ) ? 'LogDet'
2289 :     : ( $options->{ model } =~ m/Lockhart/i ) ? 'LogDet'
2290 :     : ( $options->{ model } =~ m/Steel/i ) ? 'LogDet'
2291 :    
2292 :     : ( $options->{ model } =~ m/identity/i ) ? 'identity'
2293 :     : ( $options->{ model } =~ m/similarity/i ) ? 'identity'
2294 :    
2295 :     : 'F84';
2296 :    
2297 :     return $model;
2298 : golsen 1.1 }
2299 :    
2300 :    
2301 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2302 : golsen 1.3 # process_protein_model
2303 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2304 : golsen 1.3 sub process_protein_model
2305 : golsen 1.1 {
2306 : golsen 1.3 my ( $options ) = @_;
2307 :    
2308 :     $options && ref $options eq 'HASH' && $options->{ model }
2309 :     or return 'JTT';
2310 :    
2311 :     # Allow a ridiculous list of synonyms:
2312 :    
2313 :     my $model = ( $options->{ model } =~ m/PAM/i ) ? 'PAM'
2314 :     : ( $options->{ model } =~ m/Dayhoff/i ) ? 'PAM'
2315 :     : ( $options->{ model } =~ m/Kosiol/i ) ? 'PAM'
2316 :     : ( $options->{ model } =~ m/Goldman/i ) ? 'PAM'
2317 :     : ( $options->{ model } =~ m/DCMut/i ) ? 'PAM'
2318 :    
2319 :     : ( $options->{ model } =~ m/PMB/i ) ? 'PMB'
2320 :     : ( $options->{ model } =~ m/Henikoff/i ) ? 'PMB'
2321 :     : ( $options->{ model } =~ m/Smith/i ) ? 'PMB'
2322 :     : ( $options->{ model } =~ m/Tillier/i ) ? 'PMB'
2323 :     : ( $options->{ model } =~ m/Veerassamy/i ) ? 'PMB'
2324 :    
2325 :     : ( $options->{ model } =~ m/JTT/i ) ? 'JTT'
2326 :     : ( $options->{ model } =~ m/Jones/i ) ? 'JTT'
2327 :     : ( $options->{ model } =~ m/Taylor/i ) ? 'JTT'
2328 :     : ( $options->{ model } =~ m/Thornton/i ) ? 'JTT'
2329 :    
2330 :     : 'JTT';
2331 :    
2332 :     return $model;
2333 : golsen 1.1 }
2334 :    
2335 :    
2336 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2337 : golsen 1.3 # process_tree_format( \%options )
2338 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2339 : golsen 1.3 sub process_tree_format
2340 : golsen 1.1 {
2341 : golsen 1.3 my ( $options ) = @_;
2342 :    
2343 :     $options && ref $options eq 'HASH' or return 'overbeek';
2344 :    
2345 :     # Allow a ridiculous list of synonyms:
2346 :    
2347 :     return ! $options->{ treeformat } ? 'overbeek' :
2348 :     $options->{ treeformat } =~ m/overbeek/i ? 'overbeek' :
2349 :     $options->{ treeformat } =~ m/gjo/i ? 'gjonewick' :
2350 :     $options->{ treeformat } =~ m/fig/i ? 'fig' :
2351 :     'overbeek'; # Default
2352 : golsen 1.1 }
2353 :    
2354 :    
2355 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2356 : golsen 1.3 # process_rate_hmm
2357 : golsen 1.2 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2358 : golsen 1.3 sub process_rate_hmm
2359 : golsen 1.2 {
2360 : golsen 1.3 my ( $rate_hmm ) = @_;
2361 :    
2362 :     # Verify the list:
2363 :    
2364 :     if ( ( ref( $rate_hmm ) ne 'ARRAY' )
2365 :     || ( @$rate_hmm < 2 )
2366 :     || ( @$rate_hmm > 9 )
2367 :     )
2368 :     {
2369 :     print STDERR "gjophylip rate_hmm must be a reference to 2-9 rate-probability pairs\n";
2370 :     return undef;
2371 :     }
2372 :    
2373 :     # Verify the elements of the list:
2374 :    
2375 :     my $total = 0;
2376 :     my @rate_hmm2 = grep { ( ref( $_ ) eq 'ARRAY' )
2377 :     && ( @$_ == 2 ) # [ rate, probability ]
2378 :     && ( $_->[0] >= 0 ) # Rate can be zero
2379 :     && ( $_->[1] = 0 ) # Probability cannot be zero
2380 :     && ( $total += $_->[1] ) # Cumumlative probability
2381 :     }
2382 :     @$rate_hmm;
2383 :    
2384 :     # If anything was lost, there was an error:
2385 :    
2386 :     ( @rate_hmm2 == @$rate_hmm )
2387 :     or print STDERR "gjophylip rate_hmm must be a reference to rate-probability pairs\n"
2388 :     and return undef;
2389 :    
2390 :     # Normalize the probabilities:
2391 :    
2392 :     foreach ( @rate_hmm2 ) { $_->[1] /= $total }
2393 :    
2394 :     return \@rate_hmm2;
2395 : golsen 1.2 }
2396 :    
2397 :    
2398 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2399 : golsen 1.3 # write_categories( $categories )
2400 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2401 :     sub write_categories
2402 :     {
2403 :     my $categories = shift;
2404 :     open( CATEGORIES, '>categories' ) or return 0;
2405 :     print CATEGORIES "$categories\n";
2406 :     close( CATEGORIES );
2407 :     }
2408 :    
2409 :    
2410 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2411 : golsen 1.3 # write_dist_infile( $distances )
2412 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2413 :     sub write_dist_infile
2414 :     {
2415 :     my $distances = shift;
2416 :     ( ref( $distances ) eq 'ARRAY' ) && @$distances or return 0;
2417 :     my $ntaxa = @$distances;
2418 :     my $format = "%-10s" . (" %11.6f" x $ntaxa) . "\n";
2419 :     open( DISTANCES, '>infile' ) or return 0;
2420 :     print DISTANCES "$ntaxa\n";
2421 :     foreach ( @$distances ) { printf DISTANCES $format, @$_ }
2422 :     close( DISTANCES );
2423 :     }
2424 :    
2425 :    
2426 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2427 : golsen 1.3 # write_intree( @trees )
2428 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2429 :     sub write_intree
2430 :     {
2431 :     open( INTREE, '>intree' ) or return 0;
2432 :     print INTREE scalar @_, "\n";
2433 :     foreach ( @_ ) { print INTREE gjonewicklib::strNewickTree( $_ ), "\n" }
2434 :     close( INTREE );
2435 :     }
2436 :    
2437 :    
2438 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2439 :     # write_intree_no_count
2440 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2441 :     sub write_intree_no_count
2442 :     {
2443 :     open( INTREE, '>intree' ) or return 0;
2444 :     foreach ( @_ ) { print INTREE gjonewicklib::strNewickTree( $_ ), "\n" }
2445 :     close( INTREE );
2446 :     }
2447 :    
2448 :    
2449 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2450 :     # write_seq_infile( @alignment )
2451 :     # write_seq_infile( \@alignment )
2452 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2453 :     sub write_seq_infile
2454 :     {
2455 :     @_ && $_[0] && ref $_[0] eq 'ARRAY' && @{$_[0]} && defined $_[0]->[0]
2456 :     or return undef;
2457 :     my @seq = ref $_[0]->[0] eq 'ARRAY' ? @{ $_[0] } : @_;
2458 :    
2459 :     open( INFILE, '>infile' ) or return 0;
2460 :     printf INFILE "%d %d\n", scalar @_, length( $seq[0]->[-1] );
2461 :     foreach ( @seq ) { printf INFILE "%-10s %s\n", @$_[0, -1] }
2462 :     close( INFILE );
2463 :     }
2464 :    
2465 :    
2466 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2467 :     # write_weights( $weights )
2468 : golsen 1.1 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2469 :     sub write_weights
2470 :     {
2471 :     my $weights = shift;
2472 :     open( WEIGHTS, '>weights' ) or return 0;
2473 :     print WEIGHTS "$weights\n";
2474 :     close( WEIGHTS );
2475 :     }
2476 :    
2477 :    
2478 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2479 : golsen 1.3 # read_categories
2480 :     #
2481 :     # Expect category rates on first line, followed by site categories.
2482 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2483 :     sub read_categories
2484 :     {
2485 :     $_[0] && -f $_[0]
2486 :     or print STDERR "read_categories() requires a file name.\n"
2487 :     and return undef;
2488 :     open( CATEGS, '<$_[0]' )
2489 :     or print STDERR "Could not open categories file $_[0]\n"
2490 :     and return undef;
2491 :     my @rates = split " ", <CATEGS>;
2492 :     my $categs = join( '', map { chomp; s/ \s+//g; $_ }
2493 :     <CATEGS>
2494 :     );
2495 :     close( CATEGS );
2496 :     return [ \@rates, $categs ];
2497 :     }
2498 :    
2499 :    
2500 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2501 :     # read_weights
2502 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2503 :     sub read_weights
2504 :     {
2505 :     $_[0] && -f $_[0]
2506 :     or print STDERR "read_weights() requires a file name.\n"
2507 :     and return undef;
2508 :     open( WEIGHTS, '<$_[0]' )
2509 :     or print STDERR "Could not open weights file $_[0]\n"
2510 :     and return undef;
2511 :     my $weights = join( '', map { chomp; s/ \s+//g; $_ } <WEIGHTS> );
2512 :     close( WEIGHTS );
2513 :    
2514 :     return $weights;
2515 :     }
2516 :    
2517 :    
2518 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2519 : golsen 1.1 # read_distances
2520 :     #
2521 : golsen 1.3 # @distance_matrix = read_distances() # D = outfile
2522 :     # \@distance_matrix = read_distances()
2523 :     # @distance_matrix = read_distances( $file )
2524 :     # \@distance_matrix = read_distances( $file )
2525 :     #
2526 : golsen 1.1 # @distance_matrix = ( [ id1, dist11, dist12, dist13, ... ],
2527 :     # [ id2, dist21, dist22, dist14, ... ],
2528 :     # ...
2529 :     # )
2530 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2531 :     sub read_distances
2532 :     {
2533 : golsen 1.3 my $file = $_[0] && -f $_[0] ? $_[0] : 'outfile';
2534 :    
2535 : golsen 1.1 my @distances;
2536 : golsen 1.3 open( DISTS, "<$file" ) or return ();
2537 :     my $line;
2538 :     defined( $line = <DISTS> ) or return ();
2539 : golsen 1.1 my ( $ndist ) = $line =~ m/(\d+)/;
2540 :    
2541 :     for ( my $i = 0; $i < $ndist; $i++ )
2542 :     {
2543 : golsen 1.3 defined( $line = <DISTS> ) or return ();
2544 :     chomp $line;
2545 : golsen 1.1 my ( @row ) = split " ", $line;
2546 :     while ( @row <= $ndist )
2547 :     {
2548 : golsen 1.3 defined( $line = <DISTS> ) or return ();
2549 :     chomp $line;
2550 : golsen 1.1 push @row, split " ", $line;
2551 :     }
2552 :    
2553 :     push @distances, \@row;
2554 :     }
2555 :     close( DISTS );
2556 :    
2557 :     wantarray ? @distances : \@distances;
2558 :     }
2559 :    
2560 :    
2561 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2562 :     # read_likelihoods
2563 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2564 :     sub read_likelihoods
2565 :     {
2566 : golsen 1.3 my $file = $_[0] && -f $_[0] ? $_[0] : 'outfile';
2567 :    
2568 :     open( OUTFILE, "<$file" ) or return ();
2569 : golsen 1.1 my @likelihoods = map { chomp; s/.* //; $_ }
2570 :     grep { /^Ln Likelihood/ }
2571 :     <OUTFILE>;
2572 :     close( OUTFILE );
2573 : golsen 1.3
2574 :     wantarray ? @likelihoods : \@likelihoods;
2575 : golsen 1.1 }
2576 :    
2577 :    
2578 : golsen 1.2 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2579 : golsen 1.3 # ( \%id_index, \%bs_frac ) = read_consense_partitians()
2580 :     # ( \%id_index, \%bs_frac ) = read_consense_partitians( $file )
2581 : golsen 1.2 #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2582 :     sub read_consense_partitians
2583 :     {
2584 : golsen 1.3 my $file = defined $_[0] && -f $_[0] ? $_[0] : 'outfile';
2585 :     open( OUTFILE, "<$file" ) or return ();
2586 : golsen 1.2
2587 :     local $_;
2588 : golsen 1.3 while ( defined( $_ = <OUTFILE> ) && ! /^Species in order/ ) {}
2589 : golsen 1.2 $_ or close( OUTFILE ) and return ();
2590 :    
2591 :     # Hash with zero-based offset of each taxon in the partitian set strings
2592 :    
2593 :     my %id_index;
2594 : golsen 1.3 while ( defined( $_ = <OUTFILE> ) && ! /^Sets included in the consensus tree/ )
2595 : golsen 1.2 {
2596 :     chomp;
2597 :     my ( $n, $id ) = split;
2598 :     $id_index{ $id } = $n - 1 if $n && $id;
2599 :     }
2600 :     $_ or close( OUTFILE ) and return ();
2601 :    
2602 : golsen 1.3 while ( defined( $_ = <OUTFILE> ) && ! /How many times out of \s*\d\S*\d/ ) {}
2603 : golsen 1.2 $_ or close( OUTFILE ) and return ();
2604 :     my( $ttl ) = /How many times out of \s*(\d\S*\d)/;
2605 :     $ttl or close( OUTFILE ) and return ();
2606 :    
2607 :     my %bs_frac;
2608 : golsen 1.3 while ( defined( $_ = <OUTFILE> ) && ! /^Extended majority rule/ )
2609 : golsen 1.2 {
2610 :     next if /^Set/;
2611 :     chomp;
2612 :     my ( $set, $cnt ) = m/^(\S.*\S)\s+(\d+\.\d+)\s*$/;
2613 :     $set && $cnt or next;
2614 :     $set =~ s/\s+//g;
2615 :     $bs_frac{ $set } = $cnt / $ttl;
2616 :     }
2617 :     $_ or close( OUTFILE ) and return ();
2618 :    
2619 :     close( OUTFILE );
2620 :    
2621 :     return ( \%id_index, \%bs_frac );
2622 :     }
2623 :    
2624 :    
2625 : golsen 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3