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

Annotation of /FigKernelPackages/proml.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package proml;
2 :    
3 :     #===============================================================================
4 :     # A perl interface to the proml program in the PHYLIP program package
5 :     #
6 :     # @tree_likelihood_pairs = proml( \@alignment, \%options )
7 :     # @tree_likelihood_pairs = proml( \@alignment, %options )
8 :     # @tree_likelihood_pairs = proml( \%options ) # alignment must be included as option
9 :     # @tree_likelihood_pairs = proml( %options ) # alignment must be included as option
10 :     #
11 :     # @alignment = array of id_seq pairs, or id_definition_seq triples
12 :     #
13 :     #-------------------------------------------------------------------------------
14 :     # A perl interface for using proml to estimate site-specific rates of change
15 :     #
16 :     # ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, proml_opts )
17 :     #
18 :     # $categories = [ $n_categories, [ $rates ], $site_categories ];
19 :     #
20 :     #===============================================================================
21 :     #
22 :     # A perl interface to the proml program in the PHYLIP program package
23 :     #
24 :     # @tree_likelihood_pairs = proml( \@alignment, \%options )
25 :     # @tree_likelihood_pairs = proml( \@alignment, %options )
26 :     # @tree_likelihood_pairs = proml( \%options ) # alignment must be included as option
27 :     # @tree_likelihood_pairs = proml( %options ) # alignment must be included as option
28 :     #
29 :     # @alignment = array of id_seq pairs, or id_definition_seq triples
30 :     #
31 :     # options:
32 :     #
33 :     # For proml:
34 :     # alignment => \@alignment the way to supply the alignment as an option, rather than first param
35 :     # alpha => float alpha parameter of gamma distribution (0.5 - inf)
36 :     # categories => [ n, [ rates ], site_categories ]
37 :     # coef_of_var => float 1/sqrt(alpha) for gamma distribution (D = 0)
38 :     # gamma_bins => int number of rate categories used to approximate gamma (D=5)
39 :     # global => bool global rearrangements
40 :     # invar_frac => 0 - 1 fraction of site that are invariant
41 :     # jumble_seed => odd int jumble random seed
42 :     # model => model evolution model JTT (D) | PMB | PAM
43 :     # n_jumble => int number of jumbles
44 :     # persistance => float persistance length of rate category
45 :     # rate_hmm => [ n_rates, [ rates ], [ probabilies ] ]
46 :     # rearrange => [ trees ] rearrange user trees
47 :     # slow => bool more accurate but slower search (D = 0)
48 :     # user_lengths => bool use supplied branch lengths
49 :     # user_trees => [ trees ] user trees
50 :     # weights => site_weights
51 :     #
52 :     # Other:
53 :     # keep_duplicates => bool do not remove duplicate sequences (D = false) [NOT IMPLIMENTED]
54 :     # program => program allows fully defined path
55 :     # tmp => directory directory for tmp_dir (D = /tmp or .)
56 :     # tmp_dir => directory directory for temporary files (D = $tmp/proml.$$)
57 :     # tree_format => overbeek | gjo | fig format of output tree
58 :     #
59 :     # tmp_dir is created and deleted unless its name is supplied, and it already
60 :     # exists.
61 :     #
62 :     #
63 :     # Options that do not require other data:
64 :     # G (global search toggle)
65 :     # L (user lengths toggle)
66 :     # P (JTT / PAM toggle)
67 :     # S (slow and accurate)
68 :     # U (requires intree file)
69 :     # W (requires weights file)
70 :     #
71 :     # Some option data input orders:
72 :     #
73 :     # J
74 :     # Seed
75 :     # N reps
76 :     # Y
77 :     #
78 :     # R
79 :     # Y
80 :     # Coefficient of variation
81 :     # Rate categories
82 :     # Spurious random seed
83 :     #
84 :     # R
85 :     # R
86 :     # Y
87 :     # Coefficient of variation
88 :     # Gamma rate categories + 1
89 :     # Fraction invariant
90 :     # Spurious random seed
91 :     #
92 :     # C (requires categories file)
93 :     # N cat
94 :     # Rate values (n of them)
95 :    
96 :    
97 :     use strict;
98 :     use gjonewicklib qw( gjonewick_to_overbeek
99 :     newick_is_unrooted
100 :     newick_relabel_nodes
101 :     newick_tree_length
102 :     overbeek_to_gjonewick
103 :     parse_newick_tree_str
104 :     strNewickTree
105 :     uproot_newick
106 :     );
107 :    
108 :    
109 :     sub proml
110 :     {
111 :     my $align;
112 :     if ( ref( $_[0] ) eq 'ARRAY' )
113 :     {
114 :     $align = shift @_;
115 :     ( $align && ( ref( $align ) eq 'ARRAY' ) )
116 :     || ( ( print STDERR "proml::proml() called without alignment\n" )
117 :     && ( return () )
118 :     );
119 :     }
120 :    
121 :     my %options;
122 :     if ( $_[0] )
123 :     {
124 :     %options = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
125 :     }
126 :    
127 :     #---------------------------------------------------------------------------
128 :     # Work on a copy of the alignment. Id is always first, seq is always last
129 :     #---------------------------------------------------------------------------
130 :    
131 :     $align ||= $options{ alignment } || $options{ align };
132 :    
133 :     my ( $seq, $id );
134 :     my %id;
135 :     my %local_id;
136 :     my $local_id = 'seq0000000';
137 :     my @align = map { $id = $_->[0];
138 :     $local_id++;
139 :     $id{ $local_id } = $id;
140 :     $local_id{ $id } = $local_id;
141 :     $seq = $_->[-1];
142 :     $seq =~ s/[BJOUZ]/X/gi; # Bad letters go to X
143 :     $seq =~ s/[^A-Z]/-/gi; # Anything else becomes -
144 :     [ $local_id, $seq ]
145 :     } @$align;
146 :    
147 :     #---------------------------------------------------------------------------
148 :     # Process proml options:
149 :     #---------------------------------------------------------------------------
150 :    
151 :     my $categories = $options{ categories }; # [ n_cat, [ cat_rates ], site_cats ]
152 :     if ( $categories )
153 :     {
154 :     if ( ref( $categories ) ne 'ARRAY'
155 :     || @$categories != 3
156 :     || $categories->[0] < 1 || $categories->[0] > 9
157 :     || ref( $categories->[1] ) ne 'ARRAY'
158 :     || @{$categories->[1]} != $categories->[0]
159 :     )
160 :     {
161 :     print STDERR "proml::proml categories option value must be [ n_cat, [ cat_rates ], site_cats ]\n";
162 :     return ();
163 :     }
164 :    
165 :     # Rate values cannot have very many decimal places or proml can't read it:
166 :    
167 :     @{$categories->[1]} = map { sprintf "%.6f", $_ } @{$categories->[1]};
168 :     }
169 :    
170 :     my $coef_of_var = $options{ coef_of_var }
171 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
172 :     || 0;
173 :     if ( $coef_of_var < 0 )
174 :     {
175 :     print STDERR "proml::proml coef_of_var option value must be >= 0\n";
176 :     return ();
177 :     }
178 :    
179 :     my $gamma_bins = int( $options{ gamma_bins } || ( $coef_of_var ? 5 : 2 ) );
180 :     if ( ( $gamma_bins < 2 ) || ( $gamma_bins > 9 ) )
181 :     {
182 :     print STDERR "proml::proml gamma_bins option value must be > 1 and <= 9\n";
183 :     return ();
184 :     }
185 :    
186 :     my $global = $options{ global } || 0;
187 :    
188 :     my $invar_frac = $options{ invar_frac } || 0;
189 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
190 :     {
191 :     print STDERR "proml::proml invar_frac option value must be >= 0 and < 1\n";
192 :     return ();
193 :     }
194 :    
195 :     my $n_jumble = int( $options{ n_jumble } || ( $options{ jumble_seed } ? 1 : 0) );
196 :     if ( $n_jumble < 0 )
197 :     {
198 :     print STDERR "proml::proml n_jumble option value must be >= 0\n";
199 :     return ();
200 :     }
201 :    
202 :     my $jumble_seed = int( $options{ jumble_seed } || 4 * int( 499999999 * rand() ) + 1 );
203 :     if ( ( $jumble_seed <= 0) || ( $jumble_seed % 2 != 1 ) )
204 :     {
205 :     print STDERR "proml::proml jumble_seed option value must be an odd number > 0\n";
206 :     return ();
207 :     }
208 :    
209 :     my $model = ( $options{ model } =~ m/PAM/i ) ? 'PAM'
210 :     : ( $options{ model } =~ m/Dayhoff/i ) ? 'PAM'
211 :     : ( $options{ model } =~ m/PMB/i ) ? 'PMB'
212 :     : ( $options{ model } =~ m/Henikoff/i ) ? 'PMB'
213 :     : ( $options{ model } =~ m/Tillier/i ) ? 'PMB'
214 :     : ( $options{ model } =~ m/JTT/i ) ? 'JTT'
215 :     : ( $options{ model } =~ m/Jones/i ) ? 'JTT'
216 :     : ( $options{ model } =~ m/Taylor/i ) ? 'JTT'
217 :     : ( $options{ model } =~ m/Thornton/i ) ? 'JTT'
218 :     : 'JTT';
219 :    
220 :     my $persistance = $options{ persistance } || 0;
221 :     if ( $persistance && ( $persistance <= 1 ) )
222 :     {
223 :     print STDERR "proml::proml persistance option value must be > 1\n";
224 :     return ();
225 :     }
226 :    
227 :     my $rearrange = $options{ rearrange };
228 :    
229 :     my $slow = $options{ slow };
230 :    
231 :     my $user_lengths = $options{ user_lengths };
232 :    
233 :     my $user_trees = $options{ user_trees } || $rearrange;
234 :    
235 :     if ( $user_trees )
236 :     {
237 :     if ( ( ref( $user_trees ) ne 'ARRAY' ) || ( ! @$user_trees ) )
238 :     {
239 :     $user_trees = undef; # No trees
240 :     }
241 :     elsif ( ref( $user_trees->[0] ) ne 'ARRAY' ) # First element not tree
242 :     {
243 :     print STDERR "proml::proml usertree or rearrange option value must be reference to list of trees\n";
244 :     return ();
245 :     }
246 :     }
247 :    
248 :     my $weights = $options{ weights };
249 :    
250 :    
251 :     #---------------------------------------------------------------------------
252 :     # Options that are not proml options per se:
253 :     #---------------------------------------------------------------------------
254 :    
255 :     my $program = $options{ program } || 'proml';
256 :    
257 :     my $tmp = $options{ tmp };
258 :    
259 :     my $tmp_dir = $options{ tmp_dir };
260 :    
261 :     my $tree_format = $options{ tree_format } =~ m/overbeek/i ? 'overbeek'
262 :     : $options{ tree_format } =~ m/gjo/i ? 'gjonewick'
263 :     : $options{ tree_format } =~ m/fig/i ? 'fig'
264 :     : 'overbeek'; # Default
265 :    
266 :     my $save_tmp = $tmp_dir && -d $tmp_dir;
267 :     if ( $tmp_dir )
268 :     {
269 :     if ( -d $tmp_dir ) { $save_tmp = 1 }
270 :     else { mkdir $tmp_dir }
271 :     }
272 :     else
273 :     {
274 :     $tmp = $tmp && -d $tmp ? $tmp
275 :     : -d '/tmp' ? '/tmp'
276 :     : '.';
277 : overbeek 1.2 my $int = int( 1000000000 * rand);
278 :     $tmp_dir = "$tmp/proml.$$.$int";
279 : overbeek 1.1 mkdir $tmp_dir;
280 :     }
281 :    
282 :     #---------------------------------------------------------------------------
283 :     # Prepare data:
284 :     #---------------------------------------------------------------------------
285 :     #
286 :     # For simplicity, we will convert overbeek trees to gjo newick trees.
287 :     #
288 :     # gjonewick tree node: [ \@desc, $label, $x, \@c1, \@c2, \@c3, \@c4, \@c5 ]
289 :     #
290 :     # overbeek tree node: [ Label, DistanceToParent,
291 :     # [ ParentPointer, ChildPointer1, ... ],
292 :     # [ Name1\tVal1, Name2\tVal2, ... ]
293 :     # ]
294 :     # Root node of gjonewick always has a descendent list. If the first
295 :     # field of the first tree is not an array reference, they are overbeek
296 :     # trees. Also relabel tree tips to local ids.
297 :    
298 :     my @user_trees = ();
299 :     if ( $user_trees )
300 :     {
301 :     if ( @user_trees && ( ref( $user_trees[0]->[0] ) ne 'ARRAY' ) ) # overbeek trees
302 :     {
303 :     @user_trees = map { gjonewicklib::newick_relabel_nodes( $_, \%local_id ) }
304 :     map { gjonewicklib::overbeek_to_gjonewick( $_ ) }
305 :     @$user_trees;
306 :     }
307 :     else
308 :     {
309 :     @user_trees = map { gjonewicklib::newick_relabel_nodes( $_, \%local_id ) }
310 :     @$user_trees;
311 :     }
312 :    
313 :     # Make sure trees are unrooted:
314 :    
315 :     @user_trees = map { gjonewicklib::newick_is_unrooted( $_ ) ? $_
316 :     : gjonewicklib::uproot_newick( $_ )
317 :     }
318 :     @user_trees;
319 :     }
320 :    
321 :     #---------------------------------------------------------------------------
322 :     # Write the files and run the program:
323 :     #---------------------------------------------------------------------------
324 :    
325 :     my $cwd = $ENV{ cwd } || `pwd`;
326 :     chomp $cwd;
327 :     chdir $tmp_dir;
328 :    
329 :     unlink 'outfile' if -f 'outfile'; # Just checking
330 :     unlink 'outtree' if -f 'outtree'; # ditto
331 :    
332 :     &write_infile( @align ) or print STDERR "proml::proml: Could write infile\n"
333 :     and chdir $cwd
334 :     and return ();
335 :    
336 :     # open( PROML, "| $program > /dev/null" ) or print STDERR "proml::proml: Could not open pipe to $program\n"
337 :     open( PROML, ">params" ) or print STDERR "proml::proml: Could not open pipe to $program\n"
338 :     and chdir $cwd
339 :     and return ();
340 :    
341 :    
342 :     # Start sending optoins ot program:
343 :    
344 :     if ( $categories )
345 :     {
346 :     &write_categories( $categories->[2] ) or print STDERR "proml::proml: Could not write categories\n"
347 :     and chdir $cwd
348 :     and return ();
349 :     print PROML "C\n",
350 :     "$categories->[0]\n",
351 :     join( ' ', @{ $categories->[1] } ), "\n";
352 :     }
353 :    
354 :     if ( $invar_frac || $coef_of_var )
355 :     {
356 :     print PROML "R\n";
357 :     print PROML "R\n" if $invar_frac;
358 :     print PROML "A\n", "$persistance\n" if $persistance;
359 :    
360 :     }
361 :    
362 :     print PROML "G\n" if $global;
363 :    
364 :     print PROML "J\n", "$jumble_seed\n", "$n_jumble\n" if $n_jumble;
365 :    
366 :     print PROML "P\n" if $model =~ m/PMB/i;
367 :     print PROML "P\nP\n" if $model =~ m/PAM/i;
368 :    
369 :     print PROML "S\n" if $slow;
370 :    
371 :     if ( @user_trees )
372 :     {
373 :     &write_intree( @user_trees ) or print STDERR "proml::proml: Could write intree\n"
374 :     and chdir $cwd
375 :     and return ();
376 :     print PROML "U\n";
377 :     print PROML "V\n" if $rearrange || $global;
378 :     print PROML "L\n" if $user_lengths && ! $rearrange && ! $global;
379 :     }
380 :    
381 :     if ( $weights )
382 :     {
383 :     &write_weights( $weights ) or print STDERR "proml::proml: Could write weights\n"
384 :     and chdir $cwd
385 :     and return ();
386 :     print PROML "W\n";
387 :     }
388 :    
389 :     # All the options are sent, try to lauch the run:
390 :    
391 :     print PROML "Y\n";
392 :    
393 :     # Becuase of the options interface, these values have to be supplied after
394 :     # the Y:
395 :    
396 :     if ( $invar_frac || $coef_of_var )
397 :     {
398 :     if ( $invar_frac )
399 :     {
400 :     if ( $coef_of_var ) { $gamma_bins++ if ( $gamma_bins < 9 ) }
401 :     else { $gamma_bins = 2 }
402 :     }
403 :     print PROML "$coef_of_var\n";
404 :     print PROML "$gamma_bins\n";
405 :     print PROML "$invar_frac\n" if $invar_frac;
406 :     }
407 :     elsif ( $user_trees )
408 :     {
409 :     print PROML "13\n"; # Random number seed of unknown use
410 :     }
411 :    
412 :     close PROML;
413 :    
414 :     system "$program < params > /dev/null";
415 :    
416 :     my @likelihoods = &read_outfile();
417 :    
418 :     my @trees = gjonewicklib::read_newick_trees( 'outtree' );
419 :     @trees or print STDERR "proml::proml: Could read proml outtree file\n"
420 :     and chdir $cwd
421 :     and return ();
422 :    
423 :     # We are done, go back to the original directory:
424 :    
425 :     chdir $cwd;
426 :    
427 :     # Returned trees have our labels, and branch lengths that are in % change,
428 :     # not the more usual expected number per position:
429 :    
430 : overbeek 1.4 my @trees = map { gjonewicklib::newick_relabel_nodes( $_, \%id ) } @trees;
431 : overbeek 1.1
432 :     if ( $tree_format =~ m/overbeek/i )
433 :     {
434 :     @trees = map { gjonewicklib::gjonewick_to_overbeek( $_ ) } @trees;
435 :     }
436 :    
437 :     system "/bin/rm -r $tmp_dir" if ! $save_tmp;
438 :    
439 :     return map { [ $_, shift @likelihoods ] } @trees;
440 :     }
441 :    
442 :    
443 :     #-------------------------------------------------------------------------------
444 :     # A perl interface for using proml to estimate site-specific rates of change
445 :     #
446 :     # ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, proml_opts )
447 :     #
448 :     # $categories = [ $n_categories, [ $rates ], $site_categories ];
449 :     #
450 :     # $alignment = [ [ id, def, seq ], ... ]
451 :     # or
452 :     # [ [ id, seq ], ... ]
453 :     #
454 :     # $tree = overbeek tree or gjonewick tree
455 :     #
456 :     # proml_opts is list of key value pairs, or reference to a hash
457 :     #-------------------------------------------------------------------------------
458 :    
459 :     sub estimate_protein_site_rates
460 :     {
461 :     my ( $align, $tree, @proml_opts ) = @_;
462 :    
463 :     my @align = @$align;
464 :    
465 :     # Make the tree into a gjonewick tree so that we can manipulate it:
466 :    
467 :     if ( ref( $tree->[0] ) ne 'ARRAY' ) # overbeek tree
468 :     {
469 :     $tree = gjonewicklib::overbeek_to_gjonewick( $tree );
470 :     }
471 :    
472 :     # The minimum rate will be 1/2 change per total tree branch length.
473 :     # This needs to be checked for proml. The intent is that he optimal
474 :     # rate for a site with one amino acid change is twice this value.
475 :    
476 :     my $kmin = 1 / ( gjonewicklib::newick_tree_length( $tree ) || 1 );
477 :    
478 : overbeek 1.3 # print STDERR "Length = ", gjonewicklib::newick_tree_length( $tree ), "; kmin = $kmin\n"; ## DEBUG ##
479 : overbeek 1.1
480 :     # Generate "rate variation" by rescaling the supplied tree. We could use a
481 :     # finer grain estimator, then categorize the inferred values. This might
482 :     # work slightly better (this is what DNArates currently does).
483 :    
484 :     my $f = exp( log( 2 ) / 1 ); # Interval of 2
485 :     my @rates = map { $kmin * $f**$_ } ( 0 .. 16 ); # kmin .. 65000 * kmin in 17 bins
486 :     my @cat_vals = ( 1 .. 17 );
487 :     my @trees;
488 :     my $rate;
489 :     foreach $rate ( @rates )
490 :     {
491 :     my $tr = gjonewicklib::copy_newick_tree( $tree );
492 :     gjonewicklib::newick_rescale_branches( $tr, $rate ); # Rescales in place
493 :     push @trees, $tr;
494 :     }
495 :    
496 :     # Adjust (a copy of) the proml opts:
497 :    
498 :     my %proml_opts = ( ref( $proml_opts[0] ) eq 'HASH' ) ? %{ $proml_opts[0] } : @proml_opts;
499 :     $proml_opts{ alpha } = undef;
500 :     $proml_opts{ categories } = 0;
501 :     $proml_opts{ coef_of_var } = 0;
502 :     $proml_opts{ gamma_bins } = 0;
503 :     $proml_opts{ invar_frac } = 0;
504 :     $proml_opts{ jumble_seed } = 0;
505 :     $proml_opts{ n_jumble } = 0;
506 :     $proml_opts{ rearrange } = 0;
507 :     $proml_opts{ user_lengths } = 1;
508 :     $proml_opts{ user_trees } = \@trees;
509 :     $proml_opts{ tree_format } = 'gjo';
510 :    
511 :     # Work throught the sites, finding their optimal categories:
512 :    
513 :     my @categories;
514 :     my @weights;
515 :     my $imax = length( $align[0]->[-1] );
516 :     for ( my $i = 0; $i < $imax; $i++ )
517 :     {
518 :     my $inform = 0;
519 :     my @align2 = map { my $c = substr( $_->[-1], $i, 1 );
520 :     $inform ++ if ( $c =~ m/[ACDEFGHIKLMNPQRSTVWY]/i );
521 :     [ $_->[0], $c ]
522 :     }
523 :     @align;
524 :    
525 :     # Only analyze the rate if there are 4 or more informative sequences:
526 :    
527 :     if ( $inform >= 4 )
528 :     {
529 :     my @results = proml::proml( \@align2, \%proml_opts );
530 :    
531 :     my ( $best ) = sort { $b->[1] <=> $a->[1] }
532 :     map { [ $_, @{ shift @results }[1] ] } # get the likelihoods
533 :     @cat_vals;
534 :    
535 : overbeek 1.3 # printf STDERR "%6d %2d => %12.4f\n", $i+1, @$best; ## DEBUG ##
536 : overbeek 1.1 push @categories, $best->[0];
537 :     push @weights, 1;
538 :     }
539 :     else
540 :     {
541 :     push @categories, 9;
542 :     push @weights, 0;
543 :     }
544 :     }
545 :    
546 :     # Find the minimum category value to appear:
547 :    
548 :     my ( $mincat ) = sort { $a <=> $b } @categories;
549 :     my $adjust = $mincat - 1;
550 :    
551 :     @categories = map { min( $_ - $adjust, 9 ) } @categories;
552 :     @rates = @rates[ $adjust .. ( $adjust+8 ) ];
553 :    
554 :     # Return category and weight data:
555 :    
556 :     ( [ scalar @rates, \@rates, join( '', @categories ) ], join( '', @weights ) )
557 :     }
558 :    
559 :    
560 :     sub min { $_[0] < $_[1] ? @_[0] : @_[1] }
561 :    
562 :    
563 :     #-------------------------------------------------------------------------------
564 :     # Auxiliary functions:
565 :     #-------------------------------------------------------------------------------
566 :    
567 :     sub write_infile
568 :     {
569 :     open( INFILE, '>infile' ) or return 0;
570 :     print INFILE scalar @_, ' ', length( $_[0]->[1] ), "\n";
571 :     foreach ( @_ ) { printf INFILE "%-10s %s\n", @$_ }
572 :     close( INFILE );
573 :     }
574 :    
575 :    
576 :     sub write_intree
577 :     {
578 :     open( INTREE, '>intree' ) or return 0;
579 :     print INTREE scalar @_, "\n";
580 :     foreach ( @_ ) { print INTREE gjonewicklib::strNewickTree( $_ ), "\n" }
581 :     close( INTREE );
582 :     }
583 :    
584 :    
585 :     sub write_categories
586 :     {
587 :     my $categories = shift;
588 :     open( CATEGORIES, '>categories' ) or return 0;
589 :     print CATEGORIES "$categories\n";
590 :     close( CATEGORIES );
591 :     }
592 :    
593 :    
594 :     sub write_weights
595 :     {
596 :     my $weights = shift;
597 :     open( WEIGHTS, '>weights' ) or return 0;
598 :     print WEIGHTS "$weights\n";
599 :     close( WEIGHTS );
600 :     }
601 :    
602 :    
603 :     sub read_outfile
604 :     {
605 :     my @likelihoods;
606 :     open( OUTFILE, '<outfile' ) or return ();
607 :     while ( defined( $_ = <OUTFILE> ) )
608 :     {
609 :     next if ! m/^Ln /;
610 :     chomp;
611 :     s/.* //;
612 :     push @likelihoods, $_;
613 :     }
614 :     close( OUTFILE );
615 :     return @likelihoods;
616 :     }
617 :    
618 :    
619 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3