[Bioperl-l] PAML + Codeml problem..

Ryan Golhar golharam at umdnj.edu
Thu Aug 10 14:53:51 EDT 2006


Hi Xianjun,

1.  The Bio::Seq::translate function (to my knowledge) only uses the
generic codon table.  So, you will need to translate the DNA sequence
using some other method.  In any case, even removing the *'s from the
protein sequence still leaves the stop codons in the DNA sequence which
must be removed.

2.  The checks were written to assume that the sequences provided are
full-length coding sequences.  That means the start and stop codon are
present as well.  When the translate function is called, the stop codon
is translated as a '*'.  The script initally just remove the * from the
end of the sequence and continued on.  

I added a check to see if there is a '*' in the middle of the sequence
because I found in some of my genes that there is in fact in-frame stop
codons which actually codes for selenocysteine.  I see the warning check
isn't working for some reason - odd, it worked when I wrote it.  

You can remove the *'s from the protein sequence, but you must also be
sure to remove the corresponding codons from the DNA sequence as well
before invoking run() on the Codeml pacakge.  I suppose someone could
add a check to the script to remove the in-frame stop codons.
Remember, the pairwise_kaks script is just a starting point (tutorial)
to show you how you can go about performing this type of an analysis.  

In fact, I've since switched from PAML to using a different method PBL
which a colleuge coded.  I found that PAML tends to overestimate
synonymous rates in some cases.  

Let me know if this helps.  If not, I'll try to clarify.  

Ryan

-----Original Message-----
From: Xianjun Dong [mailto:xianjun.dong at bccs.uib.no] 
Sent: Thursday, August 10, 2006 12:03 PM
To: golharam at umdnj.edu
Cc: bioperl-l at lists.open-bio.org
Subject: RE: [Bioperl-l] PAML + Codeml problem..


Hi, Ryan

Thanks for your reply!

But here I still have two questions about the sample code:
1. the translate() function of Bio::Seq object use generic codon table,
but for Mitochondrial DNA (mtDNA), we should use different codon table.
So, if we take the human transcript ENST00000361390 as example, 

>ENST00000361390 _cDNA
ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCGAA
CGAAAAATTCTAGGCTATATACAACTACGCAAAGGCCCCAACGTTGTAGGCCCCTACGGGCTACTACAACCC
TTCGCTGACGCCATAAAACTCTTCACCAAAGAGCCCCTAAAACCCGCCACATCTACCATCACCCTCTACATC
ACCGCCCCGACCTTAGCTCTCACCATCGCTCTTCTACTATGAACCCCCCTCCCCATACCCAACCCCCTGGTC
AACCTCAACCTAGGCCTCCTATTTATTCTAGCCACCTCTAGCCTAGCCGTTTACTCAATCCTCTGATCAGGG
TGAGCATCAAACTCAAACTACGCCCTGATCGGCGCACTGCGAGCAGTAGCCCAAACAATCTCATATGAAGTC
ACCCTAGCCATCATTCTACTATCAACATTACTAATAAGTGGCTCCTTTAACCTCTCCACCCTTATCACAACA
CAAGAACACCTCTGATTACTCCTGCCATCATGACCCTTGGCCATAATATGATTTATCTCCACACTAGCAGAG
ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAATACGCC
GCAGGCCCCTTCGCCCTATTCTTCATAGCCGAATACACAAACATTATTATAATAAACACCCTCACCACTACA
ATCTTCCTAGGAACAACATATGACGCACTCTCCCCTGAACTCTACACAACATATTTTGTCACCAAGACCCTA
CTTCTAACCTCCCTGTTCTTATGAATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACCTC
CTATGAAAAAACTTCCTACCACTCACCCTAGCATTACTTATATGATATGTCTCCATACCCATTACAATCTCC
AGCATTCCCCCTCAAACCTAA

After translating with above function, the amino acid sequence is like
this, which contain *(stop codon) within the sequence(also at the end of
the sequence). But actually, this is a mtDNA, if we use different codon
table, the * within the sequence will change to 'W'(Trp). (Because in
vertebrate mitochondria "AGA" and "AGG" are also stop codons, but not
"UGA", which codes for tryptophan instead.)
>ENST00000361390 aa_beforefilter
IPMANLLLLIVPILIAMAFLMLTERKILGYIQLRKGPNVVGPYGLLQPFADAIKLFTKEPLKPATSTITLYI
TAPTLALTIALLL*TPLPIPNPLVNLNLGLLFILATSSLAVYSIL*SG*ASNSNYALIGALRAVAQTISYEV
TLAIILLSTLLISGSFNLSTLITTQEHL*LLLPS*PLAII*FISTLAETNRTPFDLAEGESELVSGFNIEYA
AGPFALFFIAEYTNIIIINTLTTTIFLGTTYDALSPELYTTYFVTKTLLLTSLFL*IRTAYPRFRYDQLIHL
L*KNFLPLTLALLI*YVSIPITISSIPPQT*

2. My second question is: 
If there are * both in the middle and end of the translated sequence
(with pattern AAAAAA*AAAAAAAAAAAAAAA*AAA*), like above case, after the
two checks for stop codon, all * will be filtered out. So, when
translate back from aa_aln to dna_aln, there should be no stop codon
included. But actually, when I track the program, it display that there
are still stop codon included. Here is the DNA alignment after recalling
the aa_to_dna_aln function. How to explain this?
 
>ENST00000361390 aa_to_dna_aln
ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCGAA
CGAAAAATTCTAGGCTATATACAACTACGCAAAGGCCCCAACGTTGTAGGCCCCTACGGGCTACTACAACCC
TTCGCTGACGCCATAAAACTCTTCACCAAAGAGCCCCTAAAACCCGCCACATCTACCATCACCCTCTACATC
ACCGCCCCGACCTTAGCTCTCACCATCGCTCTTCTACTATGAACCCCCCTCCCCATACCCAACCCCCTGGTC
AACCTCAACCTAGGCCTCCTATTTATTCTAGCCACCTCTAGCCTAGCCGTTTACTCAATCCTCTGATCAGGG
TGAGCATCAAACTCAAACTACGCCCTGATCGGCGCACTGCGAGCAGTAGCCCAAACAATCTCATATGAAGTC
ACCCTAGCCATCATTCTACTATCAACATTACTAATAAGTGGCTCCTTTAACCTCTCCACCCTTATCACAACA
CAAGAACACCTCTGATTACTCCTGCCATCATGACCCTTGGCCATAATATGATTTATCTCCACACTAGCAGAG
ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAATACGCC
GCAGGCCCCTTCGCCCTATTCTTCATAGCCGAATACACAAACATTATTATAATAAACACCCTCACCACTACA
ATCTTCCTAGGAACAACATATGACGCACTCTCCCCTGAACTCTACACAACATATTTTGTCACCAAGACCCTA
CTT---CTAACCTCCCTGTTCTTATGAATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACAC
CTCCTATGAAAAAACTTCCTACCACTCACCCTAGCATTACTTATATGATATGTCTCCATACCCATT


I attached my script for two ortholog transcripts demo and the output
(including the error msg) here. Could you kindly check for me?

Thanks!

-Xianjun

/////////////////////////////////////////////////////////////////////
/////////////////////////////// output //////////////////////////////
/////////////////////////////////////////////////////////////////////

[xianjund at lauvtre kaks]$ perl calculator.pl
>ENST00000361390
ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCGAA
CGAAAAATTCTAGGCTATATACAACTACGCAAAGGCCCCAACGTTGTAGGCCCCTACGGGCTACTACAACCC
TTCGCTGACGCCATAAAACTCTTCACCAAAGAGCCCCTAAAACCCGCCACATCTACCATCACCCTCTACATC
ACCGCCCCGACCTTAGCTCTCACCATCGCTCTTCTACTATGAACCCCCCTCCCCATACCCAACCCCCTGGTC
AACCTCAACCTAGGCCTCCTATTTATTCTAGCCACCTCTAGCCTAGCCGTTTACTCAATCCTCTGATCAGGG
TGAGCATCAAACTCAAACTACGCCCTGATCGGCGCACTGCGAGCAGTAGCCCAAACAATCTCATATGAAGTC
ACCCTAGCCATCATTCTACTATCAACATTACTAATAAGTGGCTCCTTTAACCTCTCCACCCTTATCACAACA
CAAGAACACCTCTGATTACTCCTGCCATCATGACCCTTGGCCATAATATGATTTATCTCCACACTAGCAGAG
ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAATACGCC
GCAGGCCCCTTCGCCCTATTCTTCATAGCCGAATACACAAACATTATTATAATAAACACCCTCACCACTACA
ATCTTCCTAGGAACAACATATGACGCACTCTCCCCTGAACTCTACACAACATATTTTGTCACCAAGACCCTA
CTTCTAACCTCCCTGTTCTTATGAATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACCTC
CTATGAAAAAACTTCCTACCACTCACCCTAGCATTACTTATATGATATGTCTCCATACCCATTACAATCTCC
AGCATTCCCCCTCAAACCTAA
>ENSMUST00000082392
GTGTTCTTTATTAATATCCTAACACTCCTCGTCCCCATTCTAATCGCCATAGCCTTCCTAACATTAGTAGAA
CGCAAAATCTTAGGGTACATACAACTACGAAAAGGCCCTAACATTGTTGGTCCATACGGCATTTTACAACCA
TTTGCAGACGCCATAAAATTATTTATAAAAGAACCAATACGCCCTTTAACAACCTCTATATCCTTATTTATT
ATTGCACCTACCCTATCACTCACACTAGCATTAAGTCTATGAGTTCCCCTACCAATACCACACCCATTAATT
AATTTAAACCTAGGGATTTTATTTATTTTAGCAACATCTAGCCTATCAGTTTACTCCATTCTATGATCAGGA
TGAGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGTAGCCCAAACAATTTCATATGAAGTA
ACCATAGCTATTATCCTTTTATCAGTTCTATTAATAAATGGATCCTACTCTCTACAAACACTTATTACAACC
CAAGAACACATATGATTACTTCTGCCAGCCTGACCCATAGCCATAATATGATTTATCTCAACCCTAGCAGAA
ACAAACCGGGCCCCCTTCGACCTGACAGAAGGAGAATCAGAATTAGTATCAGGGTTTAACGTAGAATACGCA
GCCGGCCCATTCGCGTTATTCTTTATAGCAGAGTACACTAACATTATTCTAATAAACGCCCTAACAACTATT
ATCTTCCTAGGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACTAACTTCATAATAGAAGCTCTA
CTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATCTT
CTATGAAAAAACTTTCTACCCCTAACACTAGCATTATGTATGTGACATATTTCTTTACCAATTTTTACAGCG
GGAGTACCACCATACATATAG

Calculate the Ka/Ks for ENSG00000198888 : ENSMUSG00000064341 ...
>ENSMUST00000082392 aa_beforefilter
VFFINILTLLVPILIAIAFLTLVERKILGYIQLRKGPNIVGPYGILQPFADAIKLFIKEPIRPLTTSISLFI
IAPTLSLTLALSL*VPLPIPHPLINLNLGILFILATSSLSVYSIL*SG*ASNSKYSLFGALRAVAQTISYEV
TIAIILLSVLLINGSYSLQTLITTQEHI*LLLPA*PIAII*FISTLAETNRAPFDLTEGESELVSGFNVEYA
AGPFALFFIAEYTNIILINALTTIIFLGPLYYINLPELYSTNFIIEALLLSSTFLWIRASYPRFRYDQLIHL
L*KNFLPLTLALCM*HISLPIFTAGVPPYI*
>ENSMUST00000082392 aa_afterfilter
VFFINILTLLVPILIAIAFLTLVERKILGYIQLRKGPNIVGPYGILQPFADAIKLFIKEPIRPLTTSISLFI
IAPTLSLTLALSLVPLPIPHPLINLNLGILFILATSSLSVYSILSGASNSKYSLFGALRAVAQTISYEVTIA
IILLSVLLINGSYSLQTLITTQEHILLLPAPIAIIFISTLAETNRAPFDLTEGESELVSGFNVEYAAGPFAL
FFIAEYTNIILINALTTIIFLGPLYYINLPELYSTNFIIEALLLSSTFLWIRASYPRFRYDQLIHLLKNFLP
LTLALCMHISLPIFTAGVPPYI
>ENST00000361390 aa_beforefilter
IPMANLLLLIVPILIAMAFLMLTERKILGYIQLRKGPNVVGPYGLLQPFADAIKLFTKEPLKPATSTITLYI
TAPTLALTIALLL*TPLPIPNPLVNLNLGLLFILATSSLAVYSIL*SG*ASNSNYALIGALRAVAQTISYEV
TLAIILLSTLLISGSFNLSTLITTQEHL*LLLPS*PLAII*FISTLAETNRTPFDLAEGESELVSGFNIEYA
AGPFALFFIAEYTNIIIINTLTTTIFLGTTYDALSPELYTTYFVTKTLLLTSLFL*IRTAYPRFRYDQLIHL
L*KNFLPLTLALLI*YVSIPITISSIPPQT*
>ENST00000361390 aa_afterfilter
IPMANLLLLIVPILIAMAFLMLTERKILGYIQLRKGPNVVGPYGLLQPFADAIKLFTKEPLKPATSTITLYI
TAPTLALTIALLLTPLPIPNPLVNLNLGLLFILATSSLAVYSILSGASNSNYALIGALRAVAQTISYEVTLA
IILLSTLLISGSFNLSTLITTQEHLLLLPSPLAIIFISTLAETNRTPFDLAEGESELVSGFNIEYAAGPFAL
FFIAEYTNIIIINTLTTTIFLGTTYDALSPELYTTYFVTKTLLLTSLFLIRTAYPRFRYDQLIHLLKNFLPL
TLALLIYVSIPITISSIPPQT

Print out the DNA sequences translated back from aa_to_dna function:
>ENSMUST00000082392 aa_to_dna_aln
GTGTTCTTTATTAATATCCTAACACTCCTCGTCCCCATTCTAATCGCCATAGCCTTCCTAACATTAGTAGAA
CGCAAAATCTTAGGGTACATACAACTACGAAAAGGCCCTAACATTGTTGGTCCATACGGCATTTTACAACCA
TTTGCAGACGCCATAAAATTATTTATAAAAGAACCAATACGCCCTTTAACAACCTCTATATCCTTATTTATT
ATTGCACCTACCCTATCACTCACACTAGCATTAAGTCTATGAGTTCCCCTACCAATACCACACCCATTAATT
AATTTAAACCTAGGGATTTTATTTATTTTAGCAACATCTAGCCTATCAGTTTACTCCATTCTATGATCAGGA
TGAGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGTAGCCCAAACAATTTCATATGAAGTA
ACCATAGCTATTATCCTTTTATCAGTTCTATTAATAAATGGATCCTACTCTCTACAAACACTTATTACAACC
CAAGAACACATATGATTACTTCTGCCAGCCTGACCCATAGCCATAATATGATTTATCTCAACCCTAGCAGAA
ACAAACCGGGCCCCCTTCGACCTGACAGAAGGAGAATCAGAATTAGTATCAGGGTTTAACGTAGAATACGCA
GCCGGCCCATTCGCGTTATTCTTTATAGCAGAGTACACTAACATTATTCTAATAAACGCCCTAACAACTATT
ATCTTCCTAGGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACTAACTTCATAATAGAAGCTCTA
CTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATCTT
CTATGAAAAAACTTTCTACCCCTAACACTAGCATTATGTATGTGACATATTTCTTTACCAATTTTT
>ENST00000361390 aa_to_dna_aln
ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCGAA
CGAAAAATTCTAGGCTATATACAACTACGCAAAGGCCCCAACGTTGTAGGCCCCTACGGGCTACTACAACCC
TTCGCTGACGCCATAAAACTCTTCACCAAAGAGCCCCTAAAACCCGCCACATCTACCATCACCCTCTACATC
ACCGCCCCGACCTTAGCTCTCACCATCGCTCTTCTACTATGAACCCCCCTCCCCATACCCAACCCCCTGGTC
AACCTCAACCTAGGCCTCCTATTTATTCTAGCCACCTCTAGCCTAGCCGTTTACTCAATCCTCTGATCAGGG
TGAGCATCAAACTCAAACTACGCCCTGATCGGCGCACTGCGAGCAGTAGCCCAAACAATCTCATATGAAGTC
ACCCTAGCCATCATTCTACTATCAACATTACTAATAAGTGGCTCCTTTAACCTCTCCACCCTTATCACAACA
CAAGAACACCTCTGATTACTCCTGCCATCATGACCCTTGGCCATAATATGATTTATCTCCACACTAGCAGAG
ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAATACGCC
GCAGGCCCCTTCGCCCTATTCTTCATAGCCGAATACACAAACATTATTATAATAAACACCCTCACCACTACA
ATCTTCCTAGGAACAACATATGACGCACTCTCCCCTGAACTCTACACAACATATTTTGTCACCAAGACCCTA
CTT---CTAACCTCCCTGTTCTTATGAATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACAC
CTCCTATGAAAAAACTTCCTACCACTCACCCTAGCATTACTTATATGATATGTCTCCATACCCATT

-------------------- WARNING ---------------------
MSG: There was an error - see error_string for the program output
---------------------------------------------------

------------- 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: main::kaks_calculate calculator.pl:176
STACK: calculator.pl:116

/////////////////////////////////////////////////////////////////////
/////////////////////////////// script //////////////////////////////
/////////////////////////////////////////////////////////////////////
sub kaks_calculate
{
   my %seqs=@_;
   #my %seqs = %$seqs_ref;
   my @prots;

   my $aln_factory = Bio::Tools::Run::Alignment::Clustalw->new
('quiet'=>1);
   
   # process each sequence
   for my $seqid (keys %seqs)
   {
	my $seq = $seqs{$seqid};
	my $protein =$seq->translate();
	my $pseq = $protein->seq();
        print ">$seqid aa_beforefilter \n$pseq\n";
        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;
	print ">$seqid aa_afterfilter \n$pseq\n";
	$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();

   print "\nPrint out the DNA sequences translated back from aa_to_dna
function:\n\n";
   foreach my $s ( $dna_aln->each_seq() ) {
	 print ">".$s->display_id." aa_to_dna_aln\n".$s->seq()."\n";
   }
 
   my $kaks_factory = Bio::Tools::Run::Phylo::PAML::Codeml->new
                   ( -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";
   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 OUT join("\t", $otus[$i]->display_id,
       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";
       }
  }

}


-------------------- WARNING ---------------------
MSG: There was an error - see error_string for the program output
---------------------------------------------------

------------- 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: main::kaks_calculate calculator.pl:176
STACK: calculator.pl:116
----------------------------------------------------------------




On Mon, 2006-07-31 at 11:20 -0400, Ryan Golhar wrote: 
> Hi Xianjun,
> 
> I just did some work on this module including the example.
> 
> >> 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?
> 
> If I'm interpreting your question correctly, the stop codons in your 
> sequence occur in-frame.  This is why it is choking.
> 
> >>So, when translate back from aa_aln to dna_aln, there should be no
> stop codon included. SO, why it can not pass?
> 
> The Ka and Ks statistics are not calculated based on the protein 
> sequence, they are calculated based on the DNA sequence.  The protein 
> sequence is used to provide a alignment for the codons of the DNA 
> sequence.  Checking the protein sequence for * is easier to identify 
> in-frame stop codons than scanning the DNA sequence.
> 
> The two checks for stop codons you mentioned are to check for stop 
> codons within the sequence without worry for the last amino acid.  The

> second part remove the * at the end of the sequence (not the middle).
> 
> If you want to remove the in-frame stop codons, you can, but do so 
> before translating it to protein sequences.
> 
> Ryan
> 
> 
> -----Original Message-----
> From: bioperl-l-bounces at lists.open-bio.org
> [mailto:bioperl-l-bounces at lists.open-bio.org] On Behalf Of Xianjun 
> Dong
> Sent: Monday, July 31, 2006 7:56 AM
> To: bioperl-l at lists.open-bio.org
> Subject: [Bioperl-l] PAML + Codeml problem..
> 
> 
> 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
> /C
> odeml.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
> ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCG
> AA
>
CGAAAAATTCTAGGCTATATACAACTACGCAAAGGCCCCAACGTTGTAGGCCCCTACGGGCTACTACAACCC
>
TTCGCTGACGCCATAAAACTCTTCACCAAAGAGCCCCTAAAACCCGCCACATCTACCATCACCCTCTACATC
>
ACCGCCCCGACCTTAGCTCTCACCATCGCTCTTCTACTATGAACCCCCCTCCCCATACCCAACCCCCTGGTC
>
AACCTCAACCTAGGCCTCCTATTTATTCTAGCCACCTCTAGCCTAGCCGTTTACTCAATCCTCTGATCAGGG
>
TGAGCATCAAACTCAAACTACGCCCTGATCGGCGCACTGCGAGCAGTAGCCCAAACAATCTCATATGAAGTC
>
ACCCTAGCCATCATTCTACTATCAACATTACTAATAAGTGGCTCCTTTAACCTCTCCACCCTTATCACAACA
>
CAAGAACACCTCTGATTACTCCTGCCATCATGACCCTTGGCCATAATATGATTTATCTCCACACTAGCAGAG
>
ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAATACGCC
>
GCAGGCCCCTTCGCCCTATTCTTCATAGCCGAATACACAAACATTATTATAATAAACACCCTCACCACTACA
>
ATCTTCCTAGGAACAACATATGACGCACTCTCCCCTGAACTCTACACAACATATTTTGTCACCAAGACCCTA
>
CTTCTAACCTCCCTGTTCTTATGAATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACCTC
>
CTATGAAAAAACTTCCTACCACTCACCCTAGCATTACTTATATGATATGTCTCCATACCCATTACAATCTCC
> AGCATTCCCCCTCAAACCTAA
> >ENSMUST00000082392
> GTGTTCTTTATTAATATCCTAACACTCCTCGTCCCCATTCTAATCGCCATAGCCTTCCTAACATTAGTAG
> AA
>
CGCAAAATCTTAGGGTACATACAACTACGAAAAGGCCCTAACATTGTTGGTCCATACGGCATTTTACAACCA
>
TTTGCAGACGCCATAAAATTATTTATAAAAGAACCAATACGCCCTTTAACAACCTCTATATCCTTATTTATT
>
ATTGCACCTACCCTATCACTCACACTAGCATTAAGTCTATGAGTTCCCCTACCAATACCACACCCATTAATT
>
AATTTAAACCTAGGGATTTTATTTATTTTAGCAACATCTAGCCTATCAGTTTACTCCATTCTATGATCAGGA
>
TGAGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGTAGCCCAAACAATTTCATATGAAGTA
>
ACCATAGCTATTATCCTTTTATCAGTTCTATTAATAAATGGATCCTACTCTCTACAAACACTTATTACAACC
>
CAAGAACACATATGATTACTTCTGCCAGCCTGACCCATAGCCATAATATGATTTATCTCAACCCTAGCAGAA
>
ACAAACCGGGCCCCCTTCGACCTGACAGAAGGAGAATCAGAATTAGTATCAGGGTTTAACGTAGAATACGCA
>
GCCGGCCCATTCGCGTTATTCTTTATAGCAGAGTACACTAACATTATTCTAATAAACGCCCTAACAACTATT
>
ATCTTCCTAGGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACTAACTTCATAATAGAAGCTCTA
>
CTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATCTT
>
CTATGAAAAAACTTTCTACCCCTAACACTAGCATTATGTATGTGACATATTTCTTTACCAATTTTTACAGCG
> GGAGTACCACCATACATATAG
> 
> Sure, I checked it. There is some stop codon in it. If I replace it 
> with non-stop codon, it works.
> 
> For example,
> >ENST00000361390
> ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCcaaTCGCAATGGCATTCCcaaTGCTTACCG
> AA
>
CGAAAAATTCcaaGCTATATACAACTACGCAAAGGCCCCAACGTTGcaaGCCCCTACGGGCTACTACAACCC
>
TTCGCcaaCGCCAcaaAACTCTTCACCAAAGAGCCCCcaaAACCCGCCACATCTACCATCACCCTCTACATC
>
ACCGCCCCGACCTcaaCTCTCACCATCGCTCTTCTACTAcaaACCCCCCTCCCCATACCCAACCCCCTGGTC
>
AACCTCAACCcaaGCCTCCTATTTATTCcaaCCACCTCcaaCCcaaCCGTTTACTCAATCCTCcaaTCAGGG
>
caaGCATCAAACTCAAACTACGCCCcaaTCGGCGCACTGCGAGCAGcaaCCCAAACAATCTCATAcaaAGTC
>
ACCCcaaCCATCATTCTACTATCAACATTACcaacaaGTGGCTCCTTcaaCCTCTCCACCCTTATCACAACA
>
CAAGAACACCTCcaaTTACTCCTGCCATCAcaaCCCTTGGCCAcaaTAcaaTTTATCTCCACACcaaCAGAG
>
ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACcaaTCTCAGGCTTCAACATCGAATACGCC
>
GCAGGCCCCTTCGCCCTATTCTTCAcaaCCGAATACACAAACATTATTAcaacaaACACCCTCACCACTACA
>
ATCTTCCcaaGAACAACATAcaaCGCACTCTCCCCcaaACTCTACACAACATATTTTGTCACCAAGACCCTA
>
CTTCcaaCCTCCCTGTTCTTAcaaATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACCTC
>
CTAcaaAAAAACTTCCTACCACTCACCCcaaCATTACTTATAcaaTATGTCTCCATACCCATTACAATCTCC
> AGCATTCCCCCTCAAACCcaa
> >ENSMUST00000082392
> GTGTTCTTTATcaaTATCCcaaCACTCCTCGTCCCCATTCcaaTCGCCAcaaCCTTCCcaaCATcaacaa
> AA
>
CGCAAAATCTcaaGGTACATACAACTACGAAAAGGCCCcaaCATTGTTGGTCCATACGGCATTTTACAACCA
>
TTTGCAGACGCCAcaaAATTATTTAcaaAAGAACCAATACGCCCTTcaaCAACCTCTATATCCTTATTTATT
>
ATTGCACCTACCCTATCACTCACACcaaCATcaaGTCTAcaaGTTCCCCTACCAATACCACACCCATcaaTc
>
aaTTcaaACCcaaGGATTTTATTTATTTcaaCAACATCcaaCCTATCAGTTTACTCCATTCTAcaaTCAGGA
>
caaGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGcaaCCCAAACAATTTCATAcaaAGca
>
aCCAcaaCTATTATCCTTTTATCAGTTCTATcaacaaATGGATCCTACTCTCTACAAACACTTATTACAACC
>
CAAGAACACATAcaaTTACTTCTGCCAGCCcaaCCCAcaaCCAcaaTAcaaTTTATCTCAACCCcaaCAGAA
>
ACAAACCGGGCCCCCTTCGACCcaaCAGAAGGAGAATCAGAATcaaTATCAGGGTTcaaCGcaaAATACGCA
>
GCCGGCCCATTCGCGTTATTCTTTAcaaCAGAGTACACcaaCATTATTCcaacaaACGCCCcaaCAACTATT
>
ATCTTCCcaaGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACcaaCTTCAcaacaaAAGCTCTA
>
CTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATCTT
>
CTAcaaAAAAACTTTCTACCCCcaaCACcaaCATTATGTATGcaaCATATTTCTTTACCAATTTTTACAGCG
> GGAGTACCACCATACATAcaa
> 
> 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";
>   }
> }
> 



More information about the Bioperl-l mailing list