Differences

This shows you the differences between two versions of the page.

--- tutorials:eccb_t2_codeml [2012/09/07 13:59]
romainstuder
+++ tutorials:eccb_t2_codeml [2012/09/09 10:47] (current)
romainstuder
@@ Line 1: / Line 1: @@
 ==== CodeML/PAML ====
-Lysozyme
-In this post, I will make a short tutorial on one of my favourite programs, CodeML, which is definitely not the easiest to use.
+=== Install of PAML ===
-=== Theoretical principles ===
+Download it from here: http://abacus.gene.ucl.ac.uk/software/paml.html#download
+Windows (The bin/ folder already contains windows executables):
+<code>
+cd paml4.6/
+</code>
+Linux:
+<code>
+cd paml4.5/
+rm bin/*.exe
+cd src
+make -f Makefile
+ls -l
+rm *.o
+mv baseml basemlg codeml pamp evolver yn00 chi2 ../bin
+chmod +x ../bin/*
+</code>
+MacOSX:
+<code>
+cd paml4.6/
+rm bin/*.exe
+cd src
+cc -O2 -o baseml baseml.c tools.c -lm
+cc -O2 -o basemlg basemlg.c tools.c -lm
+cc -O2 -o codeml codeml.c tools.c -lm
+cc -O2 -o pamp pamp.c tools.c -lm
+cc -O2 -o mcmctree mcmctree.c tools.c -lm
+cc -O2 -o evolver evolver.c tools.c -lm
+cc -O2 -o yn00 yn00.c tools.c -lm
+cc -O2 -o chi2 chi2.c -lm
+ls -l
+rm *.o
+mv baseml basemlg codeml pamp evolver yn00 chi2 ../bin
+chmod +x ../bin/*
+</code>
+=== Theoretical principles of the Branch-site model ===
 The selective pressure in protein coding genes can be detected within the framework of comparative genomics. The selective pressure is assumed to be defined by the ratio (ω) dN/dS. dS represents the synonymous rate (changing  the amino acid) and  dN the non-synonymous rate (keeping the amino acid). In the absence of evolutionary pressure, the synonymous rate and the non-synonymous rate are equal, so the dN/dS ratio is equal to 1. Under purifying selection, natural selection prevents the replacement of amino acids, so the dN will be lower than the dS, and dN/dS < 1. And under positive selection, the replacement rate of amino acid is favoured by selection, and dN/dS > 1.
@@ Line 15: / Line 54: @@
 It estimates various parameters (Ts/Tv, dN/dS, branch length) on the codon (nucleotide) alignment, based on a predefined topology (phylogenetic tree).
-Different codon models exist in CodeML. The model 0 estimates a unique dN/dS ratio for the whole alignment. Not really interesting, except to define a null hypothesis to test against. The branch models estimate different dN/dS among lineages (ie ASPM, a gene expressed in the brain of primates). The site models estimate different dN/dS among sites (ie in the antigen-binding groove of the MHC). The branch-site models estimate different dN/dS among sites and among branches. It can detect episodic evolution in protein sequences, as in the interactions between chains in the avian MHC. In my opinion, this is the most powerful application and this is the one used in the Sectome database (to which I contributed during my PhD).
+Different categories of codon models exist in CodeML:
+  * The model 0 estimates a unique dN/dS ratio for the whole alignment. Not really interesting, except to define a null hypothesis to test against. The other branch models estimate different dN/dS among lineages (ie ASPM, a gene expressed in the brain of primates).
+  * The site models estimate different dN/dS among sites (ie in the antigen-binding groove of the MHC).
+  * The branch-site models estimate different dN/dS among sites and among branches. It can detect episodic evolution in protein sequences. IMHO, the most interesting model.
 First, we have to define the branch where we think that position could have occurred. We will call this branch the "foreground branch" and all other branches in the tree will be the "background" branches. The background branches share the same distribution of ω = dN/dS value among sites, whereas different values can apply to the foreground branch.
@@ Line 55: / Line 97: @@
 Let's go in details.
-File Preparation:
+=== Identification of positive seleciton in the serine/threonine-protein kinase gene family ==
-We need four files to run CodeML:
+In the evolution vertebrates, we would like to know if the branch leading to the Teleost fishes (genes A50 to A54) in the serine/threonine-protein kinase PAK 2 gene was under positive selection or not.
+And if yes, which residues were under positive selection.
-  - The multiple nucleotide (CDS) alignment, in PHYLIP format. CodeML will strictly remove any position that contains at least one gap or an unknown "N" nucleotide: TF105351.Eut.3.phy
+We need four files to run CodeML (unzip them all):
-  - The phylogenetic tree in newick format, with the branch of interest specified by "#1"(You can view it with NJplot or FigTree): TF105351.Eut.3.53876.tree
-  - A command file where all parameters to run CodeML under the alternative model are specified: TF105351.Eut.3.53876.ctl
+  - The multiple nucleotide (CDS) alignment, in PHYLIP format. CodeML will strictly remove any position that contains at least one gap or an unknown "N" nucleotide: {{:tutorials:tf105351.eut.3.phy.zip|}}
-  - A command file where all parameters to run CodeML under the null model are specified: TF105351.Eut.3.53876.fixed.ctl
+  - The phylogenetic tree in newick format, with the branch of interest specified by "#1"(You can view it with NJplot or FigTree): {{:tutorials:tf105351.eut.3.53876.tree.zip|}}
+  - A command file where all parameters to run CodeML under the alternative model are specified: {{:tutorials:tf105351.eut.3.53876.ctl.zip|}}
+  - A command file where all parameters to run CodeML under the null model are specified: {{:tutorials:tf105351.eut.3.53876.fixed.ctl.zip|}}
+The tree looks like:
+{{:tutorials:tree.png|}}
 == Execute CodeML ==
@@ Line 68: / Line 116: @@
 Run command file (alternative model):
-We estimate the Ts/Tv ratio (fix_kappa = 0) and the dN/dS (fix_omega = 0). The branch-site model is specified by setting the model parameter to 2 (different dN/dS for branches) and the NSsites value to 2 (which allows 3 categories for sites: purifying, neutral and positive selection).
+We estimate the Ts/Tv ratio (fix_kappa = 0) and the dN/dS (fix_omega = 0). The branch-site model is specified by setting these two parameters:
+    * model = 2 (different dN/dS for branches)
+    * NSsites value to 2 (which allows 3 categories for sites: purifying, neutral and positive selection).
 <code>
@@ Line 146: / Line 196: @@
 codeml ./TF105351.Eut.3.53876.ctl
 codeml ./TF105351.Eut.3.53876.fixed.ctl
-<./code>
+</code>
 In Windows:
@@ Line 154: / Line 204: @@
 </code>
+Two mlc output files are produced (as it can take time, you can download them directly in the next step).:
+  * {{:tutorials:tf105351.eut.3.53876.mlc.zip|TF105351.Eut.3.53876.mlc (alternative model)}}
+  * {{:tutorials:tf105351.eut.3.53876.fixed.mlc.zip|TF105351.Eut.3.53876.fixed.mlc (null model)}}
 === Analyse results ===
-) Assign significance of the detection of positive selection on the selected branch:
+== Step 1) Assign significance of the detection of positive selection on the selected branch: ==
-Two output files are produced:
-TF105351.Eut.3.53876.mlc (alternative model) and TF105351.Eut.3.53876.fixed.mlc (null model).
+We retieve the likelihood values lnL1 and lnL0 from TF105351.Eut.3.53876.mlc and TF105351.Eut.3.53876.fixed.mlc files, respectively.
-We retieve the likelihood values lnL1 and lnL0 from TF105351.Eut.3.53876.mlc and TF105351.Eut.3.53876.fixed.mlc files, respectively:
+We retieve the number of parameters np1 and np0 from TF105351.Eut.3.53876.mlc and TF105351.Eut.3.53876.fixed.mlc files, respectively.
-   * lnL(ntime: 41  np: 46):  -4707.210163      +0.000000  (lnL1)
+   * lnL(ntime: 41  np: 46):  **-4707.210163**      +0.000000  (lnL1)
-   * lnL(ntime: 41  np: 45):  -4710.222252      +0.000000  (lnL0)
+   * lnL(ntime: 41  np: 45):  **-4710.222252**      +0.000000  (lnL0)
-We can construct the LRT (you can use your favourise spreadsheet for that. Or even better: R):
+We can construct the LRT (you can use your favourite spreadsheet for that. Or even better with R):
 ΔLRT = 2×(lnL1 - lnL0) = 2×(-4707.210163 - (-4710.222252)) = 6.024178
@@ Line 175: / Line 227: @@
 The degree of freedom is 1 (np1 - np0 = 46 - 45).
-p-value = 0.014104 (under χ²) => significant.
+(With OpenOffice: 1-CHISQDIST(6.024178;1))
+p-value = 0.014104 (under χ²) => __significant!__
 A significant result with the branch-site codon model means that positive selection affected a subset of sites during a specific evolutionary time (also called episodic model of protein evolution).
-) If significant, we can retrieve sites under positive selection:
+== Step 2) If significant, we can retrieve sites under positive selection ==
 In the TF105351.Eut.3.53876.mlc, we can retrieve sites under positive selection using the Bayes Empirical Bayes (BEB) method:
@@ Line 192: / Line 246: @@
    * Position 159 has a very high probability (99.3%) of being under positive selection.
+You can visualise the multiple in Jalview.
+- Open Jalview
+- Load TF105351.Eut.3.phy
+- Then: Calculate->translate cDNA. (tips: by moving the pointer on the amino acids alignment, you can see the corresponding codon in the nucleotide alignment.
+{{:tutorials:tf105351.eut.3.aln.png|}}
+=== Using other models ===
+Other models can be tested by changing these parameters model and NSsites.
+== Example 1: ==
+Site model M1 (neutral):
+    * model = 0 (dN/dS doesn't vary on branches)
+    * NSsites = 1 (which allows 2 categories for sites: purifying and neutral).
+Site model M2a (positive selection):
+    * model = 0 (dN/dS doesn't vary on branches)
+    * NSsites = 2 (which allows 3 categories for sites: purifying, neutral and positive selection).
+Then we can compare M1 and M2a by the likelihood ratio test.
+== Example 2: ==
+Branch model M0:
+    * model = 0 (dN/dS doesn't vary on branches)
+    * NSsites = 0 (dN/dS doesn't vary on sites).
+Branch model M2 with different dN/dS (positive selection on selected branches):
+    * model = 2 (different dN/dS for branches)
+    * NSsites = 2 (dN/dS doesn't vary on sites).
+Then we can compare M0 and M2 by the likelihood ratio test.

Trace:

CATH Documentation

Differences

Search

Navigation

Print/export

Toolbox