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

Diff of /FigKernelPackages/proml.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2, Sat Jan 13 03:22:44 2007 UTC revision 1.9, Tue Jan 19 01:55:51 2010 UTC
# Line 15  Line 15 
15  #  #
16  #     ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, proml_opts )  #     ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, proml_opts )
17  #  #
18  #     $categories = [ $n_categories, [ $rates ], $site_categories ];  #     $categories = [ [ $rate1, ... ], $site_categories ];
19  #  #
20  #===============================================================================  #===============================================================================
21  #  #
# Line 33  Line 33 
33  #    For proml:  #    For proml:
34  #      alignment    => \@alignment    the way to supply the alignment as an option, rather than first param  #      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)  #      alpha        => float          alpha parameter of gamma distribution (0.5 - inf)
36  #      categories   => [ n, [ rates ], site_categories ]  #      categories   => [ [ rate1, ... ], site_categories ]
37  #      coef_of_var  => float          1/sqrt(alpha) for gamma distribution (D = 0)  #      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)  #      gamma_bins   => int            number of rate categories used to approximate gamma (D=5)
39  #      global       => bool           global rearrangements  #      global       => bool           global rearrangements
# Line 42  Line 42 
42  #      model        => model          evolution model JTT (D) | PMB | PAM  #      model        => model          evolution model JTT (D) | PMB | PAM
43  #      n_jumble     => int            number of jumbles  #      n_jumble     => int            number of jumbles
44  #      persistance  => float          persistance length of rate category  #      persistance  => float          persistance length of rate category
45  #      rate_hmm     => [ n_rates, [ rates ], [ probabilies ] ]  #      rate_hmm     => [ [ rate, prior_prob ] ... ]   # not implimented
46  #      rearrange    => [ trees ]      rearrange user trees  #      rearrange    => [ trees ]      rearrange user trees
47  #      slow         => bool           more accurate but slower search (D = 0)  #      slow         => bool           more accurate but slower search (D = 0)
48  #      user_lengths => bool           use supplied branch lengths  #      user_lengths => bool           use supplied branch lengths
# Line 63  Line 63 
63  #  Options that do not require other data:  #  Options that do not require other data:
64  #    G (global search toggle)  #    G (global search toggle)
65  #    L (user lengths toggle)  #    L (user lengths toggle)
66  #    P (JTT / PAM toggle)  #    P (JTT / PMB / PAM cycle)
67  #    S (slow and accurate)  #    S (slow and accurate)
68  #    U (requires intree file)  #    U (requires intree file)
69  #    W (requires weights file)  #    W (requires weights file)
# Line 94  Line 94 
94  #  Rate values (n of them)  #  Rate values (n of them)
95    
96    
97    use Data::Dumper;
98    
99  use strict;  use strict;
100  use gjonewicklib qw( gjonewick_to_overbeek  use gjonewicklib qw( gjonewick_to_overbeek
101                       newick_is_unrooted                       newick_is_unrooted
102                       newick_relabel_nodes                       newick_relabel_nodes
103                         newick_rescale_branches
104                       newick_tree_length                       newick_tree_length
105                       overbeek_to_gjonewick                       overbeek_to_gjonewick
106                       parse_newick_tree_str                       parse_newick_tree_str
# Line 105  Line 108 
108                       uproot_newick                       uproot_newick
109                     );                     );
110    
   
111  sub proml  sub proml
112  {  {
113      my $align;      my $align;
# Line 148  Line 150 
150      #  Process proml options:      #  Process proml options:
151      #---------------------------------------------------------------------------      #---------------------------------------------------------------------------
152    
153      my $categories   = $options{ categories };  # [ n_cat, [ cat_rates ], site_cats ]      #  [ [ cat_rate1, ... ], site_categories ]
154        #  Original format expected first field to be number of categories (which
155        #  is redundant).  Handling that form is what the shift if all about.
156    
157        my $categories   = $options{ categories };  # [ [ cat_rate1, ... ], site_categories ]
158      if ( $categories )      if ( $categories )
159      {      {
160          if ( ref( $categories ) ne 'ARRAY'          if ( ref( $categories ) ne 'ARRAY'
161            || @$categories != 3            || ! ( ( @$categories == 2 ) || ( ( @$categories == 3 ) && ( shift @$categories ) ) )
162            || $categories->[0] < 1 || $categories->[0] > 9            || ref( $categories->[0] ) ne 'ARRAY'
           || ref( $categories->[1] ) ne 'ARRAY'  
           || @{$categories->[1]} != $categories->[0]  
163             )             )
164          {          {
165              print STDERR "proml::proml categories option value must be [ n_cat, [ cat_rates ], site_cats ]\n";              print STDERR "proml::proml() categories option value must be [ [ cat_rate1, ... ], site_categories ]\n";
166              return ();              return ();
167          }          }
168    
169          #  Rate values cannot have very many decimal places or proml can't read it:          #  Rate values cannot have very many decimal places or proml can't read it:
170    
171          @{$categories->[1]} = map { sprintf "%.6f", $_ } @{$categories->[1]};          @{$categories->[0]} = map { sprintf "%.6f", $_ } @{$categories->[0]};
172      }      }
173    
174      my $coef_of_var  = $options{ coef_of_var }      my $coef_of_var  = $options{ coef_of_var }
# Line 172  Line 176 
176                    ||  0;                    ||  0;
177      if ( $coef_of_var < 0 )      if ( $coef_of_var < 0 )
178      {      {
179          print STDERR "proml::proml coef_of_var option value must be >= 0\n";          print STDERR "proml::proml() coef_of_var option value must be >= 0\n";
180          return ();          return ();
181      }      }
182    
183      my $gamma_bins   = int( $options{ gamma_bins } || ( $coef_of_var ? 5 : 2 ) );      my $gamma_bins   = int( $options{ gamma_bins } || ( $coef_of_var ? 5 : 2 ) );
184      if ( ( $gamma_bins < 2 )  || ( $gamma_bins > 9 ) )      if ( ( $gamma_bins < 2 )  || ( $gamma_bins > 9 ) )
185      {      {
186          print STDERR "proml::proml gamma_bins option value must be > 1 and <= 9\n";          print STDERR "proml::proml() gamma_bins option value must be > 1 and <= 9\n";
187          return ();          return ();
188      }      }
189    
# Line 188  Line 192 
192      my $invar_frac   = $options{ invar_frac } || 0;      my $invar_frac   = $options{ invar_frac } || 0;
193      if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )      if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
194      {      {
195          print STDERR "proml::proml invar_frac option value must be >= 0 and < 1\n";          print STDERR "proml::proml() invar_frac option value must be >= 0 and < 1\n";
196          return ();          return ();
197      }      }
198    
199      my $n_jumble     = int( $options{ n_jumble }    || ( $options{ jumble_seed } ? 1 : 0) );      my $n_jumble     = int( $options{ n_jumble }    || ( $options{ jumble_seed } ? 1 : 0) );
200      if ( $n_jumble < 0 )      if ( $n_jumble < 0 )
201      {      {
202          print STDERR "proml::proml n_jumble option value must be >= 0\n";          print STDERR "proml::proml() n_jumble option value must be >= 0\n";
203          return ();          return ();
204      }      }
205    
206      my $jumble_seed  = int( $options{ jumble_seed } || 4 * int( 499999999 * rand() ) + 1 );      my $jumble_seed  = int( $options{ jumble_seed } || 4 * int( 499999999 * rand() ) + 1 );
207      if ( ( $jumble_seed <= 0)  || ( $jumble_seed % 2 != 1 ) )      if ( ( $jumble_seed <= 0)  || ( $jumble_seed % 2 != 1 ) )
208      {      {
209          print STDERR "proml::proml jumble_seed option value must be an odd number > 0\n";          print STDERR "proml::proml() jumble_seed option value must be an odd number > 0\n";
210          return ();          return ();
211      }      }
212    
# Line 220  Line 224 
224      my $persistance  = $options{ persistance } || 0;      my $persistance  = $options{ persistance } || 0;
225      if ( $persistance && ( $persistance <= 1 ) )      if ( $persistance && ( $persistance <= 1 ) )
226      {      {
227          print STDERR "proml::proml persistance option value must be > 1\n";          print STDERR "proml::proml() persistance option value must be > 1\n";
228          return ();          return ();
229      }      }
230    
# Line 240  Line 244 
244          }          }
245          elsif ( ref( $user_trees->[0] ) ne 'ARRAY' )  # First element not tree          elsif ( ref( $user_trees->[0] ) ne 'ARRAY' )  # First element not tree
246          {          {
247              print STDERR "proml::proml usertree or rearrange option value must be reference to list of trees\n";              print STDERR "proml::proml() user_trees or rearrange option value must be reference to list of trees\n";
248              return ();              return ();
249          }          }
250      }      }
251    
252      my $weights      = $options{ weights };      my $weights      = $options{ weights };
253    
   
254      #---------------------------------------------------------------------------      #---------------------------------------------------------------------------
255      #  Options that are not proml options per se:      #  Options that are not proml options per se:
256      #---------------------------------------------------------------------------      #---------------------------------------------------------------------------
257    
258      my $program     = $options{ program }     || 'proml';      my $program     = $options{ program }     || 'proml';
259    
     my $tmp         = $options{ tmp };  
   
     my $tmp_dir     = $options{ tmp_dir };  
   
260      my $tree_format = $options{ tree_format } =~ m/overbeek/i ? 'overbeek'      my $tree_format = $options{ tree_format } =~ m/overbeek/i ? 'overbeek'
261                      : $options{ tree_format } =~ m/gjo/i      ? 'gjonewick'                      : $options{ tree_format } =~ m/gjo/i      ? 'gjonewick'
262                      : $options{ tree_format } =~ m/fig/i      ? 'fig'                      : $options{ tree_format } =~ m/fig/i      ? 'fig'
263                      :                                           'overbeek'; # Default                      :                                           'overbeek'; # Default
264    
265      my $save_tmp    = $tmp_dir && -d $tmp_dir;      my ( $tmp_dir, $save_tmp ) = temporary_directory( \%options );
     if ( $tmp_dir )  
     {  
         if ( -d $tmp_dir ) { $save_tmp = 1  }  
         else               { mkdir $tmp_dir }  
     }  
     else  
     {  
         $tmp = $tmp && -d  $tmp  ?  $tmp  
              :         -d '/tmp' ? '/tmp'  
              :                     '.';  
         my $int = int( 1000000000 * rand);  
         $tmp_dir = "$tmp/proml.$$.$int";  
         mkdir $tmp_dir;  
     }  
266    
267      #---------------------------------------------------------------------------      #---------------------------------------------------------------------------
268      #  Prepare data:      #  Prepare data:
# Line 293  Line 278 
278      #                            ]      #                            ]
279      #  Root node of gjonewick always has a descendent list.  If the first      #  Root node of gjonewick always has a descendent list.  If the first
280      #  field of the first tree is not an array reference, they are overbeek      #  field of the first tree is not an array reference, they are overbeek
281      #  trees.  Also relabel tree tips to local ids.      #  trees.
282    
283      my @user_trees = ();      my @user_trees = ();
284      if ( $user_trees )      if ( @$user_trees )
285      {      {
286          if ( @user_trees && ( ref( $user_trees[0]->[0] ) ne 'ARRAY' ) )  # overbeek trees          if ( ref( @$user_trees[0]->[0] ) ne 'ARRAY' )  # overbeek trees
287          {          {
288              @user_trees = map { gjonewicklib::newick_relabel_nodes( $_, \%local_id ) }              @user_trees = map { gjonewicklib::overbeek_to_gjonewick( $_ ) }
                           map { gjonewicklib::overbeek_to_gjonewick( $_ ) }  
289                            @$user_trees;                            @$user_trees;
290          }          }
291          else          else
292          {          {
293              @user_trees = map { gjonewicklib::newick_relabel_nodes( $_, \%local_id ) }              @user_trees = map { gjonewicklib::copy_newick_tree( $_ ) }
294                            @$user_trees;                            @$user_trees;
295          }          }
296    
297          # Make sure trees are unrooted:          # Relabel and make sure trees are unrooted:
298    
299          @user_trees = map { gjonewicklib::newick_is_unrooted( $_ ) ? $_          @user_trees = map { gjonewicklib::newick_is_unrooted( $_ ) ? $_
300                                                                     : gjonewicklib::uproot_newick( $_ )                                                                     : gjonewicklib::uproot_newick( $_ )
301                            }                            }
302                          map { gjonewicklib::newick_relabel_nodes( $_, \%local_id ); $_ }
303                        @user_trees;                        @user_trees;
304      }      }
305    
# Line 329  Line 314 
314      unlink 'outfile' if -f 'outfile';  # Just checking      unlink 'outfile' if -f 'outfile';  # Just checking
315      unlink 'outtree' if -f 'outtree';  # ditto      unlink 'outtree' if -f 'outtree';  # ditto
316    
317      &write_infile( @align ) or print STDERR "proml::proml: Could write infile\n"      &write_infile( @align ) or print STDERR "proml::proml: Could not write infile\n"
318                                 and chdir $cwd                                 and chdir $cwd
319                                 and return ();                                 and return ();
320    
321      # open( PROML, "| $program > /dev/null" ) or print STDERR "proml::proml: Could not open pipe to $program\n"      open( PROML, ">params" ) or print STDERR "proml::proml: Could not open command file for $program\n"
     open( PROML, ">params" ) or print STDERR "proml::proml: Could not open pipe to $program\n"  
322                                                 and chdir $cwd                                                 and chdir $cwd
323                                                 and return ();                                                 and return ();
324    
325    
326      #  Start sending optoins ot program:      #  Start writing options for program:
327    
328      if ( $categories )      if ( $categories )
329      {      {
330          &write_categories( $categories->[2] ) or print STDERR "proml::proml: Could not write categories\n"          &write_categories( $categories->[1] ) or print STDERR "proml::proml: Could not write categories\n"
331                                                   and chdir $cwd                                                   and chdir $cwd
332                                                   and return ();                                                   and return ();
333          print PROML "C\n",          print PROML "C\n",
334                      "$categories->[0]\n",                      scalar @{$categories->[0]}, "\n",
335                      join( ' ', @{ $categories->[1] } ), "\n";                      join( ' ', @{ $categories->[0] } ), "\n";
336      }      }
337    
338      if ( $invar_frac || $coef_of_var )      if ( $invar_frac || $coef_of_var )
# Line 366  Line 350 
350      print PROML "P\n"    if $model =~ m/PMB/i;      print PROML "P\n"    if $model =~ m/PMB/i;
351      print PROML "P\nP\n" if $model =~ m/PAM/i;      print PROML "P\nP\n" if $model =~ m/PAM/i;
352    
     print PROML "S\n" if $slow;  
   
353      if ( @user_trees )      if ( @user_trees )
354      {      {
355          &write_intree( @user_trees ) or print STDERR "proml::proml: Could write intree\n"          &write_intree( @user_trees ) or print STDERR "proml::proml: Could not write intree\n"
356                                          and chdir $cwd                                          and chdir $cwd
357                                          and return ();                                          and return ();
358          print PROML "U\n";          print PROML "U\n";
359          print PROML "V\n" if $rearrange || $global;          print PROML "V\n" if $rearrange || $global;
360          print PROML "L\n" if $user_lengths && ! $rearrange && ! $global;          print PROML "L\n" if $user_lengths && ! $rearrange && ! $global;
361      }      }
362        elsif ( $slow )  # Slow and user trees are mutually exclusive
363        {
364            print PROML "S\n";
365        }
366    
367      if ( $weights )      if ( $weights )
368      {      {
369          &write_weights( $weights ) or print STDERR "proml::proml: Could write weights\n"          &write_weights( $weights ) or print STDERR "proml::proml: Could not write weights\n"
370                                        and chdir $cwd                                        and chdir $cwd
371                                        and return ();                                        and return ();
372          print PROML "W\n";          print PROML "W\n";
373      }      }
374    
375      #  All the options are sent, try to lauch the run:      #  All the options are written, try to launch the run:
376    
377      print PROML "Y\n";      print PROML "Y\n";
378    
# Line 404  Line 390 
390          print PROML "$gamma_bins\n";          print PROML "$gamma_bins\n";
391          print PROML "$invar_frac\n" if $invar_frac;          print PROML "$invar_frac\n" if $invar_frac;
392      }      }
393      elsif ( $user_trees )  
394        if ( $user_trees )
395      {      {
396          print PROML "13\n";     #  Random number seed of unknown use          print PROML "13\n";     #  Random number seed of unknown use
397      }      }
# Line 416  Line 403 
403      my @likelihoods = &read_outfile();      my @likelihoods = &read_outfile();
404    
405      my @trees = gjonewicklib::read_newick_trees( 'outtree' );      my @trees = gjonewicklib::read_newick_trees( 'outtree' );
406      @trees or print STDERR "proml::proml: Could read proml outtree file\n"      @trees or print STDERR "proml::proml: Could not read proml outtree file\n"
407                and chdir $cwd                and chdir $cwd
408                and return ();                and return ();
409    
# Line 427  Line 414 
414      #  Returned trees have our labels, and branch lengths that are in % change,      #  Returned trees have our labels, and branch lengths that are in % change,
415      #  not the more usual expected number per position:      #  not the more usual expected number per position:
416    
417      my @trees = map { gjonewicklib::newick_rescale_branches( $_, 0.01 );      @trees = map { gjonewicklib::newick_relabel_nodes( $_, \%id ) } @trees;
                       gjonewicklib::newick_relabel_nodes( $_, \%id )  
                     }  
                 @trees;  
418    
419      if ( $tree_format =~ m/overbeek/i )      if ( $tree_format =~ m/overbeek/i )
420      {      {
# Line 446  Line 430 
430  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
431  #  A perl interface for using proml to estimate site-specific rates of change  #  A perl interface for using proml to estimate site-specific rates of change
432  #  #
433  #     ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, proml_opts )  #     ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree,  %proml_opts )
434    #     ( $categories, $weights ) = estimate_protein_site_rates( \@align, $tree, \%proml_opts )
435  #  #
436  #     $categories = [ $n_categories, [ $rates ], $site_categories ];  #     $categories = [ [ $rate1, ... ], $site_categories ];
437  #  #
438  #  $alignment = [ [ id, def, seq ], ... ]  #  $alignment = [ [ id, def, seq ], ... ]
439  #             or  #             or
# Line 463  Line 448 
448  {  {
449      my ( $align, $tree, @proml_opts ) = @_;      my ( $align, $tree, @proml_opts ) = @_;
450    
451      my @align = @$align;      my ( $seq, $id );
452        my %local_id;
453        my $local_id = 'seq0000000';
454        my @align = map { $id = $_->[0];
455                          $local_id{ $id } = ++$local_id;
456                          $seq = $_->[-1];
457                          $seq =~ s/[BJOUZ]/X/gi;  # Bad letters go to X
458                          $seq =~ s/[^A-Z]/-/gi;   # Anything else becomes -
459                          [ $local_id, $seq ]
460                        } @$align;
461    
462      #  Make the tree into a gjonewick tree so that we can manipulate it:      #  Make the tree a gjonewick tree, uproot it, and change to the local ids.
463    
464      if ( ref( $tree->[0] ) ne 'ARRAY' )   # overbeek tree      if ( ref( $tree->[0] ) ne 'ARRAY' )   # overbeek tree
465      {      {
466          $tree = gjonewicklib::overbeek_to_gjonewick( $tree );          $tree = gjonewicklib::overbeek_to_gjonewick( $tree );
467      }      }
468        else
469        {
470            $tree = gjonewicklib::copy_newick_tree( $tree );
471        }
472    
473        $tree = gjonewicklib::uproot_newick( $tree ) if ! gjonewicklib::newick_is_unrooted( $tree );
474    
475        gjonewicklib::newick_relabel_nodes( $tree, \%local_id );
476    
477      #  The minimum rate will be 1/2 change per total tree branch length.      #  The minimum rate will be 1/2 change per total tree branch length.
478      #  This needs to be checked for proml.  The intent is that he optimal      #  This needs to be checked for proml.  The intent is that he optimal
# Line 478  Line 480 
480    
481      my $kmin = 1 / ( gjonewicklib::newick_tree_length( $tree ) || 1 );      my $kmin = 1 / ( gjonewicklib::newick_tree_length( $tree ) || 1 );
482    
     print STDERR "Length = ", gjonewicklib::newick_tree_length( $tree ), "; kmin = $kmin\n"; ## DEBUG ##  
   
483      #  Generate "rate variation" by rescaling the supplied tree.  We could use a      #  Generate "rate variation" by rescaling the supplied tree.  We could use a
484      #  finer grain estimator, then categorize the inferred values.  This might      #  finer grain estimator, then categorize the inferred values.  This might
485      #  work slightly better (this is what DNArates currently does).      #  work slightly better (this is what DNArates currently does).
# Line 499  Line 499 
499      #  Adjust (a copy of) the proml opts:      #  Adjust (a copy of) the proml opts:
500    
501      my %proml_opts = ( ref( $proml_opts[0] ) eq 'HASH' ) ? %{ $proml_opts[0] } : @proml_opts;      my %proml_opts = ( ref( $proml_opts[0] ) eq 'HASH' ) ? %{ $proml_opts[0] } : @proml_opts;
502      $proml_opts{ alpha        } =  undef;  
     $proml_opts{ categories   } =  0;  
     $proml_opts{ coef_of_var  } =  0;  
     $proml_opts{ gamma_bins   } =  0;  
     $proml_opts{ invar_frac   } =  0;  
     $proml_opts{ jumble_seed  } =  0;  
     $proml_opts{ n_jumble     } =  0;  
     $proml_opts{ rearrange    } =  0;  
503      $proml_opts{ user_lengths } =  1;      $proml_opts{ user_lengths } =  1;
504      $proml_opts{ user_trees   } = \@trees;      $proml_opts{ user_trees   } = \@trees;
505      $proml_opts{ tree_format  } = 'gjo';      $proml_opts{ tree_format  } = 'gjo';
506    
507      #  Work throught the sites, finding their optimal categories:      delete $proml_opts{ alpha       } if exists $proml_opts{ alpha       };
508        delete $proml_opts{ categories  } if exists $proml_opts{ categories  };
509        delete $proml_opts{ coef_of_var } if exists $proml_opts{ coef_of_var };
510        delete $proml_opts{ gamma_bins  } if exists $proml_opts{ gamma_bins  };
511        delete $proml_opts{ invar_frac  } if exists $proml_opts{ invar_frac  };
512        delete $proml_opts{ jumble_seed } if exists $proml_opts{ jumble_seed };
513        delete $proml_opts{ n_jumble    } if exists $proml_opts{ n_jumble    };
514        delete $proml_opts{ rearrange   } if exists $proml_opts{ rearrange   };
515    
516        #  Work throught the sites, finding their optimal rates/categories:
517    
518      my @categories;      my @categories;
519      my @weights;      my @weights;
# Line 535  Line 537 
537                             map  { [ $_, @{ shift @results }[1] ] }  # get the likelihoods                             map  { [ $_, @{ shift @results }[1] ] }  # get the likelihoods
538                             @cat_vals;                             @cat_vals;
539    
540              printf STDERR "%6d  %2d => %12.4f\n", $i+1, @$best; ## DEBUG ##  #           printf STDERR "%6d  %2d => %12.4f\n", $i+1, @$best; ## DEBUG ##
541              push @categories, $best->[0];              push @categories, $best->[0];
542              push @weights,    1;              push @weights,    1;
543          }          }
# Line 556  Line 558 
558    
559      #  Return category and weight data:      #  Return category and weight data:
560    
561      ( [ scalar @rates, \@rates, join( '', @categories ) ], join( '', @weights ) )      ( [ \@rates, join( '', @categories ) ], join( '', @weights ) )
562  }  }
563    
564    
 sub min { $_[0] < $_[1] ? @_[0] : @_[1] }  
   
   
565  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
566  #  Auxiliary functions:  #  Auxiliary functions:
567  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
568    
569    sub min { $_[0] < $_[1] ? $_[0] : $_[1] }
570    
571    
572  sub write_infile  sub write_infile
573  {  {
574      open( INFILE, '>infile' ) or return 0;      open( INFILE, '>infile' ) or return 0;
# Line 605  Line 607 
607    
608  sub read_outfile  sub read_outfile
609  {  {
     my @likelihoods;  
610      open( OUTFILE, '<outfile' ) or return ();      open( OUTFILE, '<outfile' ) or return ();
611      while ( defined( $_ = <OUTFILE> ) )      my @likelihoods = map  { chomp; s/.* //; $_ }
612      {                        grep { /^Ln Likelihood/ }
613          next if ! m/^Ln /;                        <OUTFILE>;
         chomp;  
         s/.* //;  
         push @likelihoods, $_;  
     }  
614      close( OUTFILE );      close( OUTFILE );
615      return @likelihoods;      return @likelihoods;
616  }  }
617    
618    
619    #-------------------------------------------------------------------------------
620    #  The SEED has things in special places.  Be aware of them if running SEED.
621    #-------------------------------------------------------------------------------
622    
623    my $tmp;        #  The default place for temp files or directories
624    my $ext_bin;    #  FIG path to external binaries
625    
626    eval { require FIG_Config;
627           $tmp     = $FIG_Config::temp;
628           $ext_bin = $FIG_Config::ext_bin;
629         };
630    
631    
632    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
633    #  $program = executable( $program, \%options )
634    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
635    sub executable
636    {
637        my ( $program, $options ) = @_;
638        return undef if ! $program;
639    
640        $options ||= {};
641    
642        return $options->{ $program } ? $options->{ $program } # explicit?
643             : $ext_bin               ? "$ext_bin/$program"    # SEED?
644             :                          $program;
645    }
646    
647    
648    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
649    #  Find or create a temporary directory.
650    #  If it does not exist, create it.  If it exists, mark it for saving.
651    #
652    #    $tmp_dir              = temporary_directory( \%options )
653    #  ( $tmp_dir, $save_tmp ) = temporary_directory( \%options )
654    #
655    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
656    
657    sub temporary_directory
658    {
659        my $options = shift || {};
660    
661        my $tmp_dir  = $options->{ tmpdir }  || $options->{ tmp_dir };
662        my $save_tmp = $options->{ savetmp } || $options->{ save_tmp } || '';
663    
664        if ( $tmp_dir )
665        {
666            if ( -d $tmp_dir ) { $options->{ savetmp } = $save_tmp = 1 }
667        }
668        else
669        {
670            if ( $options->{ tmp }  && -d  $options->{ tmp } )
671            {
672                $tmp = $options->{ tmp };
673            }
674            elsif ( ! $tmp || ! -d $tmp )
675            {
676                $options->{ tmp } = $tmp = -d '/tmp' ? '/tmp' : '.';
677            }
678    
679            $tmp_dir = sprintf( "$tmp/" . __PACKAGE__ . "_tmp.%05d.%09d", $$, int(1000000000*rand) );
680            $options->{ tmpdir } = $tmp_dir;
681        }
682    
683        if ( $tmp_dir && ! -d $tmp_dir )
684        {
685            mkdir $tmp_dir;
686            if ( ! -d $tmp_dir )
687            {
688                print STDERR __PACKAGE__ . "::temporary_directory could not create '$tmp_dir'\n";
689                $options->{ tmpdir } = $tmp_dir = undef;
690            }
691        }
692    
693        return wantarray ? ( $tmp_dir, $save_tmp ) : $tmp_dir;
694    }
695    
696    
697  1;  1;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3