[Bioperl-l] PAML + Codeml problem..
Ryan Golhar
golharam at umdnj.edu
Tue Aug 15 00:24:42 UTC 2006
Hi Xianjun,
1. It looks like a number of people have already responded to this one.
My two cents - I think if you created a custom codon table and used
that, it would be okay, but that doesn't solve your underlying problem.
PAML will still complain because it does not know selenocysteine. It
will see the generic stop codon and halt because of that.
Unfortunately, there is nothing bioperl can do about it as PAML is an
external package to bioperl.
2. You can remove the stop codons from the CDS sequence. Then, when
you run PAML it should be okay. The filter procedure in the code is not
designed to remove the in-frame stop codons, its is only designed to
warn you about its presence. It will only remove the trailing stop
codon as PAML seems to be okay with this.
3. BioPerl (as I've learned) has Ka/Ks calculation routines built in
the DNAStatistics module. You could use those. The PBL tool I am using
was written by a colleuege whom I will need to get permission from to
send to you first.
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, August 14, 2006 11:57 AM
To: golharam at umdnj.edu
Cc: bioperl-l at lists.open-bio.org
Subject: Re: [Bioperl-l] PAML + Codeml problem..
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
> ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCG
> AA
>
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
> IPMANLLLLIVPILIAMAFLMLTERKILGYIQLRKGPNVVGPYGLLQPFADAIKLFTKEPLKPATSTITL
> YI
>
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
> ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCG
> AA
>
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
> 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
>
> Calculate the Ka/Ks for ENSG00000198888 : ENSMUSG00000064341 ...
> >ENSMUST00000082392 aa_beforefilter
> VFFINILTLLVPILIAIAFLTLVERKILGYIQLRKGPNIVGPYGILQPFADAIKLFIKEPIRPLTTSISL
> FI
>
IAPTLSLTLALSL*VPLPIPHPLINLNLGILFILATSSLSVYSIL*SG*ASNSKYSLFGALRAVAQTISYEV
>
TIAIILLSVLLINGSYSLQTLITTQEHI*LLLPA*PIAII*FISTLAETNRAPFDLTEGESELVSGFNVEYA
>
AGPFALFFIAEYTNIILINALTTIIFLGPLYYINLPELYSTNFIIEALLLSSTFLWIRASYPRFRYDQLIHL
> L*KNFLPLTLALCM*HISLPIFTAGVPPYI*
> >ENSMUST00000082392 aa_afterfilter
> VFFINILTLLVPILIAIAFLTLVERKILGYIQLRKGPNIVGPYGILQPFADAIKLFIKEPIRPLTTSISL
> FI
>
IAPTLSLTLALSLVPLPIPHPLINLNLGILFILATSSLSVYSILSGASNSKYSLFGALRAVAQTISYEVTIA
>
IILLSVLLINGSYSLQTLITTQEHILLLPAPIAIIFISTLAETNRAPFDLTEGESELVSGFNVEYAAGPFAL
>
FFIAEYTNIILINALTTIIFLGPLYYINLPELYSTNFIIEALLLSSTFLWIRASYPRFRYDQLIHLLKNFLP
> LTLALCMHISLPIFTAGVPPYI
> >ENST00000361390 aa_beforefilter
> IPMANLLLLIVPILIAMAFLMLTERKILGYIQLRKGPNVVGPYGLLQPFADAIKLFTKEPLKPATSTITL
> YI
>
TAPTLALTIALLL*TPLPIPNPLVNLNLGLLFILATSSLAVYSIL*SG*ASNSNYALIGALRAVAQTISYEV
>
TLAIILLSTLLISGSFNLSTLITTQEHL*LLLPS*PLAII*FISTLAETNRTPFDLAEGESELVSGFNIEYA
>
AGPFALFFIAEYTNIIIINTLTTTIFLGTTYDALSPELYTTYFVTKTLLLTSLFL*IRTAYPRFRYDQLIHL
> L*KNFLPLTLALLI*YVSIPITISSIPPQT*
> >ENST00000361390 aa_afterfilter
> IPMANLLLLIVPILIAMAFLMLTERKILGYIQLRKGPNVVGPYGLLQPFADAIKLFTKEPLKPATSTITL
> YI
>
TAPTLALTIALLLTPLPIPNPLVNLNLGLLFILATSSLAVYSILSGASNSNYALIGALRAVAQTISYEVTLA
>
IILLSTLLISGSFNLSTLITTQEHLLLLPSPLAIIFISTLAETNRTPFDLAEGESELVSGFNIEYAAGPFAL
>
FFIAEYTNIIIINTLTTTIFLGTTYDALSPELYTTYFVTKTLLLTSLFLIRTAYPRFRYDQLIHLLKNFLPL
> TLALLIYVSIPITISSIPPQT
>
> Print out the DNA sequences translated back from aa_to_dna function:
> >ENSMUST00000082392 aa_to_dna_aln
> GTGTTCTTTATTAATATCCTAACACTCCTCGTCCCCATTCTAATCGCCATAGCCTTCCTAACATTAGTAG
> AA
>
CGCAAAATCTTAGGGTACATACAACTACGAAAAGGCCCTAACATTGTTGGTCCATACGGCATTTTACAACCA
>
TTTGCAGACGCCATAAAATTATTTATAAAAGAACCAATACGCCCTTTAACAACCTCTATATCCTTATTTATT
>
ATTGCACCTACCCTATCACTCACACTAGCATTAAGTCTATGAGTTCCCCTACCAATACCACACCCATTAATT
>
AATTTAAACCTAGGGATTTTATTTATTTTAGCAACATCTAGCCTATCAGTTTACTCCATTCTATGATCAGGA
>
TGAGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGTAGCCCAAACAATTTCATATGAAGTA
>
ACCATAGCTATTATCCTTTTATCAGTTCTATTAATAAATGGATCCTACTCTCTACAAACACTTATTACAACC
>
CAAGAACACATATGATTACTTCTGCCAGCCTGACCCATAGCCATAATATGATTTATCTCAACCCTAGCAGAA
>
ACAAACCGGGCCCCCTTCGACCTGACAGAAGGAGAATCAGAATTAGTATCAGGGTTTAACGTAGAATACGCA
>
GCCGGCCCATTCGCGTTATTCTTTATAGCAGAGTACACTAACATTATTCTAATAAACGCCCTAACAACTATT
>
ATCTTCCTAGGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACTAACTTCATAATAGAAGCTCTA
>
CTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATCTT
> CTATGAAAAAACTTTCTACCCCTAACACTAGCATTATGTATGTGACATATTTCTTTACCAATTTTT
> >ENST00000361390 aa_to_dna_aln
> ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTACCG
> AA
>
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
> > ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCTAATCGCAATGGCATTCCTAATGCTTAC
> > CG
> > AA
> >
> CGAAAAATTCTAGGCTATATACAACTACGCAAAGGCCCCAACGTTGTAGGCCCCTACGGGCTACTACAAC
> CC
> >
> TTCGCTGACGCCATAAAACTCTTCACCAAAGAGCCCCTAAAACCCGCCACATCTACCATCACCCTCTACA
> TC
> >
> ACCGCCCCGACCTTAGCTCTCACCATCGCTCTTCTACTATGAACCCCCCTCCCCATACCCAACCCCCTGG
> TC
> >
> AACCTCAACCTAGGCCTCCTATTTATTCTAGCCACCTCTAGCCTAGCCGTTTACTCAATCCTCTGATCAG
> GG
> >
> TGAGCATCAAACTCAAACTACGCCCTGATCGGCGCACTGCGAGCAGTAGCCCAAACAATCTCATATGAAG
> TC
> >
> ACCCTAGCCATCATTCTACTATCAACATTACTAATAAGTGGCTCCTTTAACCTCTCCACCCTTATCACAA
> CA
> >
> CAAGAACACCTCTGATTACTCCTGCCATCATGACCCTTGGCCATAATATGATTTATCTCCACACTAGCAG
> AG
> >
> ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACTAGTCTCAGGCTTCAACATCGAATACG
> CC
> >
> GCAGGCCCCTTCGCCCTATTCTTCATAGCCGAATACACAAACATTATTATAATAAACACCCTCACCACTA
> CA
> >
> ATCTTCCTAGGAACAACATATGACGCACTCTCCCCTGAACTCTACACAACATATTTTGTCACCAAGACCC
> TA
> >
> CTTCTAACCTCCCTGTTCTTATGAATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACC
> TC
> >
> CTATGAAAAAACTTCCTACCACTCACCCTAGCATTACTTATATGATATGTCTCCATACCCATTACAATCT
> CC
> > AGCATTCCCCCTCAAACCTAA
> > >ENSMUST00000082392
> > GTGTTCTTTATTAATATCCTAACACTCCTCGTCCCCATTCTAATCGCCATAGCCTTCCTAACATTAGT
> > AG
> > AA
> >
> CGCAAAATCTTAGGGTACATACAACTACGAAAAGGCCCTAACATTGTTGGTCCATACGGCATTTTACAAC
> CA
> >
> TTTGCAGACGCCATAAAATTATTTATAAAAGAACCAATACGCCCTTTAACAACCTCTATATCCTTATTTA
> TT
> >
> ATTGCACCTACCCTATCACTCACACTAGCATTAAGTCTATGAGTTCCCCTACCAATACCACACCCATTAA
> TT
> >
> AATTTAAACCTAGGGATTTTATTTATTTTAGCAACATCTAGCCTATCAGTTTACTCCATTCTATGATCAG
> GA
> >
> TGAGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGTAGCCCAAACAATTTCATATGAAG
> TA
> >
> ACCATAGCTATTATCCTTTTATCAGTTCTATTAATAAATGGATCCTACTCTCTACAAACACTTATTACAA
> CC
> >
> CAAGAACACATATGATTACTTCTGCCAGCCTGACCCATAGCCATAATATGATTTATCTCAACCCTAGCAG
> AA
> >
> ACAAACCGGGCCCCCTTCGACCTGACAGAAGGAGAATCAGAATTAGTATCAGGGTTTAACGTAGAATACG
> CA
> >
> GCCGGCCCATTCGCGTTATTCTTTATAGCAGAGTACACTAACATTATTCTAATAAACGCCCTAACAACTA
> TT
> >
> ATCTTCCTAGGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACTAACTTCATAATAGAAGCTC
> TA
> >
> CTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATC
> TT
> >
> CTATGAAAAAACTTTCTACCCCTAACACTAGCATTATGTATGTGACATATTTCTTTACCAATTTTTACAG
> CG
> > 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
> > ATACCCATGGCCAACCTCCTACTCCTCATTGTACCCATTCcaaTCGCAATGGCATTCCcaaTGCTTAC
> > CG
> > AA
> >
> CGAAAAATTCcaaGCTATATACAACTACGCAAAGGCCCCAACGTTGcaaGCCCCTACGGGCTACTACAAC
> CC
> >
> TTCGCcaaCGCCAcaaAACTCTTCACCAAAGAGCCCCcaaAACCCGCCACATCTACCATCACCCTCTACA
> TC
> >
> ACCGCCCCGACCTcaaCTCTCACCATCGCTCTTCTACTAcaaACCCCCCTCCCCATACCCAACCCCCTGG
> TC
> >
> AACCTCAACCcaaGCCTCCTATTTATTCcaaCCACCTCcaaCCcaaCCGTTTACTCAATCCTCcaaTCAG
> GG
> >
> caaGCATCAAACTCAAACTACGCCCcaaTCGGCGCACTGCGAGCAGcaaCCCAAACAATCTCATAcaaAG
> TC
> >
> ACCCcaaCCATCATTCTACTATCAACATTACcaacaaGTGGCTCCTTcaaCCTCTCCACCCTTATCACAA
> CA
> >
> CAAGAACACCTCcaaTTACTCCTGCCATCAcaaCCCTTGGCCAcaaTAcaaTTTATCTCCACACcaaCAG
> AG
> >
> ACCAACCGAACCCCCTTCGACCTTGCCGAAGGGGAGTCCGAACcaaTCTCAGGCTTCAACATCGAATACG
> CC
> >
> GCAGGCCCCTTCGCCCTATTCTTCAcaaCCGAATACACAAACATTATTAcaacaaACACCCTCACCACTA
> CA
> >
> ATCTTCCcaaGAACAACATAcaaCGCACTCTCCCCcaaACTCTACACAACATATTTTGTCACCAAGACCC
> TA
> >
> CTTCcaaCCTCCCTGTTCTTAcaaATTCGAACAGCATACCCCCGATTCCGCTACGACCAACTCATACACC
> TC
> >
> CTAcaaAAAAACTTCCTACCACTCACCCcaaCATTACTTATAcaaTATGTCTCCATACCCATTACAATCT
> CC
> > AGCATTCCCCCTCAAACCcaa
> > >ENSMUST00000082392
> > GTGTTCTTTATcaaTATCCcaaCACTCCTCGTCCCCATTCcaaTCGCCAcaaCCTTCCcaaCATcaac
> > aa
> > AA
> >
> CGCAAAATCTcaaGGTACATACAACTACGAAAAGGCCCcaaCATTGTTGGTCCATACGGCATTTTACAAC
> CA
> >
> TTTGCAGACGCCAcaaAATTATTTAcaaAAGAACCAATACGCCCTTcaaCAACCTCTATATCCTTATTTA
> TT
> >
> ATTGCACCTACCCTATCACTCACACcaaCATcaaGTCTAcaaGTTCCCCTACCAATACCACACCCATcaa
> Tc
> >
> aaTTcaaACCcaaGGATTTTATTTATTTcaaCAACATCcaaCCTATCAGTTTACTCCATTCTAcaaTCAG
> GA
> >
> caaGCCTCAAACTCCAAATACTCACTATTCGGAGCTTTACGAGCCGcaaCCCAAACAATTTCATAcaaAG
> ca
> >
> aCCAcaaCTATTATCCTTTTATCAGTTCTATcaacaaATGGATCCTACTCTCTACAAACACTTATTACAA
> CC
> >
> CAAGAACACATAcaaTTACTTCTGCCAGCCcaaCCCAcaaCCAcaaTAcaaTTTATCTCAACCCcaaCAG
> AA
> >
> ACAAACCGGGCCCCCTTCGACCcaaCAGAAGGAGAATCAGAATcaaTATCAGGGTTcaaCGcaaAATACG
> CA
> >
> GCCGGCCCATTCGCGTTATTCTTTAcaaCAGAGTACACcaaCATTATTCcaacaaACGCCCcaaCAACTA
> TT
> >
> ATCTTCCcaaGACCCCTATACTATATCAATTTACCAGAACTCTACTCAACcaaCTTCAcaacaaAAGCTC
> TA
> >
> CTACTATCATCAACATTCCTATGGATCCGAGCATCTTATCCACGCTTCCGTTACGATCAACTTATACATC
> TT
> >
> CTAcaaAAAAACTTTCTACCCCcaaCACcaaCATTATGTATGcaaCATATTTCTTTACCAATTTTTACAG
> CG
> > 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";
> > }
> > }
> >
>
_______________________________________________
Bioperl-l mailing list
Bioperl-l at lists.open-bio.org
http://lists.open-bio.org/mailman/listinfo/bioperl-l
More information about the Bioperl-l
mailing list