use Bio::Tools::Run::Phylo::PAML::Codeml;
  use Bio::AlignIO;

  my $alignio = Bio::AlignIO->new(-format => 'phylip',
  			         -file   => 't/data/gf-s85.phylip');

  my $aln = $alignio->next_aln;

  my $codeml = Bio::Tools::Run::Phylo::PAML::Codeml->new();
  $codeml->alignment($aln);
  my ($rc,$parser) = $codeml->run();
  my $result = $parser->next_result;
  my $MLmatrix = $result->get_MLmatrix();
  print "Ka = ", $MLmatrix->[0]->[1]->{'dN'},"\n";
  print "Ks = ", $MLmatrix->[0]->[1]->{'dS'},"\n";
  print "Ka/Ks = ", $MLmatrix->[0]->[1]->{'omega'},"\n";

 
    $MINNAMELEN = 25;
    $PROGRAMNAME = 'codeml' . ($^O =~ /mswin/i ?'.exe':'');
    if( defined $ENV{'PAMLDIR'} ) {
	$PROGRAM = Bio::Root::IO->catfile($ENV{'PAMLDIR'},$PROGRAMNAME). ($^O =~ /mswin/i ?'.exe':'');;
    }
   
    # valid values for parameters, the default one is always
    # the first one in the array
    # much of the documentation here is lifted directly from the codeml.ctl
    # example file provided with the package
    %VALIDVALUES = ( 
		     'outfile' => 'mlc',
		     'noisy'   => [ 0..3,9],
		     'verbose' => [ 1,0,2], # 0:concise, 1:detailed, 2:too much

                     # (runmode) 0:user tree, 1:semi-autmatic, 2:automatic
		     #           3:stepwise addition, 4,5:PerturbationNNI
		     #           -2:pairwise
		     'runmode' => [ -2, 0..5], 

		     'seqtype' => [ 1..3], # 1:codons, 2:AAs, 3:codons->AAs

		     'CodonFreq' => [ 2, 0,1,3,4,5,6,7], # 0:1/61 each, 1:F1X4, 
		                                # 2:F3X4, 3:codon table

		     # (aaDist) 0:equal, +:geometric, -:linear, 
		     #          1-6:G1974,Miyata, c,p,v,a
		     'aaDist'  => [ 0,'+','-', 1..6], 

                     # (aaRatefile) only used for aa seqs 
		     # with model=empirical(_F)
		     # default is usually 'wag.dat', also
		     # dayhoff.dat, jones.dat, mtmam.dat, or your own
		     'aaRatefile' => 'wag.dat', 

		     # (model) models for codons 
		     # 0: one, 1:b, 2:2 or more dN/dS ratios for branches
		     'model'    => [0..3,7], 
		     
		     # (NSsites) number of S sites
		     # 0: one w;1:neutral;2:selection; 3:discrete;4:freqs;
                     # 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
                     # 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
                     # 13:3normal>0
		     'NSsites'  => [0..13],

		     # (icode) genetic code
		     # 0:universal code
		     # 1:mamalian mt
		     # 2:yeast mt
		     # 3:mold mt,
		     # 4:invertebrate mt
		     # 5:ciliate nuclear
		     # 6:echinoderm mt
		     # 7:euplotid mt
		     # 8:alternative yeast nu.
		     # 9:ascidian mt
		     #10:blepharisma nu
		     # these correspond to 1-11 in the genbank transl table
		     
		     'icode'    => [ 0..10], 
		     
		     'Mgene'    => [0,1], # 0:rates, 1:separate
		     
		     'fix_kappa'=> [0,1], # 0:estimate kappa, 1:fix kappa
		     'kappa'    => '2',   # initial or fixed kappa
		     'fix_omega'=> [0,1], # 0: estimate omega, 1: fix omega
		     'omega'    => '1', # initial or fixed omega for 
		                          # codons or codon-base AAs
		     'fix_alpha'=> [1,0], # 0: estimate gamma shape param
		                          # 1: fix it at alpha
		     'alpha'    => '0.', # initial or fixed alpha
		                        # 0: infinity (constant rate)
		     'Malpha'   => [0,1], # different alphas for genes
		     'ncatG'    => [1..10], # number of categories in 
		                        # dG of NSsites models

		     # (clock)
		     # 0: no clock, 1: global clock, 2: local clock
		     # 3: TipDate
		     'clock'    => [0..3],
		     # (getSE) Standard Error: 
		     # 0:don't want them, 1: want S.E.
		     'getSE'    => [0,1], 
		     # (RateAncestor) 
		     # 0,1,2 rates (alpha>0) or
		     # ancestral states (1 or 2)
		     'RateAncestor' => [1,0,2], 
		     'Small_Diff'    => '.5e-6',
                     # (cleandata) remove sites with ambiguity data
		     # 1: yes, 0:no 
		     'cleandata'     => [0,1], 
		     # this is the number of datasets in
		     # the file - we would need to change
		     # our api to allow >1 alignment object
		     # to be referenced at time
		     'ndata'         => 1,
		     # (method) 
		     # 0: simultaneous,1: 1 branch at a time
		     'method'        => [0,1], 

		     # allow branch lengths to be fixed
		     # 0 ignore
		     # -1 use random starting points
		     # 1 use the branch lengths in initial ML iteration
		     # 2 branch lengths are fixed
		     'fix_blength'   => [0,-1,1,2],
		     );

        return 'codeml';

        return Bio::Root::IO->catfile($ENV{PAMLDIR}) if $ENV{PAMLDIR};

  my($class,@args) = @_;

  my $self = $class->SUPER::new(@args);
  $self->{'_branchLengths'} = 0;
  my ($aln, $tree, $st, $params, $exe, 
      $ubl) = $self->_rearrange([qw(ALIGNMENT TREE SAVE_TEMPFILES 
				    PARAMS EXECUTABLE BRANCHLENGTHS)],
				    @args);
  defined $aln && $self->alignment($aln);
  defined $tree && $self->tree($tree, branchLengths => ($ubl || 0) );
  defined $st  && $self->save_tempfiles($st);
  defined $exe && $self->executable($exe);

  $self->set_default_parameters();
  if( defined $params ) {
      if( ref($params) !~ /HASH/i ) { 
	  $self->warn("Must provide a valid hash ref for parameter -FLAGS");
      } else {
	  map { $self->set_parameter($_, $$params{$_}) } keys %$params;
      }
  }
  return $self;

   my ($self,$aln,$tree) = @_;
   unless ( $self->save_tempfiles ) {
       # brush so we don't get plaque buildup ;)
       $self->cleanup();
   }
   $tree = $self->tree unless $tree;
   $aln  = $self->alignment unless $aln;
   if( ! $aln ) { 
       $self->warn("must have supplied a valid alignment file in order to run codeml");
       return 0;
   }
   my ($tempdir) = $self->tempdir();
   my ($tempseqFH,$tempseqfile);
   if( ! ref($aln) && -e $aln ) { 
       $tempseqfile = $aln;
   } else { 
       ($tempseqFH,$tempseqfile) = $self->io->tempfile
	   ('-dir' => $tempdir, 
	    UNLINK => ($self->save_tempfiles ? 0 : 1));
       my $alnout = Bio::AlignIO->new('-format'      => 'phylip',
				     '-fh'          => $tempseqFH,
                                     '-interleaved' => 0,
                                     '-idlength'    => $MINNAMELEN > $aln->maxdisplayname_length() ? $MINNAMELEN : $aln->maxdisplayname_length() +1);
       
       $alnout->write_aln($aln);
       $alnout->close();
       undef $alnout;   
       close($tempseqFH);
   }
   # now let's print the codeml.ctl file.
   # many of the these programs are finicky about what the filename is 
   # and won't even run without the properly named file.  Ack
   
   my $codeml_ctl = "$tempdir/codeml.ctl";
   open(CODEML, ">$codeml_ctl") or $self->throw("cannot open $codeml_ctl for writing");
   print CODEML "seqfile = $tempseqfile\n";
   my $outfile = $self->outfile_name;
   print CODEML "outfile = $outfile\n";

   if( $tree ) {
       my ($temptreeFH,$temptreefile);
       if( ! ref($tree) && -e $tree ) { 
	   $temptreefile = $tree;
       } else { 
	   ($temptreeFH,$temptreefile) = $self->io->tempfile
	       ('-dir' => $tempdir, 
		UNLINK => ($self->save_tempfiles ? 0 : 1));
	   
	   my $treeout = Bio::TreeIO->new('-format' => 'newick',
					  '-fh'     => $temptreeFH);
	   $treeout->write_tree($tree);
	   $treeout->close();
	   close($temptreeFH);
       }
       print CODEML "treefile = $temptreefile\n";
   }

   my %params = $self->get_parameters;
   while( my ($param,$val) = each %params ) {
       next if $param eq 'outfile';
       print CODEML "$param = $val\n";
   }
   close(CODEML);
#    my ($rc,$parser) = (1);
#    {
#        my $cwd = cwd();
#        my $exit_status;
#        chdir($tempdir);
#    }
   return $tempdir;

   my ($self) = shift;;
   my $outfile = $self->outfile_name;
   my $tmpdir = $self->prepare(@_);
   
   my ($rc,$parser) = (1);
   {
       my $cwd = cwd();
       my $exit_status;
       chdir($tmpdir);
       my $codemlexe = $self->executable();
       $self->throw("unable to find or run executable for 'codeml'") unless $codemlexe && -e $codemlexe && -x _;
       my $run;
       if( $self->{'_branchLengths'} ) { 
	   open($run, "echo $self->{'_branchLengths'} | $codemlexe |") or $self->throw("Cannot open exe $codemlexe");
       } else {
	   open($run, "$codemlexe |") or $self->throw("Cannot open exe $codemlexe");
       }
       my @output = <$run>;
       $exit_status = close($run);
       $self->error_string(join('',@output));
       if( (grep { /\berr(or)?: /io } @output)  || !$exit_status) {
	   $self->warn("There was an error - see error_string for the program output");
	   $rc = 0;
       }
       eval {
	   $parser = Bio::Tools::Phylo::PAML->new(-file => "$tmpdir/$outfile",
						  -verbose => $self->verbose,
						  -dir => "$tmpdir");
       };
       if( $@ ) {
	   $self->warn($self->error_string);
       }
       chdir($cwd);
   }
   return ($rc,$parser);

   my ($self,$value) = @_;
   if( defined $value) {
      $self->{'error_string'} = $value;
    }
    return $self->{'error_string'};

   my ($self,$aln) = @_;

   if( defined $aln ) { 
       if( -e $aln ) { 
	   $self->{'_alignment'} = $aln;
       } elsif( !ref($aln) || ! $aln->isa('Bio::Align::AlignI') ) { 
	   $self->warn("Must specify a valid Bio::Align::AlignI object to the alignment function not $aln");
	   return undef;
       } else {
	   $self->{'_alignment'} = $aln;
       }
   }
   return  $self->{'_alignment'};

   my ($self, $tree, %params) = @_;
   if( defined $tree ) { 
       if( ! ref($tree) || ! $tree->isa('Bio::Tree::TreeI') ) { 
	   $self->warn("Must specify a valid Bio::Tree::TreeI object to the alignment function");
       }
       $self->{'_tree'} = $tree;
       if ( defined $params{'_branchLengths'} ) {
	 my $ubl = $params{'_branchLengths'};
	 if ($ubl !~ m/^(0|1|2)$/) {
	   $self->throw("The branchLengths parameter to tree() must be 0 (ignore), 1 (initial values) or 2 (fixed values) only");
	 }
	 $self->{'_branchLengths'} = $ubl;
       }
   }
   return $self->{'_tree'};

   my ($self) = @_;
   # we're returning a copy of this
   return %{ $self->{'_codemlparams'} };

   my ($self,$param,$value) = @_;
   unless (defined $self->{'no_param_checks'} && $self->{'no_param_checks'} == 1) {
       if ( ! defined $VALIDVALUES{$param} ) { 
           $self->warn("unknown parameter $param will not be set unless you force by setting no_param_checks to true");
           return 0;
       } 
       if ( ref( $VALIDVALUES{$param}) =~ /ARRAY/i &&
            scalar @{$VALIDVALUES{$param}} > 0 ) {
       
           unless ( grep { $value eq $_ } @{ $VALIDVALUES{$param} } ) {
               $self->warn("parameter $param specified value $value is not recognized, please see the documentation and the code for this module or set the no_param_checks to a true value");
               return 0;
           }
       }
   }
   $self->{'_codemlparams'}->{$param} = $value;
   return 1;

   my ($self,$keepold) = @_;
   $keepold = 0 unless defined $keepold;
   
   while( my ($param,$val) = each %VALIDVALUES ) {
       # skip if we want to keep old values and it is already set
       next if( defined $self->{'_codemlparams'}->{$param} && $keepold);
       if(ref($val)=~/ARRAY/i ) {
	   $self->{'_codemlparams'}->{$param} = $val->[0];
       }  else { 
	   $self->{'_codemlparams'}->{$param} = $val;
       }
   }

   my ($self,$value) = @_;
   if( defined $value) {
      $self->{'no_param_checks'} = $value;
    }
    return $self->{'no_param_checks'};

    my $self = shift;
    if( @_ ) {
	return $self->{'_codemlparams'}->{'outfile'} = shift @_;
    }
    unless (defined $self->{'_codemlparams'}->{'outfile'}) {
        $self->{'_codemlparams'}->{'outfile'} = 'mlc';
    }
    return $self->{'_codemlparams'}->{'outfile'};

    my $self= shift;
    unless ( $self->save_tempfiles ) {
	$self->cleanup();
    }
    $self->SUPER::DESTROY();
}

1;

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  http://redmine.open-bio.org/projects/bioperl/

 Title   : save_tempfiles
 Usage   : $obj->save_tempfiles($newval)
 Function: 
 Returns : value of save_tempfiles
 Args    : newvalue (optional)

 Title   : tempdir
 Usage   : my $tmpdir = $self->tempdir();
 Function: Retrieve a temporary directory name (which is created)
 Returns : string which is the name of the temporary directory
 Args    : none

 Title   : cleanup
 Usage   : $codeml->cleanup();
 Function: Will cleanup the tempdir directory after a PAML run
 Returns : none
 Args    : none

 Title   : io
 Usage   : $obj->io($newval)
 Function:  Gets a Bio::Root::IO object
 Returns : Bio::Root::IO
 Args    : none