codeml的理论和实践

最近被数据分析虐的死去活来，好难。

参考文献：

Likelihood ratio tests for detecting positiveselection and application to primate Lysozyme Evolution

User guide of Paml version 4.9

理论部分

dS ：synonymous substitution 同义替换

dN ：non-synonymous substitution 非同义替换

同义替换指：由于密码子具有简并性，密码子中一些核苷酸的替换并不会导致编码氨基酸的变化。非同义替换指：密码子中核苷酸的替换导致了编码氨基酸的变化。

dN/dS：非同义替换和同义替换速率的比值

Darwinian selection may have a non-synonymous/synonymousrate ratio (dN/dS) that is different fromthose of other lineages or greater than one. (positive selection)

基因的变异只有表达出来（编码蛋白发生改变），自然选择才能起作用。

如果一个分支的基因变异受到自然选择的正选择（positive selection）（正选择指的是自然选择支持变异，变异是有利的），那么该分支的dN/dS应该与其他分支不同，或者大于1。

需要注意的是某基因在某一分支的dN/dS大于1来证明该基因在该分支上受正选择是很保守的，dN/dS很少会大于1，如果某一个分支的dN/dS值显著大于其他分支，也能说明问题。

如果分支的某基因选择上中性（neutral selection），即无害也无利，那么理论上dN/dS=1

但大多数情况下，dN/dS<1，即非同义 替换率小于同义替换率

The basic model for the likelihood analysisis a version of the codon-substitution model of Goldman and Yang (1994) andaccounts for the genetic code structure, transition/transversion rate bias, anddifferent base frequencies at codon positions.

codeml里用于最大似然分析的基础模型（basal model）基于一个密码子替换模型，该模型考虑了遗传序列的结构、转换/颠换率的偏差、以及密码子不同位置替换的基础频率。

从密码子i到密码子j的瞬时替换率为：

K代表转换率/颠换率（transition/transversion）

ω 代表dN/dS

The ω ratio is a measure of natural selection acting on the protein. Simplisitically, values of ω<1, =1, and >1 means negative purifying sleection, neutral selection, and positive selection.

πj 代表密码子j的平衡频率，由密码子三个位点的核苷酸频率估算得到

由这一个模型衍生出不同的模型，这些模型具有不同水平的dN/dS异质性。

However, the ratio averaged over all sites and all lineages is alomost never >1, since positive selection is unlikely to affect over all sites and all lineages over prolonged time. Thus, interest has been focused on detecting positive selection that affects only some lineages or some sites.

分支模型branch model

The branch models allow the ω ratio to vary among branches in the phylogeny and are useful for detecting positive selection acting on particular lineages.

分支模型允许ω在系统发育树上的不同分支上存在变异，适用于检验编码基因在特定分支是否存在正选择。

最简单的模型是one-ratio模型：假设所有的分支都有共同的dN/dS

另一个极端是Free-ratio模型：每个分支都有独立的dN/dS

还有很多中间类型的模型，比如two-ratio模型：假设感兴趣的分支具有的非同义替换和同义替换比率（ω1）和背景的非同义替换和同义替换比率（ω0）不同。

The above models can be compared using thelikelihood ratio test to examine interesting hypothesis.

上述模型可通过似然比检验进行比较，验证感兴趣的假说。

比如one-ratio模型和two-ratio模型就可以一起比较，验证某一分支的dN/dS是否和其他分支不同，这里two-ratio模型就是备择假设(alernative hypothesis),one-ratio模型就是零假设（null hypothesis）

位点模型site models

下面不会讲到位点模型的实际操作，这里简单介绍一下。

The site models allow the ω ratio to vary among sites (among codons or amino acids in the protein).

位点模型允许不同位点（密码子之间或蛋白的氨基酸之间）上的ω存在差异。

分支-位点模型brach-site models

The branch-site models allow ω to vary among sites in the protein and across branches on the tree and aim to detect positive selection affecting a few sites along particular lineages.

分支-位点模型允许编码基因的不同位点在系统发育树的不同分支的ω存在变异,目的是检验正选择是否作用于特定分支的一些位点。

Easycodeml里的分支位点模型用到的是下面这篇文章里提到的Test2：

Evaluation of an improved branch-site likelihoos method for detecting positive selection at the molecular level.

Test2, also known the branch-site test of positive selection, compares the modified model A with the corresponding null model with ω2=1 fixed (fix_omege=1 and omega=1)

Test2里用于似然比检验的模型就是model A和对应的null model。

先来介绍下model A：

modelA在似然比检验中是备择假设。

在modelA里，系统发育树被预先分为前景枝和背景枝，只有前景枝才经历正选择。模型假设有四个类型的位点（site class）：

class 0：包括在整个系统发育树中保守的密码子，dN/dS估算值:0<ω0<1

class 1:包括在整个系统发育树中演化上中性的密码子，dN/dS估算值:ω1=1

class 2a and 2b:包括在背景枝保守(2a，0<ω0<1)，在前景枝受正选择（dN/dS估算值:ω2>1）的密码子

class 2b:包括在背景枝演化上中性（2b，ω1=1），在前景枝受正选择（dN/dS估算值:ω2>1）的密码子

Test2中，和modelA进行似然比检验的零模型（null model）就是将modelA中ω2值固定，ω2=1。

在modelA里面，一些位点（2a和2b）的前景枝的ω值预设为ω2>1,因此Test2也是对前景枝某些位点是否经历正选择的直接检验。

进行似然比检验的方法是卡方分布。

实操部分

实操部分介绍两个方法做分支模型的检验，

一个是windows系统下用可视化软件EasyCodeML，

一个是在linux系统下用paml的codeml。

1 EasyCodeML

更详细的介绍可以参考高芳銮老师的这篇文章：

https://blog.sciencenet.cn/blog-460481-1163040.html

EasyCodeML：A visualtool for analysis of selection using CodeML

Paml: a program package for phylogeneticanalysis by maximum likelihood

Paml 是使用最大似然法进行系统发育分析的程序包，包括的程序有baseml、codeml、evolver、basemlg等等。

Paml：http://abacus.gene.ucl.ac.uk/software/paml.html

而EasyCodeML是使用codeml进行选择分析的可视化工具，也就是说EasyCodeML可以实现Paml里面codeml程序的功能，而且提供可视化操作界面。

EasyCodeML：https://github.com/BioEasy/EasyCodeML/

下面是利用EasyCodeML中的分支模型（branch model）检验感兴趣的分支是否具有和其他分支不同的非同义替换和同义替换比率(dN/dS)的实践。

Easycodeml

1 ) 安装软件

软件的下载地址：https://github.com/BioEasy/EasyCodeML/

下载的压缩包解压之后，就可以直接双击打开其中的应用程序使用。

2 ) 准备文件

需要准备的文件包括序列文件和树文件，如果是在paml里运算，还需要一个控制文件，但在EasyCodeML里不需要

序列文件就是不同分类群的要检测的基因序列，进行translation align（这种比对方法是先将碱基翻译成密码子，然后将密码子作为一个整体比对）。比对后，删去终止密码子。如果有gap，我采用的方法是直接删去比对矩阵中有gap的整列。

上方工具栏Tools-Seqformat Convertor:将比对序列的.fasta或其他格式的文件转化为.paml文件

直接将文件拖入original file框中，右边三角下来菜单中选择转换前的文件格式，生成的.paml文件会出现在转换前的文件所在的文件夹中。

树文件就是要研究的分类群的系统发育树文件，需要newick格式。分析时用到的树的拓扑结构，而不用枝长信息。要求树的拓扑结构能反映真实的系统发育关系，因此这个系统发育树不一定是基于要研究的基因建立的树。

The tree topology, but not the branch length is used to fit different models.

需要注意的是序列文件和树文件中的分类单元名称要一致，不要出现空格或分号。

3 ) 运行EasyCodeML

runing model

preset是预设模式，适合不太熟悉分析的新手。

custom适合熟悉分析，希望自己设置参数的用户。

Setup

Model Selection 根据目的选择，本例这里选branch model

working directionary 指定生成文件所在路径

If the option “Clean data” was selected, all sites with ambiguity characters and alignment  gaps will be removed from the aligned sequence file. Note that alignment gaps are treated as missing data.

clean data 如果选择clean data选项，对比序列中存在的gaps和ambiguity characters会被移除掉。

Aligned sequence file in palm format 指定paml格式的序列文件所在的路径，

icode下拉选项中可以选择编码基因的序列类型，比如这里用的是昆虫的线粒体基因，选择4 invertebrate mt.

Tree File in Newick Format:指定newick格式的树文件所在的路径，点击Check，检查树文件和序列文件的分类群名称是否一致，点击label对树文件中感兴趣的分支进行标记：

选中感兴趣的分支，标记成功后该分支会变黄，并带有#1，也就是前景枝。这里只有一个前景枝，如果有两个，选中第二个分支，点击2nd，以此类推。其他的分支称为背景枝

在自己的树文件所在的文件夹下，会自动生成一个在原树文件名后带后缀labeled的树文件（如果有需要，可用于paml程序包的分析）

最后点击Run CodeML

运算过程可以在Runing status中查看，如果遇到错误也可以根据报错结果解决。

4 ）结果解读

程序运行结束后，会弹出提示任务完成的小窗口。

Summary of results处点击export，给出结果(excel表格)的输出路径。结果呈现如下图：

在整个运算过程中，用到了两个branch model，

一个是model0，也就是之前提到过的one-ratio模型：假设所有的分支都有共同的dN/dS

一个是two ratio model 2:假设感兴趣的分支和其他分支有不同的dN/dS

在two ratio model 2中，ω1代表的是前景枝的dN/dS的估算值，ω0代表的是背景支的dN/dS的估算值

在Model 0中，ω代表的是在假设整个树只有一个dN/dS的前提下，dN/dS的估算值

likelihood ratio test结果：

lnL代表的是模型的最大似然值的log值：

以上图的结果为例，Two-ratio model的lnL值为-5820.843753，one-ratio model的lnL值为 -5822.141900，二者差值的二倍（1.20x2）与自由度为1（df=1）的卡方分布进行比较。

自由度的值为两个模型参数个数（np, number of parameters）的差值.

p-value<0.05（显著significant）或 p-value<0.01（极显著extremely significant），说明相比于one-ratio model，Two-ratio model更符合数据，选择备择假设（Two-ratio model），

原始结果文件.mlc在指定的工作路径下查看。

2 paml

以下操作都是在linux系统下进行：

1）安装paml

conda安装：

conda install -c bioconda palm

2) 准备文件

序列文件的准备可以参考paml的案例文件，这是我用的序列文件（我偷懒了，用Easycodeml里面的Seqformat Convertor转化.fasta文件为.paml文件，再将后缀改为.nuc，顺便说下linux系统批量改文件后缀的命令：rename .paml .nuc  *.paml）。

最上方的数字：59代表分类单元数，120代表比对序列长度（核苷酸长度）

树文件的话需要准备两个，都为.newick文件，一个是没有标注感兴趣分支的，一个是标注感兴趣分支的（我又偷懒了，标注的树文件也是通过easycodeml标注的）

控制文件的内容如下，本例中我们需要两份控制文件，分别是one ratio model的控制文件和two ratio model的控制文件。

最上方的seqfile、treefile、outfile分别是告诉程序你的序列文件、树文件和最后结果文件的绝对路径

需要注意的几个参数是：

runmode = 0，意思是所用的树文件是用户提供的树。（因为基于所分析基因构建的系统发育树往往不能反映真实的系统发育关系）

icode = 4 代表分析的invertebrate mt（无脊椎动物的线粒体基因）

fix_kappa = 0， 代表kappa根据数据进行估算（kappa to be estimated） ,kappa代表转换率/颠换率（transition/transversion）

mode = 0， 代表的是该控制文件是运算one ratio model的控制文件，即假设所有的分支都有共同的dN/dS

mode = 2，代表该控制文件是运算具有多个ratio model的控制文件，即假设分支上具有两个及以上的dN/dS

3）运行paml中的Codmel程序

我是在控制文件所在的路径下运行该程序。

#查看codeml的路径：

which codeml

#运行程序：one ratio model

~/miniconda3/envs/paml/bin/codeml ATP8_M0_codeml.ctl

在上面的命令中，~/miniconda3/envs/paml/bin/codeml是codeml程序路径，ATP8_M0_codeml.ctl是控制文件

#运行程序：two ratio model

~/miniconda3/envs/paml/bin/codeml ATP8_BM_codeml.ctl

运算结束后，结果会出现在控制文件指定的输出路径下

4）结果解读

在结果文件中找到以下数据(这里以one ratio model 运算结果为例)：

lnL value of one ratio model
omage of one ratio model

模型的lnL值，omega值（ω）。自己可以进行最大似然比检验，或者用Easycodeml上方工具栏中Tools下面的最大似然比检验工具LRT's calculator

其解读与上述Easycodeml部分相同。