Bio::Tools::Phylo::PAML HOWTO

The paml package consists of many different executable programs, but the BioPerl Bio::Tools::Phylo::PAML object (hereafter referred to as simply the paml 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).

While the general maximum likelihood approach used by the paml 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:

The upshot of these comments is that one paml 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 paml 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 paml module treat these sources of multiple results?

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

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 ...
}
	

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 paml 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" paml 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.

Bioperl also has facilities for running paml 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 paml 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.

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 paml executing module which sets up the running parameters and converts the alignment to the appropriate format.

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 = $seqio->next_seq ) {
    $seqs{$seq->display_id} = $seq;
    # translate them into protein
    my $protein = $seq->translate();
    my $pseq = $protein->seq();
    if( $pseq =~ /\*/ && 
	$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 < 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 = &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 < (scalar @otus -1) ; $i++) {
    for( my $j = $i+1; $j < (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";
    }
}