[Bioperl-l] PAML + Codeml problem..

[Xianjun Dong]

Hi, 

I have a problem during running the Codeml Wiki-HOWTO code:

Here is the error message:
////////////////////////////////////////////////////////////////
[xianjund at lauvtre kaks]$ perl paml.pl test.fa

-------------------- WARNING ---------------------
MSG: There was an error - see error_string for the program output
STACK
Bio::Tools::Run::Phylo::PAML::Codeml::run /Home/extern/xianjund/src/bioperl/bioperl-run/Bio/Tools/Run/Phylo/PAML/Codeml.pm:581
STACK toplevel paml.pl:61

------------- EXCEPTION: Bio::Root::NotImplemented -------------
MSG: Unknown format of PAML output
STACK: Error::throw
STACK:
Bio::Root::Root::throw /usr/lib/perl5/site_perl/5.8.5/Bio/Root/Root.pm:328
STACK:
Bio::Tools::Phylo::PAML::_parse_summary /usr/lib/perl5/site_perl/5.8.5/Bio/Tools/Phylo/PAML.pm:359
STACK:
Bio::Tools::Phylo::PAML::next_result /usr/lib/perl5/site_perl/5.8.5/Bio/Tools/Phylo/PAML.pm:224
STACK: paml.pl:62
----------------------------------------------------------------
////////////////////////////////////////////////////////////////

My test sequence is:
>ENST00000361390
ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCGAACGAAAAATTCTAGGCTATATACAACTACGCAAAGGCCCCAACGTTGTAGGCCCCTACGGGCTACTACAACCCTTCGCTGACGCCATAAAACTCTTCACCAAAGAGCCCCTAAAACCCGCCACATCTACCATCACCCTCTACATCACCGCCCCGACCTTAGCTCTCACCATCGCTCTTCTACTATGAACCCCCCTCCCCATACCCAACCCCCTGGTCAACCTCAACCTAGGCCTCCTATTTATTCTAGCCACCTCTAGCCTAGCCGTTTACTCAATCCTCTGATCAGGGTGAGCATCAAACTCAAACTACGCCCTGATCGGCGCACTGCGAGCAGTAGCCCAAACAATCTCATATGAAGTCACCCTAGCCATCATTCTACTATCAACATTACTAATAAGTGGCTCCTTTAACCTCTCCACCCTTATCACAACACAAGAACACCTCTGATTACTCCTGCCATCATGACCCTTGGCCATAATATGATTTATCTCCACACTAGCAGAGACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAATACGCCGCAGGCCCCTTCGCCCTATTCTTCATAGCCGAATACACAAACATTATTATAATAAACACCCTCACCACTACAATCTTCCTAGGAACAACATATGACGCACTCTCCCCTGAACTCTACACAACATATTTTGTCACCAAGACCCTACTTCTAACCTCCCTGTTCTTATGAATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACCTCCTATGAAAAAACTTCCTACCACTCACCCTAGCATTACTTATATGATATGTCTCCATACCCATTACAATCTCCAGCATTCCCCCTCAAACCTAA
>ENSMUST00000082392
GTGTTCTTTATTAATATCCTAACACTCCTCGTCCCCATTCTAATCGCCATAGCCTTCCTAACATTAGTAGAACGCAAAATCTTAGGGTACATACAACTACGAAAAGGCCCTAACATTGTTGGTCCATACGGCATTTTACAACCATTTGCAGACGCCATAAAATTATTTATAAAAGAACCAATACGCCCTTTAACAACCTCTATATCCTTATTTATTATTGCACCTACCCTATCACTCACACTAGCATTAAGTCTATGAGTTCCCCTACCAATACCACACCCATTAATTAATTTAAACCTAGGGATTTTATTTATTTTAGCAACATCTAGCCTATCAGTTTACTCCATTCTATGATCAGGATGAGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGTAGCCCAAACAATTTCATATGAAGTAACCATAGCTATTATCCTTTTATCAGTTCTATTAATAAATGGATCCTACTCTCTACAAACACTTATTACAACCCAAGAACACATATGATTACTTCTGCCAGCCTGACCCATAGCCATAATATGATTTATCTCAACCCTAGCAGAAACAAACCGGGCCCCCTTCGACCTGACAGAAGGAGAATCAGAATTAGTATCAGGGTTTAACGTAGAATACGCAGCCGGCCCATTCGCGTTATTCTTTATAGCAGAGTACACTAACATTATTCTAATAAACGCCCTAACAACTATTATCTTCCTAGGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACTAACTTCATAATAGAAGCTCTACTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATCTTCTATGAAAAAACTTTCTACCCCTAACACTAGCATTATGTATGTGACATATTTCTTTACCAATTTTTACAGCGGGAGTACCACCATACATATAG

Sure, I checked it. There is some stop codon in it. If I replace it with
non-stop codon, it works.

For example, 
>ENST00000361390
ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCcaaTCGCAATGGCATTCCcaaTGCTTACCGAACGAAAAATTCcaaGCTATATACAACTACGCAAAGGCCCCAACGTTGcaaGCCCCTACGGGCTACTACAACCCTTCGCcaaCGCCAcaaAACTCTTCACCAAAGAGCCCCcaaAACCCGCCACATCTACCATCACCCTCTACATCACCGCCCCGACCTcaaCTCTCACCATCGCTCTTCTACTAcaaACCCCCCTCCCCATACCCAACCCCCTGGTCAACCTCAACCcaaGCCTCCTATTTATTCcaaCCACCTCcaaCCcaaCCGTTTACTCAATCCTCcaaTCAGGGcaaGCATCAAACTCAAACTACGCCCcaaTCGGCGCACTGCGAGCAGcaaCCCAAACAATCTCATAcaaAGTCACCCcaaCCATCATTCTACTATCAACATTACcaacaaGTGGCTCCTTcaaCCTCTCCACCCTTATCACAACACAAGAACACCTCcaaTTACTCCTGCCATCAcaaCCCTTGGCCAcaaTAcaaTTTATCTCCACACcaaCAGAGACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACcaaTCTCAGGCTTCAACATCGAATACGCCGCAGGCCCCTTCGCCCTATTCTTCAcaaCCGAATACACAAACATTATTAcaacaaACACCCTCACCACTACAATCTTCCcaaGAACAACATAcaaCGCACTCTCCCCcaaACTCTACACAACATATTTTGTCACCAAGACCCTACTTCcaaCCTCCCTGTTCTTAcaaATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACCTCCTAcaaAAAAACTTCCTACCACTCACCCcaaCATTACTTATAcaaTATGTCTCCATACCCATTACAATCTCCAGCATTCCCCCTCAAACCcaa
>ENSMUST00000082392
GTGTTCTTTATcaaTATCCcaaCACTCCTCGTCCCCATTCcaaTCGCCAcaaCCTTCCcaaCATcaacaaAACGCAAAATCTcaaGGTACATACAACTACGAAAAGGCCCcaaCATTGTTGGTCCATACGGCATTTTACAACCATTTGCAGACGCCAcaaAATTATTTAcaaAAGAACCAATACGCCCTTcaaCAACCTCTATATCCTTATTTATTATTGCACCTACCCTATCACTCACACcaaCATcaaGTCTAcaaGTTCCCCTACCAATACCACACCCATcaaTcaaTTcaaACCcaaGGATTTTATTTATTTcaaCAACATCcaaCCTATCAGTTTACTCCATTCTAcaaTCAGGAcaaGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGcaaCCCAAACAATTTCATAcaaAGcaaCCAcaaCTATTATCCTTTTATCAGTTCTATcaacaaATGGATCCTACTCTCTACAAACACTTATTACAACCCAAGAACACATAcaaTTACTTCTGCCAGCCcaaCCCAcaaCCAcaaTAcaaTTTATCTCAACCCcaaCAGAAACAAACCGGGCCCCCTTCGACCcaaCAGAAGGAGAATCAGAATcaaTATCAGGGTTcaaCGcaaAATACGCAGCCGGCCCATTCGCGTTATTCTTTAcaaCAGAGTACACcaaCATTATTCcaacaaACGCCCcaaCAACTATTATCTTCCcaaGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACcaaCTTCAcaacaaAAGCTCTACTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATCTTCTAcaaAAAAACTTTCTACCCCcaaCACcaaCATTATGTATGcaaCATATTTCTTTACCAATTTTTACAGCGGGAGTACCACCATACATAcaa

But my question is: it does not occur in the codon position (say, the
third codon's position is not a times of 3). Why it effect the result?

And also there is code to filter out the stop codon in the sample code
(as the following shown)
///////////////////////////////
    if( $pseq =~ /\*/ &&
        $pseq !~ /\*$/ ) {
          warn("provided a CDS sequence with a stop codon, PAML will
choke!");
          exit(0);
    }
    # Tcoffee can't handle '*' even if it is trailing
    $pseq =~ s/\*//g;
/////////////////////////////

So, when translate back from aa_aln to dna_aln, there should be no stop
codon included. SO, why it can not pass? 

Thanks for answer!

P.S: attach my code here:
/////////////////////////////////////////////////////////
#!/usr/bin/perl -w
use strict;
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 = Bio::Tools::Run::Alignment::Clustalw->new('quiet'=>1);
my $seqdata = shift || 'test.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 CDS 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 CDS sequences to proceed");
    exit(0);
}
 
# open(OUT, ">align_output.txt") ||  die("cannot open output
align_output for writing");
# Align the sequences with clustalw
my $aa_aln = $aln_factory->align(\@prots);
# project the protein alignment back to CDS coordinates
my $dna_aln = aa_to_dna_aln($aa_aln, \%seqs);
 
my @each = $dna_aln->each_seq();
 
my $kaks_factory = Bio::Tools::Run::Phylo::PAML::Codeml->new
                   ( -params => { 'runmode' => -2,
                                  'seqtype' => 1,
                                },
		     -save_tempfiles => 1,
		     -verbose => 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 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 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";
  }
}

-- 
Xianjun Dong 
PhD Student 
Computational Biology Unit 
Bergen Center for Computational Science 
University of Bergen
Høyteknologisenteret, Thormøhlensgate 55 
N-5008 Bergen,Norway. 
Webpage: http://www.ii.uib.no/~xianjund/ 
MSN: sterding at hotmail.com
Phone No: +47 - 55584354 (office)
          +47 - 47361688 (mobile)
Fax No: +47 - 55584295