File Coverage

blib/lib/Bio/Tools/Run/Phylo/PAML/Codeml.pm

Criterion	Covered	Total	%
statement	122	179	68.1
branch	37	94	39.3
condition	6	32	18.7
subroutine	23	25	92.0
pod	13	13	100.0
total	201	343	58.6

package Bio::Tools::Run::Phylo::PAML::Codeml;
$Bio::Tools::Run::Phylo::PAML::Codeml::VERSION = '1.7.2';
use utf8;
use strict;
use warnings;

use vars qw(@ISA %VALIDVALUES $MINNAMELEN $PROGRAMNAME $PROGRAM);
use Bio::Root::Root;
use Bio::AlignIO;
use Bio::TreeIO;
use Bio::Tools::Run::WrapperBase;
use Bio::Tools::Phylo::PAML;
use Cwd;

@ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase);

# ABSTRACT: Wrapper aroud the PAML program codeml
# AUTHOR: Jason Stajich
# OWNER: Jason Stajich
# LICENSE: Perl_5



BEGIN {

    $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],
             );
}


sub program_name {
        return 'codeml';
}


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



sub new {
  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;
}



sub prepare{
   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;
}




sub run {
   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);
}


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

}


sub alignment{
   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'};
}


sub tree {
   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'};
}


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



sub set_parameter{
   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;
}


sub set_default_parameters{
   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;
       }
   }
}




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




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




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

1;

__END__