[Bioperl-l] PAML + Codeml problem..

Xianjun Dong xianjun.dong at bccs.uib.no
Mon Aug 14 15:57:14 UTC 2006


Hi, Ryan and all other helpers,

I finally could run my script and solved the problem of codonTable. (I
checked the DNA type -- mtDNA or nucleotide DNA -- first before I call
translate). Thanks a lot for your help. 

But I still have some questions:
1. For the case which in-frame stop codon codes for selenocysteine('U'),
like the transcript ENSMUST00000094469, it should be translated into
'U', not '*' since the IUPAC/IUBMB has officially recommended it. But
when I use the codontable_id=1(generic codon table), it still was '*'.
Is it because the package(Bio::Tools::CodonTable) is not so updated as
the IUPAC rules?

2. Ryan, I still want to confirm one point for your sample code: Can I
just directly remove the in-frame stop codons (both in the middle and in
the tail) from the CDS sequence, and then get dna_aln by Clustalw, and
then invoke run() on the Codeml package? I don't think the filter
procedure in the sample code works very well.

3. What's more, there are two ways to get Ka/Ks through the PAML
package:
my $yn = new Bio::Tools::Run::Phylo::PAML::Yn00();  
and
my $kaks_factory = Bio::Tools::Run::Phylo::PAML::Codeml->new
                   ( -params => { 'runmode' => -2,
                                  'seqtype' => 1,
                                } );

I checked both PODs for this two modules. The default setting for Yn00()
should be same as the above Codeml setting. But the Ks output for the
same sequences is much different.

For example, here is the output for the sequences below:

[xianjund at lauvtre kaks]$ perl paml.pl seq.fa
Yn00: 
Ka = 0.6267 Ks = 0.9160 Ka/Ks = 0.6841
Codeml:
SEQ1    SEQ2    Ka      Ks      Ka/Ks   PROT_PERCENTID  CDNA_PERCENTID
ENSMUST00000094469      ENST00000361918 0.7419  1.6483  0.4501  47.62
55.16

Sequences are here:
>ENSMUST00000094469
ATGAGCATCCTACTGTCGCCGCCGTCGCTGCTGCTGCTTCTTGCAGCCCTTGTGGCTCCA
GCCACCTCCACCACCAACTACCGACCGGATTGGAACCGTCTTCGAGGCCTGGCCAGGGGG
CGGGTGGAGACCTGTGGAGGACAGTTGAATCGCCTAAAGGAGGTGAAGGCCTTTGTCAAA
GAAGCTCAGGTGCCCCCCGAGTACCTGTGGGCGCCCGCTAAGCCCCCCGAGGAAGCTTCA
GAACACGACTGGCTGTGA
>ENST00000361918
ATGAGCCTCCTGTTGCCTCCGCTGGCGCTGCTGCTGCTTCTCGCGGCGCTTGTGGCCCCA
GAGCTCGTGCTGCTGGGCCGCCGCTACGAGGAACTAGAGCGCATCCCACTCAGTGAAATG
ACCCGCGAAGAGATCAATGCGCTAGTGCAGGAGCTCGGCTTCTACCGCAAGGCGGCGCCC
GACGCGCAGGTGCCCCCCGAGTACGTGTGGGCGCCCGCGAAGCCCCCAGAGGAAACTTCG
GACCACGCTGACCTGTAG

4. BTW, could you share your method PBL with me? I want to learn more on
how to overcome the overestimate synonymous rates cases.

Thanks!

-Xianjun



On Thu, 2006-08-10 at 14:53 -0400, Ryan Golhar wrote':
> 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