Codeml: 计算dN/dS

0. Note

对于两个物种来源的两条同源基因的CDS序列，我们分别标记为CDS-1和CDS-2。对于这两条CDS，仅有互相之间的dN/dS，也就是 dN/dS between these two CDS。在计算之前，要把两条CDS序列进行对齐，然后再计算，因此CDS-1作为参照求得的CDS-2的dN/dS与CDS-2作为参照求得的CDS-1的dN/dS完全相同。根本原因是，非同义突变在无论哪条CDS做参照的情况下都是非同义突变，而同义突变在无论哪条CDS做参照的情况下也都是同义突变。

0. Paper

https://academic.oup.com/mbe/article/40/4/msad041/7140562?login=true

0. HitHub

https://github.com/abacus-gene/paml-tutorial

1. 官网

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

2. 官方文档

Manual:

http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf
https://blog.csdn.net/winglyx/article/details/6731787

Q&A:

http://abacus.gene.ucl.ac.uk/software/pamlFAQs.pdf

3. 教程

https://zhuanlan.zhihu.com/p/34463038
https://zhuanlan.zhihu.com/p/105159767
https://www.jianshu.com/p/372cf1f02d80
https://www.plob.org/article/21094.html
https://www.plob.org/article/6943.html
https://www.jianshu.com/p/1c4e40d075c6

4. 比较两个物种的一组同源基因

主要参考：https://zhuanlan.zhihu.com/p/34463038

需要的文件：

sequence.cds: 两个物种各自一个 CDS 序列

sequence.protein: 两个物种各自一个 protein 序列

tree.txt: 两个物种的 进化树 文件；如果是两个物种，则只需 (Species_name_1, Species_name_2) 即可

control.cnt: codeml 需要的控制文件

需要的软件：

clustalw2 (clustal 包): 用来比对氨基酸序列

pal2nal.pl (pal2nal 包): 用来比对核酸序列，并将序列比对结果导出为 PAML 需要的格式

codeml (paml 包): 用来进行 dN/dS 的计算

流程：

• 1. 对比此组同源基因的 protein 序列：clustalw2 sequence.protein
     此过程产生两个文件：sequence.aln + sequence.dnd
     其中 sequence.aln 需要在下一步进行使用
• 2. 对比此组同源基因的 cds 序列并输出为 PAML 要求的格式：
     pal2nal.pl sequence.aln sequence.cds -output paml -nogap > Two.codon
• 3. 进行比对：codeml control.cnt
     要求：Two.codon, tree.txt 和 control.cnt 三个文件需要放在同一目录的同一层

control.cnt : codeml 使用的控制文件

seqfile = Two.codon : 作为输入的比对后的序列文件
treefile = tree.txt : 进化树文件
outfile = output : 输出内容存放的文件；输出内容会放入和此控制文件同一层目录

noisy = 0 * 0,1,2,3,9: how much rubbish on the screen
verbose = 0 * 1: detailed output, 0: concise output
runmode = -2 * 0: user tree; 1: semi-automatic; 2: automatic
* 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise : 重点参数

cleandata = 1 * "I added on 07/07/2004" Mikita Suyama

seqtype = 1 * 1:codons; 2:AAs; 3:codons-->AAs
CodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
model = 0 * models for codons: * 0:one, 1:b, 2:2 or more dN/dS ratios for branches : 重点参数

NSsites = 0 * dN/dS among sites. 0:no variation, 1:neutral, 2:positive : 重点参数
icode = 0 * 0:standard genetic code; 1:mammalian mt; 2-10:see below
Mgene = 0 * 0:rates, 1:separate; 2:pi, 3:kappa, 4:all

fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated
kappa = 2 * initial or fixed kappa
fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate
omega = 1 * initial or fixed omega, for codons or codon-transltd AAs

fix_alpha = 1 * 0: estimate gamma shape parameter; 1: fix it at alpha
alpha = .0 * initial or fixed alpha, 0:infinity (constant rate)
Malpha = 0 * different alphas for genes
ncatG = 4 * # of categories in the dG or AdG models of rates

clock = 0 * 0: no clock, unrooted tree, 1: clock, rooted tree
getSE = 0 * 0: don't want them, 1: want S.E.s of estimates
RateAncestor = 0 * (1/0): rates (alpha>0) or ancestral states (alpha=0)
method = 0 * 0: simultaneous; 1: one branch at a time

重点参数讲解：

runmode=: 决定比对(所用树)的模式

0: 使用用户输入的 tree 文件
-2: 成对比对，当仅有两个物种时使用此模式；本示例适合此参数

model=: 决定 dN/dS 计算的模式

0: 计算所有物种互相比对后的 均值
2: 对每个分支都计算一个值；此时需要指出 目标分支 和 背景分支

NSsites=: 一般选 0

参数选择：

对应 M0 模式，即所有位点和所有谱系的 ω 比为 1

runmode= -2
model= 0
NSsites= 0

结果文件解读：

其中的 ω 表示所有分支共同 (平均) 的 dN/dS
https://blog.sciencenet.cn/home.php?mod=space&uid=3434047&do=blog&id=1389140

5. 计算两个物种的多组同源基因

https://zhuanlan.zhihu.com/p/105159767

6. 计算各个分支独特

https://blog.sciencenet.cn/blog-3433349-1241310.html

omega 0/ ω0 表示背景 dN/dS

7. Model 选择

https://www.jianshu.com/p/1c4e40d075c6

image.png

8. 计算多个物种的多组同源基因

https://zhuanlan.zhihu.com/p/105159767

tree.txt 文件的构建和相应 codeml.ctl 文件的书写

https://github.com/abacus-gene/paml-tutorial/tree/main/positive-selection/02_extra_analyses

https://academic.oup.com/mbe/article/40/4/msad041/7140562?login=true