[Bioperl-l] PAML::Codeml outputs unstable value, why?
Albert Vilella
avilella at gmail.com
Tue May 29 13:02:44 UTC 2007
codeml in PAML can give different results in cases where the optimization
reaches different local maxima depending on the different starting points of
each run (seed values). So, at least for some methods and options, this
instability is inherent to the underlying algorithm.
Even more, for some methods and options, it is even recommended in PAML
documentation to run the same data more than once, to see if the results are
the same, which would be a good indication that the model is robust given
the data.
Maybe PAML's author can give a more specific answer for your data at:
http://www.rannala.org/gsf/viewforum.php?f=1
Cheers,
Albert.
On 5/29/07, Dong Xianjun <Xianjun.Dong at bccs.uib.no> wrote:
>
> HI, dear all, //sorry for duplicated msg for *Jason* and *Neil*
>
> I'm bothering by two problems when I use PAML module to calculate Ka/Ks
> for my sequences. Could you help me?
>
> 1. Codeml could produce different Ka/Ks value if I run it twice. I check
> it both in command line and in Perl wrapper of
> Bio::Tools::Run::Phylo::PAML::Codeml;
>
> The input sequences are:
> >seq1
> TCTCTCTGGCCCAAAATCCGGGTTCCATTAAAAGTTGTGAGGACTGCTGAAAACAAGTTAAGTAACCGTTTCTTCCCTTATGATGAAATCGAGACAGAAGCTGTTCTGGCCATTGATGATGATATCATTATGCTGACCTCTGACGAGCTGCAATTTGGTTATGAG
> >seq2
>
> TCACTGTGGCCCAAAGTCGCAGTGCCTCTTAAAGTGGTCCGCACCAAAGAAAACAAGCTCAGCAATCGATTCTTTCCGTTTGATGAGATCGAGACAGAAGCTGTCCTGGCCATTGACGATGACATCATCATGTTAACCTCAGATGAGCTACAGTTTGGATATGAG
>
> For command-line program, I used Codeml in PAML3.14, with specifications
> in codeml.ctl (runmode = -2, seqtype = 1). I tried to run the program four
> times. The output are like below (from the output file). We could see that
> they are different from each other. they should be same or slightly
> different. Right? But they are NOT. Weird!
>
> ----------------------------------------------------------------------------------------------------------------------------------
> t=11.5447 S= 42.4 N= 122.6 dN/dS= 0.0035 dN= 0.0522 dS=14.8339
> t= 9.4132 S= 41.8 N= 123.2 dN/dS= 0.0041 dN= 0.0507 dS=12.2349
> t=11.6305 S= 42.2 N= 122.8 dN/dS= 0.0034 dN= 0.0510 dS=14.9961
> t= 7.7879 S= 41.4 N= 123.6 dN/dS= 0.0050 dN= 0.0505 dS=10.1852
>
> ----------------------------------------------------------------------------------------------------------------------------------
> I found the same problem when I use the Perl Wrapper of
> Bio::Tools::Run::Phylo::PAML::Codeml; (I attached my Perl script here,
> similar to the one in BioPerl HOWTO).
>
> 2. Another strange thing is, if I switch to use program YN00 in the
> package of PAML, the output are stable. However, it's much different from
> Codeml. (see below)
>
> ----------------------------------------------------------------------------------------------------------------------------------
> seq. seq. S N t kappa omega dN +- SE
> dS +- SE
> 2 1 40.4 124.6 1.7452 1.3163 0.0378 0.0804 +- 0.0265
> 2.1300 +- 1.2272
>
> ----------------------------------------------------------------------------------------------------------------------------------
> Why like this? Which one I should believe?
>
>
> Is there any guy who would kindly help me to run the perl script (twice to
> check whether they are different)? or help to run the codeml in command
> line?
> I don't know whether there is anyone noticed this before, or because of
> the wrong version of PAML.
>
> Regards,
>
> Xianjun
>
>
>
> Himanshu Ardawatia wrote:
>
> #!/usr/bin/perl
>
> use strict;
> use warnings;
>
>
> use Bio::Tools::Run::Phylo::PAML::Codeml;
> use Bio::Tools::Run::Alignment::Clustalw;
>
> # for projecting alignments from protein to R/DNA space
> use Bio::Align::Utilities qw(aa_to_dna_aln);
>
> # for input of the sequence data
> use Bio::SeqIO;
> use Bio::AlignIO;
>
> my $aln_factory = new Bio::Tools::Run::Alignment::Clustalw();
>
> #my $seqdata = 'chuck.fa';
> my $seqdata = 'xianjun.fa ';
>
> my $seqIO = new Bio::SeqIO(-file => $seqdata,
> -format => 'fasta');
> my %seqs;
> my @prots;
>
> my $output;
> # process each sequence
> while( my $seq = $seqIO->next_seq ) {
> $seqs{$seq->display_id} = $seq;
> # translate them into protein
> my $protein = $seq->translate();
> my $pseq = $protein->seq();
> if( $pseq =~ /\*/ &&
> $pseq !~ /\*$/ ) {
> warn("provided a cDNA sequence with a stop codon, PAML will choke!");
> exit(0);
> }
> # Tcoffee can't handle '*' even if it is trailing
> $pseq =~ s/\*//g;
> $protein->seq($pseq);
> push @prots, $protein;
> }
>
> if( @prots < 2 ) {
> warn("Need at least 2 cDNA sequences to proceed");
> exit(0);
> }
>
> open(OUT, ">align_output.txt") ||
> die("cannot open output $output for writing");
> # Align the sequences with clustalw
>
> my $aa_aln = $aln_factory->align(\@prots);
>
> # project the protein alignment back to cDNA coordinates
> my $dna_aln = &aa_to_dna_aln($aa_aln, \%seqs);
>
> my @each = $dna_aln->each_seq();
>
> my $kaks_factory = new Bio::Tools::Run::Phylo::PAML::Codeml
> ( -params => { 'runmode' => -2,
> 'seqtype' => 1,
> 'model' => 1,
> }
> );
>
> # set the alignment object
> $kaks_factory->alignment($dna_aln);
>
> # run the KaKs analysis
> my ($rc,$parser) = $kaks_factory->run();
> my $result = $parser->next_result;
> my $MLmatrix = $result->get_MLmatrix();
>
> my @otus = $result->get_seqs();
> # this gives us a mapping from the PAML order of sequences back to
> # the input order (since names get truncated)
> my @pos = map {
> my $c= 1;
> foreach my $s ( @each ) {
> last if( $s->display_id eq $_->display_id );
> $c++;
> }
> $c;
> } @otus;
>
> print OUT join("\t", qw(SEQ1 SEQ2 Ka Ks Ka/Ks PROT_PERCENTID
> CDNA_PERCENTID)), "\n";
> for( my $i = 0; $i < (scalar @otus -1) ; $i++) {
> for( my $j = $i+1; $j < (scalar @otus); $j++ ) {
> my $sub_aa_aln = $aa_aln->select_noncont($pos[$i],$pos[$j]);
> my $sub_dna_aln = $dna_aln->select_noncont($pos[$i],$pos[$j]);
> print OUT join("\t",
> $otus[$i]->display_id,
> $otus[$j]->display_id,$MLmatrix->[$i]->[$j]->{'dN'},
> $MLmatrix->[$i]->[$j]->{'dS'},
> $MLmatrix->[$i]->[$j]->{'omega'},
> sprintf("%.2f",$sub_aa_aln->percentage_identity),
> sprintf("%.2f",$sub_dna_aln->percentage_identity),
> ), "\n";
> }
> }
>
>
> On 5/29/07, Himanshu Ardawatia <himanshu.ardawatia at bccs.uib.no > wrote:
> >
> > Hi Xianjun,
> >
> > I recognize this script. But it was a bit cumbersom to use this as many
> > things are done in the script (like multiple alignment, aa to dna alignment
> > and ka/ks calculation) so one does not have real control on these different
> > aspect.
> > I do not remeber getting different Ka/Ks in different runs though. But I
> > remeber that one I ran the script with different versions of clustalw and it
> > REALLY gave different results !! So please make sure if the clustalw
> > versions are the same in all your runs. Best is to use the latest version.
> >
> > Finally I wrote my simple script which would generate a codeml.ctl file
> > for each set of sequences and run the codeml based on that and then more on.
> > Disadvantage of this can be that some files keep getting over-written (like
> > the one which have their names hard-coded in codeml program) and if one
> > needs those files as well then one needs to run the codeml cycles for each
> > set of sequences in different directories.
> >
> > One advantage of this kind of script is that you can use whichever
> > alignment program you want to use and so on....But then its also extra steps
> > of yourself doing multiple alignment and aa to dna alignment etc....
> >
> > Does it make sense? If you still get different outputs with same version
> > of clustalw then I can sit with you and look at things together. Or else try
> > the script method which I mentioned.
> >
> > Cheers and Fu
> > Himanshu
> > \\
> > On 5/28/07, Dong Xianjun < Xianjun.Dong at bccs.uib.no> wrote:
> > >
> > > HI, Himanshu
> > >
> > > I am sure you did some work in Ka/Ks calculation. Here I have a
> > > question
> > > bothering me; the output for Bio::Tools::Run::Phylo::PAML::Codeml is
> > > not
> > > stable(different for each runtime), and also different from the output
> > >
> > > with modeul of Bio::Tools::Run::Phylo::PAML::Yn00.
> > >
> > > Here I attached the script. Could you help to have a look and try to
> > > run
> > > the script? How is your way to calculate the Kaks ratio?
> > >
> > > Thanks
> > >
> > > --
> > > ---------------------------
> > > Sterding (Xianjun) Dong
> > > PhD student, Boris Lenhard's group
> > > Bergen Center of Computational Science
> > > Bergen University, Norway
> > > Mobile: 0047-47361688
> > > Telephone: 0047-55276381
> > > Skype: xianjun.dong
> > >
> > >
> > >
> >
>
> --
> ---------------------------
> Sterding (Xianjun) Dong
> PhD student, Boris Lenhard's group
> Bergen Center of Computational Science
> Bergen University, Norway
> Mobile: 0047-47361688
> Telephone: 0047-55276381
> Skype: xianjun.dong
>
>
> _______________________________________________
> 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