Bio::EnsEMBL::Pipeline::GeneDuplication PAML
Included librariesPackage variablesGeneral documentationMethods
Toolbar
WebCvsRaw content
Package variables
No package variables defined.
Included modules
Bio::EnsEMBL::Pipeline::RunnableI
Bio::EnsEMBL::Utils::Argument qw ( rearrange )
Bio::EnsEMBL::Utils::Exception qw ( warning throw )
Bio::Tools::Phylo::PAML
Cwd
Inherit
Bio::EnsEMBL::Pipeline::RunnableI
Synopsis
No synopsis!
Description
No description!
Methods
DESTROY
No description
Code
_aadist
No description
Code
_aaratefile
No description
Code
_aligned_seqs
No description
Code
_alpha
No description
Code
_clock
No description
Code
_codeml_executable
No description
Code
_codonfreq
No description
Code
_config_file
No description
Code
_fix_alpha
No description
Code
_fix_kappa
No description
Code
_fix_omega
No description
Code
_getse
No description
Code
_icode
No description
Code
_kappa
No description
Code
_malpha
No description
Code
_mgene
No description
Code
_model
No description
Code
_ncatg
No description
Code
_noisy
No description
Code
_nssites
No description
Code
_omega
No description
Code
_outfile
No description
Code
_pre_copying
No description
Code
_rateancestor
No description
Code
_runmode
No description
Code
_seqfile
No description
Code
_seqtype
No description
Code
_small_diff
No description
Code
_treefile
No description
Code
_verbose
No description
Code
_write_config_file
No description
Code
_write_seqs
No description
Code
new
No description
Code
run_codeml
No description
Code
work_dir
No description
Code
Methods description
None available.
Methods code
DESTROYdescriptionprevnextTop
sub DESTROY {
  my $self = shift;

  unlink $self->_config_file if $self->_config_file; 
  unlink $self->_seqfile if $self->_seqfile;

  unlink ($self->work_dir . "/rub",    $self->work_dir . "/rst",
	  $self->work_dir . "/rst1",   $self->work_dir . "/lnf",
	  $self->work_dir . "/2ML.dN", $self->work_dir . "/2ML.dS",
	  $self->work_dir . "/2ML.t",  $self->work_dir . "/2NG.dN", 
	  $self->work_dir . "/2NG.dS", $self->work_dir . "/2NG.t") 
    if $self->_outfile;

  unlink $self->_outfile if $self->_outfile;

  warning ("Paml execution resulted in a core file being dumped.")
    if (unlink $self->work_dir . "/core");

  warning ("Could not remove working directory [". $self->work_dir ."].\n$!") 
    if (! rmdir $self->work_dir);
}
_aadistdescriptionprevnextTop
sub _aadist {
  my ($self) = shift;

  if (@_) {
    $self->{_aadist} = shift;
  }

  unless (defined $self->{_aadist}){
#    warning ("Aadist not set.  Defaulting to option 1: equal distances.");
return 0; } return $self->{_aadist}
}
_aaratefiledescriptionprevnextTop
sub _aaratefile {
  my ($self) = shift;

  if (@_) {
    $self->{_aaratefile} = shift;
  }

  if (!defined $self->{_aaratefile}){
    if ($self->_seqtype == 2 && $self->_model > 1) {
#      warning ("Aaratefile not set.  Defaulting to file wag.dat.");
} return 'wag.dat'; } return $self->{_aaratefile}
}
_aligned_seqsdescriptionprevnextTop
sub _aligned_seqs {
  my $self = shift;

  if (@_){
    $self->{_aligned_seqs} = shift;

    throw ("No sequences to set.")
      unless $self->{_aligned_seqs} and 
	$self->{_aligned_seqs}->[0]->isa("Bio::Seq");
  }

  return $self->{_aligned_seqs}
}
_alphadescriptionprevnextTop
sub _alpha {
  my ($self) = shift;

  if (@_) {
    $self->{_alpha} = shift;
  }

  if ($self->_fix_alpha != 0 && !defined $self->{_alpha}){
    throw ("Youve asked for alpha to be fixed, but\n" .
		 "have not specified a value for it to be\n" .
		 "fixed to.  Try setting '-alpha' => x or\n" . 
		 "whatever");
  } 

  if (! defined $self->{_alpha}){
    return 0;
  }

  return $self->{_alpha}
}
_clockdescriptionprevnextTop
sub _clock {
  my ($self) = shift;

  if (@_) {
    $self->{_clock} = shift;
  }

  if (!defined $self->{_clock}){
    # What other parameters are set with clock - should check these before throwing a warning.
# warning ("Clock parameter not set. Defaulting to option 0: no clock.");
return 0 } return $self->{_clock}
}
_codeml_executabledescriptionprevnextTop
sub _codeml_executable {
  my $self = shift;

  if (@_) {
    $self->{_codeml_executable} = shift;
  }

  unless ($self->{_codeml_executable}) {
    print STDERR "Explicit location of codeml executable not set.  Assuming\n" .
      "that codeml is in your path.";

    $self->{_codeml_executable} = 'codeml';
  }

  return $self->{_codeml_executable}
}
_codonfreqdescriptionprevnextTop
sub _codonfreq {
  my ($self) = shift;

  if (@_) {
    $self->{_codonfreq} = shift;
  }

  if ($self->_seqtype == 2 && !defined $self->{_codonfreq}){
#    warning ("Codonfreq not set.  Defaulting to option 2: F3X4 frequencies.");
return 2; } if (! $self->{_codonfreq}){ return 2; } return $self->{_codonfreq}
}
_config_filedescriptionprevnextTop
sub _config_file {
  my ($self) = shift;


  if (@_){
    $self->{_config_file} = shift;
  }

#print "Config file for PAML : " . $self->{_config_file} . "\n";
return $self->{_config_file} } ### Config options ###
}
_fix_alphadescriptionprevnextTop
sub _fix_alpha {
  my ($self) = shift;

  if (@_) {
    $self->{_fix_alpha} = shift;
  }

  unless (defined $self->{_fix_alpha}){
#    print STDERR "Defaulting to unfixed alpha - will estimate alpha.\n";
return 0 # default to unfixed alpha
} return $self->{_fix_alpha}
}
_fix_kappadescriptionprevnextTop
sub _fix_kappa {
  my ($self) = shift;

  if (@_) {
    $self->{_fix_kappa} = shift;
  }

  unless (defined $self->{_fix_kappa}){
#    print STDERR "Defaulting to unfixed kappa - will estimate kappa.\n";
return 0 # default to unfixed kappa
} return $self->{_fix_kappa}
}
_fix_omegadescriptionprevnextTop
sub _fix_omega {
  my ($self) = shift;

  if (@_) {
    $self->{_fix_omega} = shift;
  }

  unless (defined $self->{_fix_omega}){
#    print STDERR "Defaulting to unfixed omega - will estimate omega.\n";
return 0 # default to unfixed omega
} return $self->{_fix_omega}
}
_getsedescriptionprevnextTop
sub _getse {
  my ($self) = shift;

  if (@_) {
    $self->{_getse} = shift;
  }

  if (!defined $self->{_getse}){
#    print STDERR "By default, not estimating SE.\n";
return 0 } return $self->{_getse}
}
_icodedescriptionprevnextTop
sub _icode {
  my ($self) = shift;

  if (@_) {
    $self->{_icode} = shift;
  }

  unless (defined $self->{_icode}){
    warning ("PAML: Defaulting to Universal Genetic Code.\n");
    return 0 # Universal code
} return $self->{_icode}
}
_kappadescriptionprevnextTop
sub _kappa {
  my ($self) = shift;

  if (@_) {
    $self->{_kappa} = shift;
  }

  if ($self->_fix_kappa != 0 && !defined $self->{_kappa}){
    throw ("Youve asked for kappa to be fixed, but\n" .
		 "have not specified a value for it to be\n" .
		 "fixed to.  Try setting '-kappa' => x or\n" . 
		 "whatever");
  } 

  if (! $self->{_kappa}){
    return 2;
  }

  return $self->{_kappa}
}
_malphadescriptionprevnextTop
sub _malpha {
  my ($self) = shift;

  if (@_) {
    $self->{_malpha} = shift;
  }

  if (!defined $self->{_malpha}){
    # silently default to 0.  Ive never used this - what 
# other parameters does it depend on?
return 0 } return $self->{_malpha}
}
_mgenedescriptionprevnextTop
sub _mgene {
  my ($self) = shift;

  if (@_) {
    $self->{_mgene} = shift;
  }

  unless (defined $self->{_mgene}){
#    print STDERR "Defaulting to mgene = 0\n";
return 0 # default to separate
} return $self->{_mgene}
}
_modeldescriptionprevnextTop
sub _model {
  my ($self, $set, $value) = @_;

  if ($set) {
    $self->{_model} = $value + 1; # Add one to cope with 0 as a valid option.
} unless (defined $self->{_model}){ throw ("Model has not been set."); } return $self->{_model} - 1 # Remove one to obtain correct value.
}
_ncatgdescriptionprevnextTop
sub _ncatg {
  my ($self) = shift;

  if (@_) {
    $self->{_ncatg} = shift;
  }

  if (!defined $self->{_ncatg}){
    # silently default to 10.  Ive never used this - what 
# other parameters does it depend on?
return 10 } return $self->{_ncatg}
}
_noisydescriptionprevnextTop
sub _noisy {
  my ($self) = shift;

  if (@_) {
    $self->{_noisy} = shift;
  }

  if (!$self->{_noisy}){
    return 0
  } else {
    return $self->{_noisy}
  }
}
_nssitesdescriptionprevnextTop
sub _nssites {
  my ($self) = shift;

  if (@_) {
    $self->{_nssites} = shift;
  }

  unless (defined $self->{_nssites}){
#    warning ("NSsites has not been set.  Defaulting to 0.");
} return $self->{_nssites}
}
_omegadescriptionprevnextTop
sub _omega {
  my ($self) = shift;

  if (@_) {
    $self->{_omega} = shift;
  }

  if ($self->_fix_omega != 0 && !defined $self->{_omega}){
    throw ("Youve asked for omega to be fixed, but\n" .
		 "have not specified a value for it to be\n" .
		 "fixed to.  Try setting '-omega' => x or\n" . 
		 "whatever");
  } 

  if (! defined $self->{_omega}){
    return 0;
  }

  return $self->{_omega}
}
_outfiledescriptionprevnextTop
sub _outfile {
  my ($self) = shift;

  if (@_) {
    $self->{_outfile} = shift;
  }

  unless ($self->{_outfile}) {
    $self->{_outfile} = $self->work_dir . 'paml_out_' . time;
  }

  return $self->{_outfile}
}
_pre_copyingdescriptionprevnextTop
sub _pre_copying {
  my $self = shift;

  #system("cp " . $self->_executable . " " . $self->work_dir);
# Ultimately all the datafiles need to be copied into here.
return 1;
}
_rateancestordescriptionprevnextTop
sub _rateancestor {
  my ($self) = shift;

  if (@_) {
    $self->{_rateancestor} = shift;
  }

  if (!defined $self->{_rateancestor}){
#    print STDERR "By default, RateAncestor option set to 0 - not" .
# "\nestimating ancestral states.\n";
return 0 } return $self->{_rateancestor} } # What is this value for?
}
_runmodedescriptionprevnextTop
sub _runmode {
  my ($self, $set, $value) = @_;

  if (defined $set) {
    $self->{_runmode} = $value + 1;  # Add one to cope with 0 as a valid option.
} unless (defined $self->{_runmode}){ throw ("Runmode has not been set."); } return $self->{_runmode} - 1 # Remember to remove one to return correct option.
}
_seqfiledescriptionprevnextTop
sub _seqfile {
  my ($self) = shift;

  if (@_) {
    $self->{_seqfile} = shift;
  }

  unless ($self->{_seqfile}) {
    $self->{_seqfile} = $self->work_dir . 'paml_in_' . time;
  }

  return $self->{_seqfile}
}
_seqtypedescriptionprevnextTop
sub _seqtype {
  my ($self) = shift;

  if (@_) {
    $self->{_seqtype} = shift;
  }

  unless (defined $self->{_seqtype}){
    throw ("Seqtype has not been set.");
  } 

  return $self->{_seqtype}
}
_small_diffdescriptionprevnextTop
sub _small_diff {
  my ($self) = shift;

  if (@_) {
    $self->{_small_diff} = shift;
  }

  if (!defined $self->{_small_diff}){
    # Dont know what this does - default set directly from sample file.
return '0.5e-6' } return $self->{_small_diff}
}
_treefiledescriptionprevnextTop
sub _treefile {
  my ($self) = shift;

  if (@_) {
    $self->{_treefile} = shift;
  }

  if (!$self->{_treefile} && $self->_runmode == 0) {
    throw ("Input tree file has not been defined.");
  }

  if (!$self->{_treefile}) {
    $self->{_treefile} = '';
  }

  return $self->{_treefile}
}
_verbosedescriptionprevnextTop
sub _verbose {
  my ($self) = shift;

  if (@_) {
    $self->{_verbose} = shift;
  }

  if (!$self->{_verbose}){
    return 0
  } else {
    return $self->{_verbose}
  }
}
_write_config_filedescriptionprevnextTop
sub _write_config_file {
  my ($self) = @_;

  my $config_string;

  $config_string .= '     seqfile = ' . $self->_seqfile      . "\n";
  $config_string .= '    treefile = ' . $self->_treefile     . "\n";
  $config_string .= '     outfile = ' . $self->_outfile      . "\n";
  $config_string .= '       noisy = ' . $self->_noisy        . "\n";
  $config_string .= '     verbose = ' . $self->_verbose      . "\n";
  $config_string .= '     runmode = ' . $self->_runmode      . "\n";
  $config_string .= '     seqtype = ' . $self->_seqtype      . "\n";
  $config_string .= '   CodonFreq = ' . $self->_codonfreq    . "\n";
  $config_string .= '      aaDist = ' . $self->_aadist       . "\n";
  $config_string .= '  aaRatefile = ' . $self->_aaratefile   . "\n";
  $config_string .= '       model = ' . $self->_model        . "\n";
  $config_string .= '     NSsites = ' . $self->_nssites      . "\n";
  $config_string .= '       icode = ' . $self->_icode        . "\n";
  $config_string .= '       Mgene = ' . $self->_mgene        . "\n";
  $config_string .= '   fix_kappa = ' . $self->_fix_kappa    . "\n";
  $config_string .= '       kappa = ' . $self->_kappa        . "\n";
  $config_string .= '   fix_omega = ' . $self->_fix_omega    . "\n";
  $config_string .= '       omega = ' . $self->_omega        . "\n";
  $config_string .= '   fix_alpha = ' . $self->_fix_alpha    . "\n";
  $config_string .= '       alpha = ' . $self->_alpha        . "\n";
  $config_string .= '      Malpha = ' . $self->_malpha       . "\n";
  $config_string .= '       ncatG = ' . $self->_ncatg        . "\n";
  $config_string .= '       clock = ' . $self->_clock        . "\n";
  $config_string .= '       getSE = ' . $self->_getse        . "\n";
  $config_string .= 'RateAncestor = ' . $self->_rateancestor . "\n";
  $config_string .= '  Small_Diff = ' . $self->_small_diff   . "\n";

  my $file_with_path = $self->work_dir . "paml_codeml.ctl";

  open(CONFIG, ">$file_with_path") or die "Cant open config file for writing.";

  print CONFIG $config_string;

  close(CONFIG);

  $self->_config_file($file_with_path);

  return 1
}


### Files and executables ###
}
_write_seqsdescriptionprevnextTop
sub _write_seqs {
  my ($self, $filename) =  @_;

  throw ("No sequences to write.")
    unless $self->_aligned_seqs;

  system("rm -f $filename");

  open(OUT, ">$filename") or die "Cant write to file [$filename]\n";

  print OUT scalar @{$self->_aligned_seqs} . " " . $self->_aligned_seqs->[0]->length . "\n";

  foreach my $aligned_seq (@{$self->_aligned_seqs}){
    print OUT $aligned_seq->display_id . "\n" . $aligned_seq->seq . "\n";
  }

  close(OUT);

  return 1
}
newdescriptionprevnextTop
sub new {
  my ($class, @args) = @_;

  my $self = bless {},$class;

  my ($executable,
      $work_dir,
      $aligned_seqs,
      $seqfile,
      $treefile,
      $outfile,
      $noisy,
      $verbose,
      $runmode,
      $seqtype,
      $codonfreq,
      $aadist,
      $aaratefile,
      $model,
      $nssites,
      $icode,
      $mgene,
      $fix_kappa,
      $kappa,
      $fix_omega,
      $omega,
      $fix_alpha,
      $alpha,
      $malpha,
      $ncatg,
      $clock,
      $getse,
      $rateancestor,
      $small_diff) = rearrange([qw(EXECUTABLE
				   WORK_DIR
				   ALIGNED_SEQS
				   SEQFILE
				   TREEFILE
				   OUTFILE
				   NOISY
				   VERBOSE
				   RUNMODE
				   SEQTYPE
				   CODONFREQ
				   AADIST
				   AARATEFILE
				   MODEL
				   NSSITES
				   ICODE
				   MGENE
				   FIX_KAPPA
				   KAPPA
				   FIX_OMEGA
				   OMEGA
				   FIX_ALPHA
				   ALPHA
				   MALPHA
				   NCATG
				   CLOCK
				   GETSE
				   RATEANCESTOR
				   SMALL_DIFF)],@args);

  throw ("Have sequence input and input file - which do I use?") 
    if (defined $aligned_seqs && defined $seqfile);
  
  $self->_aligned_seqs($aligned_seqs) if $aligned_seqs;

  $self->_codeml_executable($executable) if $executable;
  $self->work_dir($work_dir)             if $work_dir;
  $self->_seqfile($seqfile)              if $seqfile;
  $self->_treefile($treefile)            if $treefile;

  $self->_pre_copying; # Set up our work dir the way we like it.
$self->_outfile($outfile) if $outfile; $self->_noisy($noisy) if $noisy; $self->_verbose($verbose) if $verbose; $self->_runmode('set', $runmode); $self->_seqtype($seqtype) if $seqtype; $self->_codonfreq($codonfreq) if $codonfreq; $self->_aadist($aadist) if $aadist; $self->_aaratefile($aaratefile) if $aaratefile; $self->_model('set', $model); $self->_nssites($nssites); $self->_icode($icode); $self->_treefile($mgene) if $mgene; $self->_treefile($fix_kappa) if $fix_kappa; $self->_treefile($kappa) if $kappa; $self->_treefile($fix_omega) if $fix_omega; $self->_treefile($omega) if $omega; $self->_treefile($fix_alpha) if $fix_alpha; $self->_treefile($alpha) if $alpha; $self->_treefile($malpha) if $malpha; $self->_treefile($ncatg) if $ncatg; $self->_treefile($clock) if $clock; $self->_treefile($getse) if $getse; $self->_treefile($rateancestor) if $rateancestor; $self->_treefile($small_diff) if $small_diff; return $self
}
run_codemldescriptionprevnextTop
sub run_codeml {
  my $self = shift;

  # Directory changing shenanigans to make PAML behave!
my $user_dir = cwd(); # Move to our working dir
warning ($!) if (! chdir $self->work_dir); $self->_write_seqs($self->_seqfile); $self->_write_config_file; my $command = $self->_codeml_executable . " " . $self->_config_file; print STDERR $command . "\n"; eval { system($command) }; if ($@ or -e $self->work_dir . "/core") { throw ("Something went wrong when codeml was executed.\n" . $@); } # Move back to original dir
print STDERR $! if (! chdir $user_dir); my $parser = Bio::Tools::Phylo::PAML->new('-file' => $self->_outfile, '-dir' => $self->work_dir); return $parser; } # Make a temporary playground for PAML and some of its anti-social ways.
# Execute application from here such that all unwanted temp files can be
# neatly deleted. Also, can copy assorted external files, like grantham.dat
# here before execution. Give temp dir a unique name to avoid conflicts
# from multiple simultaneous jobs.
}
work_dirdescriptionprevnextTop
sub work_dir {
  my ($self) = shift;

  unless ($self->{_work_dir}) {

    my $work_dir;

    if (@_) {
      my $arg_dir = shift;
      $work_dir = $arg_dir . '/' . 'paml_temp_' . time . '/';
    } else {
      $work_dir = '/tmp/paml_temp_' . time . '/';
    }

    unless (-d $work_dir){
      mkdir $work_dir;
    }

    $self->{_work_dir} = $work_dir;
  }

  return $self->{_work_dir}
}
General documentation
No general documentation available.