<!DOCTYPE article PUBLIC "-//OASIS//DTD DocBook V4.1//EN" [
<!ENTITY paml "<acronym>PAML</acronym>" >
]>
<!-- $Id: PAML.sgml,v 1.4 2003/09/17 12:46:51 bosborne Exp $ -->
<article>
<articleinfo>

<title>Bio::Tools::Phylo::PAML HOWTO</title>
<author>
	<firstname>Aaron</firstname>
	<surname>Mackey</surname>
	<affiliation>
		<orgname>
			<ulink url="http://www.virginia.edu">University of Virginia</ulink>
		</orgname>
		<address>
			<email>amackey@virginia.edu</email>
		</address>
	</affiliation>
</author>
<author>
	<firstname>Jason</firstname>
	<surname>Stajich</surname>
      <affiliation>
	<orgname><ulink url="http://www.duke.edu">Duke
	University</ulink></orgname>
	<orgdiv><ulink url="http://upg.duke.edu">University Program in Genetics</ulink></orgdiv>
	<orgdiv><ulink url="http://cgt.genetics.duke.edu">Center for
	Genome Technology</ulink></orgdiv>
	<address>
	  Duke University Medical Center
	  <pob>Box 3568</pob>
	  <city>Durham, </city>
	  <state>North Carolina</state> 
	  <postcode>27710-3568</postcode> 
	  <country>USA</country>
	  <email>jason-at-bioperl.org</email>
        </address>
	</affiliation>
</author>

<pubdate>2002-08-01</pubdate>
 
<revhistory>
	<revision>
		<revnumber>0.1</revnumber>
		<date>2002-08-01</date>
		<authorinitials>ajm</authorinitials>
		<revremark>first draft</revremark>
	</revision>
	<revision>
	        <revnumber>0.2</revnumber>
		<date>2003-03-01</date>
		<authorinitials>jes</authorinitials>
		<revremark>Added pairwise Ka,Ks example code and running code</revremark>
	</revision>
</revhistory>

<legalnotice>
<para>
	This document is copyright Aaron Mackey, 2002.  For
	reproduction other than personal use please contact me at
	amackey@virginia.edu
</para>
</legalnotice>

<abstract>
<para>
	<emphasis>paml</emphasis> is a package of C programs that implement
	<emphasis>P</emphasis>hylogenetic
	<emphasis>A</emphasis>nalyses using
	<emphasis>M</emphasis>aximum
	<emphasis>L</emphasis>ikelihood, written by Dr. Ziheng Yang,
	University College London.  These programs implement a wide
	variety of models to explore the evolutionary relationships
	between sequences at either the protein, codon or raw DNA
	level.  This document's aim is to explore and document how the
	BioPerl <emphasis>paml</emphasis> parser and result objects "work".
</para>
</abstract>

</articleinfo>

<section id="background">
<title>Background</title>
<para>
	The <emphasis>paml</emphasis> package consists of many different executable
	programs, but the BioPerl Bio::Tools::Phylo::PAML object
	(hereafter referred to as simply <emphasis>the</emphasis> <emphasis>paml</emphasis>
	object) focuses on dealing with the output of the main
	analysis programs "baseml", "codeml" (sometimes called "aaml")
	and "codemlsites" (a batch version of "codeml").  All of these
	programs use maximum likelihood methods to fit a mathematical
	model of evolution to sequence data provided by the user.  The
	main difference between these programs is the type of sequence
	on which they operate (baseml for raw DNA, codeml for DNA
	organized as codons, aaml for amino acids).
</para>

<para>
	While the general maximum likelihood approach used by the <emphasis>paml</emphasis>
	programs is the same for all of them, the specific
	evolutionary models available for each sequence type vary
	greatly, as do the parameters specific to each model.  The
	programs function in a handful of disparate modes, each
	requiring slight variations of inputs that can possibly
	include:

	<orderedlist>
	<listitem>
	<para>
	multiply-aligned sequences. representing 1 or more distinct
	genes [ <emphasis>paml</emphasis> parameter <emphasis>Mgene</emphasis> = 1 ], in 1
	or more distinct datasets [ <emphasis>paml</emphasis> <emphasis>ndata</emphasis>
	&gt; 1 ])
	</para>
	</listitem>

	<listitem>
	<para>
	a user-provided tree topology (or multiple tree topologies to
	be evaluated and contrasted)
	</para>
	</listitem>

	<listitem>
	<para>
	a set of instructions in a control file that specify the model
	(or models) to be used, various options to specify how to
	handle the sequence data (e.g. whether to dismiss columns with
	gaps or not [ <emphasis>cleandata</emphasis> = 1 ]), initial
	or fixed values for model parameters, and the filenames for
	other input data.
	</para>
	</listitem>
	</orderedlist>

	The output from <emphasis>paml</emphasis> is directed to multiple "targets": data
	is written to the user-specified primary output file
	(conventionally named with an .mlc extension), as well as
	various accessory files with fixed names
	(e.g. <filename>2ML.t</filename>, <filename>2ML.dN</filename>,
	<filename>2ML.dS</filename> for pairwise Maximum Likelihood
	calculations) that appear in the same directory that the
	output file is found.
</para>

<para>
	The upshot of these comments is that one <emphasis>paml</emphasis> program "run"
	can potentially generate results for many genes, many
	datasets, many tree toplogies and many evolutionary models,
	spread across multiple output files.  Currently, the <emphasis>paml</emphasis>
	programs deal with the various categories of multiple analyses
	in the following "top-down" order: datasets, genes, models,
	tree topologies.  So how shall the BioPerl <emphasis>paml</emphasis> module treat
	these sources of multiple results?
</para>
</section>

<section id="results">
<title>Accessing <emphasis>paml</emphasis> results</title>
<para>
	The BioPerl <emphasis>paml</emphasis> result parser takes the view that a distinct
	"recordset" or single, top-level PAML::Result object
	represents a single <emphasis>dataset</emphasis>.  Each
	PAML::Result object may therefore contain data from multiple
	genes, models, and/or tree topologies.  To parse the output
	from a multiple-dataset <emphasis>paml</emphasis> run, the familiar "next_result"
	iterator common to other BioPerl modules is invoked.

<example id="nextres">
<title>Iterating over results with next_result</title>
<programlisting>
use Bio::Tools::Phylo::PAML;

my $parser = new Bio::Tools::Phylo::PAML (-file => "./output.mlc",
                                          -dir  => "./",
                                          -ctlf => "./codeml.ctl");

while(my $result = $parser->next_result) {
    # do something with the results from this dataset ...
}
</programlisting>
</example>

	In this example, we've created a new top-level <emphasis>paml</emphasis> parser,
	specifying <emphasis>paml</emphasis>'s primary output file, the directory in which
	any other accessory files may be found, and the control file.
	We then trigger the parser to begin parsing the data,
	returning a new PAML::Result object for each dataset found in
	the output.
</para>

<para>
	The PAML::Result object provides access to the wide variety of
	data found in the output files.  The specific kinds of data
	available depends on which <emphasis>paml</emphasis> analysis program was run,
	and the mode and models employed.  Generally, these include a
	recapitulation of the input sequences and their multiple
	alignment (which may differ slightly from the original input
	sequences due to the data "cleansing" <emphasis>paml</emphasis> performs),
	descriptive statistics of the input sequences (e.g. codon
	usage tables, nucleotide or amino acid composition), pairwise
	Nei and Gojobori (NG) calculation matrices (for codon models),
	fitted model parameter values (including branch-specific
	parameters associated with any provided tree topology),
	reconstructed ancestral sequences (again, associated with an
	accompanying tree topology), or statistical comparisons of
	multiple tree topologies.
</para>

</section>

<section id="running">
<title>Running PAML from within Bioperl</title>
<para>

       Bioperl also has facilities for runnning <emphasis>paml</emphasis> from within a
       Perl script.  This allows you to compute Ka and Ks estimations
       from within an analysis pipeline.  The following section will
       describe the process of getting data into Bioperl, running the
       alignment process, and setting up a <emphasis>paml</emphasis> process.  This code
       is focusing on estimations of all the pairwise Ka and Ks values
       however it can be used to easily compute more sophisticated
       questions about variable rates, etc. 
              
</para>
<para>

       This code below is an excerpt from 
       scripts/utilities/pairwise_kaks.PLS which will calculate all
       pairwise Ka,Ks values for a set of cDNA sequences stored in a
       file.  It will first translate the cDNA into protein and align
       the protein sequences.  This is a simple way to insure gaps
       only occur at codon boundaries and amino acid substitution
       rates are applied when calculating the MSA.  The protein
       alignment is them projected back into cDNA coordinates using a
       method called aa_to_dna_aln.  Finally the cDNA alignment is
       provided to a <emphasis>paml</emphasis> executing module which sets up the running
       parameters and converts the alignment to the appropriate format.

 
<programlisting>

use Bio::Tools::Run::Phylo::PAML::Codeml;
use Bio::Tools::Run::Alignment::Clustalw;

# for projecting alignments from protein to R/DNA space
use Bio::Align::Utilities qw(aa_to_dna_aln);

# for input of the sequence data
use Bio::SeqIO;
use Bio::AlignIO;

my $aln_factory = new Bio::Tools::Run::Alignment::Clustalw();

my $seqdata = 'cdna.fa';

my $seqIO = new Bio::SeqIO(-file   => $seqdata,
                           -format => 'fasta');
my %seqs;
my @prots;
# process each sequence
while( my $seq = $seqin->next_seq ) {
    $seqs{$seq->display_id} = $seq;
    # translate them into protein
    my $protein = $seq->translate();
    my $pseq = $protein->seq();
    if( $pseq =~ /\*/ &amp;&amp; 
	$pseq !~ /\*$/ ) {
	warn("provided a cDNA sequence with a stop codon, PAML will choke!");
	exit(0);
    }
    # Tcoffee can't handle '*' even if it is trailing
    $pseq =~ s/\*//g;
    $protein->seq($pseq);
    push @prots, $protein;
}

if( @prots &lt; 2 ) {
    warn("Need at least 2 cDNA sequences to proceed");
    exit(0);
}

open(OUT, ">align_output.txt") || 
      die("cannot open output $output for writing");
# Align the sequences with clustalw
my $aa_aln = $aln_factory->align(\@prots);
# project the protein alignment back to cDNA coordinates
my $dna_aln = &amp;aa_to_dna_aln($aa_aln, \%seqs);

my @each = $dna_aln->each_seq();

my $kaks_factory = new Bio::Tools::Run::Phylo::PAML::Codeml
                  ( -params => { 'runmode' => -2,
		                 'seqtype' => 1,
		               }
	          );

# set the alignment object
$kaks_factory->alignment($dna_aln);

# run the KaKs analysis
my ($rc,$parser) = $kaks_factory->run();
my $result = $parser->next_result;
my $MLmatrix = $result->get_MLmatrix();

my @otus = $result->get_seqs();
# this gives us a mapping from the PAML order of sequences back to 
# the input order (since names get truncated)
my @pos = map { 
    my $c= 1;
    foreach my $s ( @each ) {
	last if( $s->display_id eq $_->display_id );
	$c++;
    }
    $c; 
} @otus; 

print OUT join("\t", qw(SEQ1 SEQ2 Ka Ks Ka/Ks PROT_PERCENTID CDNA_PERCENTID)), "\n";
for( my $i = 0; $i &lt; (scalar @otus -1) ; $i++) {
    for( my $j = $i+1; $j &lt; (scalar @otus); $j++ ) {
	my $sub_aa_aln  = $aa_aln->select_noncont($pos[$i],$pos[$j]);
	my $sub_dna_aln = $dna_aln->select_noncont($pos[$i],$pos[$j]);
	print OUT join("\t",  
		       $otus[$i]->display_id,
		       $otus[$j]->display_id,$MLmatrix->[$i]->[$j]->{'dN'},
		       $MLmatrix->[$i]->[$j]->{'dS'},
		       $MLmatrix->[$i]->[$j]->{'omega'},
		       sprintf("%.2f",$sub_aa_aln->percentage_identity),
		       sprintf("%.2f",$sub_dna_aln->percentage_identity),
		       ), "\n";
    }
}

</programlisting>
</para>
<para>

</para>

</section>

</article>
