[Bioperl-l] PAML::Codeml outputs unstable value, why?
Albert Vilella
avilella at gmail.com
Tue May 29 13:45:28 UTC 2007
On 5/29/07, Dong Xianjun <Xianjun.Dong at bccs.uib.no> wrote:
>
> Thanks for information, Albert.
>
> But still in two questions:
> Albert Vilella wrote:
>
> 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.
>
> 1. How to set the initial value in order to get a reasonable estimation?
> Do you have some experience for that?
>
People usually change the initial omega in the conf. For example, 3 runs
with 0.001, 1 and 5.
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.
>
> 2. Is there a recommend way to test the significance if the results are
> different? For example, in my case, dS could range from 10.1852 to 14.9961for the four runtime. If that means the model is not robust(how to check
> this?), should I change to use another model?
>
I would prefer PAML's author to answer this question :)
How could YN00 reach stable result? (Is it because YN00 does not require
> initial value for optimization?) Why could YN00 produce so different result
> from Codeml? (for YN00, dS=2.1300 with SE=1.2272; for Codeml, dS=
> 10.1852-14.9961)
>
I think Yn00 is less prone to give different local maxima than some codeml
models, but then, codeml is better in giving true positives in cases where
yn00 will give false negatives...
Maybe PAML's author can give a more specific answer for your data at:
> http://www.rannala.org/gsf/viewforum.php?f=1
>
>
> Actually I already post my question in the author's forum. Let's wait and
> see.
>
Yes, I would wait for his answers, which should be way more reliable than
mine :)
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.ctlfile 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
> >
> >
>
> --
> ---------------------------
> 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
>
>
More information about the Bioperl-l
mailing list