<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html;charset=ISO-8859-1" http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000000">
HI, dear all, //sorry for duplicated msg for <i>Jason</i> and <i>Neil</i><br>
<br>
I'm bothering by two problems when I use PAML module to calculate Ka/Ks
for my sequences. Could you help me?<br>
<br>
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;<br>
<br>
The input sequences are:<br>
>seq1
TCTCTCTGGCCCAAAATCCGGGTTCCATTAAAAGTTGTGAGGACTGCTGAAAACAAGTTAAGTAACCGTTTCTTCCCTTATGATGAAATCGAGACAGAAGCTGTTCTGGCCATTGATGATGATATCATTATGCTGACCTCTGACGAGCTGCAATTTGGTTATGAG<br>
>seq2<br>
TCACTGTGGCCCAAAGTCGCAGTGCCTCTTAAAGTGGTCCGCACCAAAGAAAACAAGCTCAGCAATCGATTCTTTCCGTTTGATGAGATCGAGACAGAAGCTGTCCTGGCCATTGACGATGACATCATCATGTTAACCTCAGATGAGCTACAGTTTGGATATGAG<br>
<br>
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!<br>
----------------------------------------------------------------------------------------------------------------------------------<br>
<tt>t=11.5447 S= 42.4 N= 122.6 dN/dS= 0.0035 dN= 0.0522
dS=14.8339<br>
t= 9.4132 S= 41.8 N= 123.2 dN/dS= 0.0041 dN= 0.0507 dS=12.2349<br>
t=11.6305 S= 42.2 N= 122.8 dN/dS= 0.0034 dN= 0.0510 dS=14.9961<br>
t= 7.7879 S= 41.4 N= 123.6 dN/dS= 0.0050 dN= 0.0505 dS=10.1852<br>
</tt>----------------------------------------------------------------------------------------------------------------------------------<br>
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).<br>
<br>
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) <br>
----------------------------------------------------------------------------------------------------------------------------------<br>
<tt>seq. seq. S N t kappa omega dN +- SE
dS +- SE<br>
2 1 40.4 124.6 1.7452 1.3163 0.0378 0.0804 +- 0.0265
2.1300 +- 1.2272</tt><br>
----------------------------------------------------------------------------------------------------------------------------------<br>
Why like this? Which one I should believe?<br>
<br>
<br>
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?<br>
I don't know whether there is anyone noticed this before, or because of
the wrong version of PAML.<br>
<br>
Regards,<br>
<br>
Xianjun<br>
<br>
<br>
<br>
Himanshu Ardawatia wrote:
<blockquote
cite="mid62d36e2b0705290158h1c85362cp824778ca5ecc8645@mail.gmail.com"
type="cite">#!/usr/bin/perl<br>
<br>
use strict;<br>
use warnings;<br>
<br>
<br>
use Bio::Tools::Run::Phylo::PAML::Codeml;<br>
use Bio::Tools::Run::Alignment::Clustalw;<br>
<br>
# for projecting alignments from protein to R/DNA space<br>
use Bio::Align::Utilities qw(aa_to_dna_aln); <br>
<br>
# for input of the sequence data<br>
use Bio::SeqIO;<br>
use Bio::AlignIO;<br>
<br>
my $aln_factory = new Bio::Tools::Run::Alignment::Clustalw();<br>
<br>
#my $seqdata = 'chuck.fa';<br>
my $seqdata = 'xianjun.fa
';<br>
<br>
my $seqIO = new Bio::SeqIO(-file => $seqdata,<br>
-format => 'fasta');<br>
my %seqs;<br>
my @prots;<br>
<br>
my $output;<br>
# process each sequence<br>
while( my $seq = $seqIO->next_seq ) { <br>
$seqs{$seq->display_id} = $seq;<br>
# translate them into protein<br>
my $protein = $seq->translate();<br>
my $pseq = $protein->seq();<br>
if( $pseq =~ /\*/ && <br>
$pseq !~ /\*$/ ) { <br>
warn("provided a cDNA sequence with a stop codon, PAML will
choke!");<br>
exit(0);<br>
}<br>
# Tcoffee can't handle '*' even if it is trailing<br>
$pseq =~ s/\*//g;<br>
$protein->seq($pseq); <br>
push @prots, $protein;<br>
}<br>
<br>
if( @prots < 2 ) {<br>
warn("Need at least 2 cDNA sequences to proceed");<br>
exit(0);<br>
}<br>
<br>
open(OUT, ">align_output.txt") ||<br>
die("cannot open output $output for writing"); <br>
# Align the sequences with clustalw<br>
<br>
my $aa_aln = $aln_factory->align(\@prots);<br>
<br>
# project the protein alignment back to cDNA coordinates<br>
my $dna_aln = &aa_to_dna_aln($aa_aln, \%seqs);<br>
<br>
my @each = $dna_aln->each_seq(); <br>
<br>
my $kaks_factory = new Bio::Tools::Run::Phylo::PAML::Codeml<br>
( -params => { 'runmode' => -2,<br>
'seqtype' => 1,<br>
'model' => 1, <br>
}<br>
);<br>
<br>
# set the alignment object<br>
$kaks_factory->alignment($dna_aln);<br>
<br>
# run the KaKs analysis<br>
my ($rc,$parser) = $kaks_factory->run();<br>
my $result = $parser->next_result; <br>
my $MLmatrix = $result->get_MLmatrix();<br>
<br>
my @otus = $result->get_seqs();<br>
# this gives us a mapping from the PAML order of sequences back to <br>
# the input order (since names get truncated)<br>
my @pos = map { <br>
my $c= 1;<br>
foreach my $s ( @each ) {<br>
last if( $s->display_id eq $_->display_id );<br>
$c++;<br>
}<br>
$c; <br>
} @otus; <br>
<br>
print OUT join("\t", qw(SEQ1 SEQ2 Ka Ks Ka/Ks PROT_PERCENTID
CDNA_PERCENTID)), "\n"; <br>
for( my $i = 0; $i < (scalar @otus -1) ; $i++) {<br>
for( my $j = $i+1; $j < (scalar @otus); $j++ ) {<br>
my $sub_aa_aln = $aa_aln->select_noncont($pos[$i],$pos[$j]);<br>
my $sub_dna_aln = $dna_aln->select_noncont($pos[$i],$pos[$j]); <br>
print OUT join("\t", <br>
$otus[$i]->display_id,<br>
$otus[$j]->display_id,$MLmatrix->[$i]->[$j]->{'dN'},<br>
$MLmatrix->[$i]->[$j]->{'dS'}, <br>
$MLmatrix->[$i]->[$j]->{'omega'},<br>
sprintf("%.2f",$sub_aa_aln->percentage_identity),<br>
sprintf("%.2f",$sub_dna_aln->percentage_identity), <br>
), "\n";<br>
}<br>
}<br>
<br>
<br>
<div><span class="gmail_quote">On 5/29/07, <b
class="gmail_sendername">Himanshu Ardawatia</b> <<a
href="mailto:himanshu.ardawatia@bccs.uib.no">himanshu.ardawatia@bccs.uib.no
</a>> wrote:</span>
<blockquote class="gmail_quote"
style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Hi
Xianjun,<br>
<br>
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. <br>
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. <br>
<br>
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. <br>
<br>
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.... <br>
<br>
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.<br>
<br>
Cheers and Fu<br>
Himanshu <br>
\\
<div><span class="e" id="q_112d6f0a9baae3a7_1"><br>
<div><span class="gmail_quote">On 5/28/07, <b
class="gmail_sendername">Dong Xianjun</b> <<a
href="mailto:Xianjun.Dong@bccs.uib.no" target="_blank"
onclick="return top.js.OpenExtLink(window,event,this)">
Xianjun.Dong@bccs.uib.no</a>> wrote:</span>
<blockquote class="gmail_quote"
style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">HI,
Himanshu<br>
<br>
I am sure you did some work in Ka/Ks calculation. Here I have a question<br>
bothering me; the output for Bio::Tools::Run::Phylo::PAML::Codeml is not<br>
stable(different for each runtime), and also different from the output <br>
with modeul of Bio::Tools::Run::Phylo::PAML::Yn00.<br>
<br>
Here I attached the script. Could you help to have a look and try to run<br>
the script? How is your way to calculate the Kaks ratio?<br>
<br>
Thanks<br>
<br>
--<br>
---------------------------<br>
Sterding (Xianjun) Dong<br>
PhD student, Boris Lenhard's group<br>
Bergen Center of Computational Science<br>
Bergen University, Norway<br>
Mobile: 0047-47361688<br>
Telephone: 0047-55276381 <br>
Skype: xianjun.dong<br>
<br>
<br>
</blockquote>
</div>
<br>
</span></div>
</blockquote>
</div>
<br>
</blockquote>
<br>
<pre class="moz-signature" cols="72">--
---------------------------
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
</pre>
</body>
</html>