PAML是一个程序包，用于通过最大似然法分析DNA或蛋白质序列的系统发育关系。它由Ziheng Yang维护，并根据GNU GPL v3分发。它为UNIX/Linux/Mac OSX提供了ANSI C源代码，并为MS Windows提供了可执行文件。PAML并不擅长构建树。当您使用其他程序（如PAUP*、PHYLIP、MOLPHY、PhyML、RaxML等）构建了树之后，它可用于估计参数和检验假设以研究进化过程。

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

一，PAML可以做什么？
1.系统发育树的检验和比较
2.复杂替代模型中参数的估计
3.对不同假设模型进行似然率假设性检验
4.基于全局或局部分子钟模型估算分歧时间
5.使用核酸，氨基酸和密码子模型重构祖先序列
6.Monte Carlo模拟生成核酸，密码子和氨基酸序列数据集
7.估算同义和非同义替代率和检测蛋白编码DNA序列中的正选择
8.贝叶斯估算物种的分歧时间来合并化石尺度上的不确定性

二.为什么要计算dN/dS?
蛋白编码基因的自然选择压力水平可以通过dN/dS值的大小来衡量，其中，dS代表同义替换率（synonymous rate），dN代表非同义替换率（non-synonymous rate）。在没有受到选择压力时，同义替换率和非同义替换率相等，此时dN/dS = 1；当受到负选择或净化选择压力时，自然选择会阻止氨基酸发生改变，同义替换率会大于非同义替换率，即dN/dS < 1;当受到正选择压力时，氨基酸的置换率会受自然选择的青睐，即蛋白功能可能会发生适应性改变，此时dN/dS > 1。

三，什么是同义替换和非同义替换？
同义替换 (synonymous):大多数是中性的，既没有正向选择也没有负向选择（因为虽然核苷酸变了，但是氨基酸没有发生改变）
非同义替换（non-synonymous）：大多数是吃亏的，会对个体产生一定不利的影响（氨基酸一变，蛋白质可能就变了，结构功能等就都不一样了）

四，其他背景知识
ω是dN/dS的估计值，通过最大似然法计算得到
零假设：所有枝系有共同的ω
备择假设：一些枝系有不同的ω

1.fasta文件：fasta格式的多序列比对文件，要经过mafft，prank等软件进行比对，在经过gblocks,trimal等质量过滤软件处理后获得，序列长度必须为3的倍数且不允许存在终止密码子。（前期已经使用ParaAT脚本处理完成）

2.tree文件：nwk格式的树文件，主要包括分析涉及物种的系统发育关系，物种名必须与序列比对里面的物种名保持一致，否则会报错；树文件要以分号；结尾。

示例文件：

(JCY4_JCY4_GM004662:0,((PA4_PA4_052656075.1:0.0029115,PA4_PA4_052656077.1:0.560755):0.604171,JCYB8_JCYB8_GM004310:0):0);

3.codeml.ctl: paml软件的参数控制文件，通过修改该文件参数，实现模型切换

codeml.ctl文件详细说明：

输入输出参数：
seqfile: 输入的多序列比对文件路径
treefile: 输入的树文件路径
outfile: 输出文件路径
noisy: 输出到屏幕上的信息量等级：0,1,2,3,9
verbose: 输出到结果文件中的信息量等级：0，精简模式结果（推荐）；1，输出详细信息，包括碱基序列；2，输出更多信息
getSE: 是否计算并获得各参数的标准误：0，不需要；1，需要
RateAncestor: 是否计算序列中每个位点的替换率：0，不需要；1，需要；当设置成1时，会在结果文件rates中给出各个位点的碱基替换率；同时也会进行祖先序列的重构，在结果文件rst中有体现

数据使用说明参数：
runmode: 程序获取进化树拓扑结构的方法：0，直接从输入文件中获取进化树拓扑结构（推荐）；1，程序以输入文件中的多分叉树作为起始树，并使用星状分解法搜索最佳树；2，程序直接以星状树作为起始树，并使用星状分解法搜索最佳树；3，使用逐步添加法搜索最佳树；4，使用简约法获取起始树，再使用邻近分支交换寻找最优树；5，以输入文件中的树作为起始树，再使用邻近分支交换寻找最优树；-2，对密码子序列进行两两比较并使用ML方法计算DnDs,或对蛋白序列进行两两比较计算ML距离，而不进行其他参数（枝长和omega等）的计算。
fix_blength: 程序处理输入树文件中枝长数据的方法：0，忽略输入树文件中的枝长信息；-1，使用一个随机起始进行计算；1，以输入的枝长信息作为起始值进行ML迭代分析，此时需要注意输入的枝长信息是碱基替换率，而不是碱基替换数；2，不需要使用ML迭代计算枝长，直接使用输入的枝长信息，需要注意，若枝长信息和序列信息不吻合可能导致程序崩溃；3，让ML计算出的枝长和输入的枝长呈正比。
seqtype: 输入的多序列比对数据的类型：1，密码子数据； 2，氨基酸数据； 3，输入数据虽然为密码子序列，但先转换为氨基酸序列后再进行分析。
CodonFreq: 密码子频率的计算方法：0，表示除三种终止密码子频率为0，其余密码子频率全部为1/61; 1,程序分别计算三个密码子位点的四种碱基频率，再算四种碱基频率在三个位点上的算术平均值，计算61种密码子频率时，使用该均值进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；2，程序分别计算三个密码子位点的四种碱基频率，计算61种密码子频率时，使用相应位点的碱基频率进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；3，直接使用观测到的各密码子的总的频数、所有密码值的总数，得到所有密码子的频率。一般选择第三种方法，设置CodonFreq的值为2,。第一种方法让所有密码子频率相等，不符合实际情况；第二种方法没有考虑到三种位点各碱基频率的差异，得到的密码子频率不准确，；第三种方法能比较准确计算出各密码子的频率；第四种方法由输入数据得到真实的各密码子频率，但有时候有些非终止密码子的频率为0，可能不利于后续计算。
cleandata: 是否移除不明确字符（N、？、W、R和Y等)或含以后gap的列后再进行数据分析;0，不需要，但在序列两两比较的时候，还是会去除后进行比较；1，需要，**有可能会导致运行失败**。
ndata: 输入的多序列比对的数据个数。
clock: 进化树中各分支上的变异速率是否一致，服从分子钟理论：0，变异速率不一致，不服从分子钟理论（当输入数据中物种相差较远时，各分支变异速率不一致）；1，所有分枝具有相同的变异速率，全局上服从分子钟理论；2，进化树局部符和分子钟理论，程序认为除了指定的分支具有不同的进化速率，其它分枝都具有相同的进化速率，这要求输入的进化树信息中使用#来指定分支；3，对多基因数据进行联合分析。
Mgene: 是否有多个基因的多序列比对信息输入，以及多各基因之间的参数是否一致：0，输入的多序列比对文件中仅包含一个基因时，或多个基因具有相同的Kappa和Pi参数；1，输入文件包含多个基因，这些基因之间是相互独立的（这些基因之间具有不同的Kappa和Pi值，且其进化树的枝长也不相关）；2，输入文件包含多个基因，这些基因具有相同的Kappa值，不同的Pi值；3，输入文件包含多个基因，这些基因具有相同的Pi值，不同的Kappa值；4，输入文件包含多个基因，这些基因具有不同的Kappa值和不同的Pi值。当值是2、3或4时，多个基因的进化树枝长虽然长度不一样，但是呈正比关系。这点和参数值等于1是不一样的。
icode: 设置遗传密码。其值1-10和NCBI的1-11遗传密码规则对应：0，表示通用的遗传密码。
Small_Diff: 设置一个很小的值，一般位于1e-8到1e-5之间。推荐检测设置不同的值在比较其结果，需要该参数值对结果没什么影响。

位点替换模型参数：
model: 若输入数据是密码子序列，该参数用于设置branch models，即进化树各分枝的omega值的分布：0，进化树上所有分枝的omega值一致；1，对每个分枝单独进行omega计算；2，设置多类omega值，根据树文件中对分枝的编号信息来确定类别，具有相同编号的分枝具有相同的omega值，没有编号的分枝具有相同的omega值，程序分别计算各编号和没有编号的omega值。
          若输入数据是蛋白序列，或数据是密码子序列且seqtype值是3时，该参数用于设置氨基酸替换模型：0，Poisson；1，氨基酸替换率和氨基酸的观测频率成正比；2，从aaRatefile参数指定的文件路径中读取氨基酸替换率信息，这些信息是根据经验获得，并记录到后缀为.dat的配置文件中。这些经验模型(Empirical Models)文件位于PAML软件安装目录中的dat目录下，例如，   Dayhoff(dayhoff.dat)、WAG(wag.dat)、LG(lg.dat)、mtMAN(mtman.dat)和mtREV24(mtREV24.dat)等；
aaRatefile = dat/wag.dat: 当对蛋白数据进行分析，且model = 2时，该参数生效，用于设置氨基酸替换模型。
aaDist: 设置氨基酸之间的距离。
NSsites: 输入数据时密码子序列时生效，用于设置site model，即序列各位点的omega值的分布：0，所有位点具有相同的omega值；1，各位点上的omega值小于1或等于1（服从中性进化neutral）；2，各位点上的omega值小于1、等于1或大于1（选择性进化selection）；3，discrete；4，freq；5:gamma；6，2gamma；7，beta；8，beta&w；9，beta&gamma；10，beta&gamma+1；11，beta&normal>1；12，0&2normal>1；13，3normal>0。
            可以一次输入多个模型进行计算并比较，其结果输出的rst文件中。
fix_alpha: 序列中不同的位点具有不同的碱基替换率，服从discrete-gamma分布，该模型通过alpha（shape参数）和ncatG参数控制。该参数设置是否给定一个alpha值：0，使用ML方法对alpha值进行计算；1，使用下一个参数设置一个固定的alpha值。
               对于密码子序列，当NSsites参数值不为0或model不为0时，推荐设置fix_alpha = 1且alpha = 0，即不设置alpha值，认为位点间的变异速率一致，否则程序报错。若设置了alpha值，则程序认为不同密码子位点的变异速率不均匀，且同时所有位点的    omega值一致，当然各分枝的omega值也会一致，这时要求NSsites和model参数值都设置为0（这一般不是我们需要的分析，它不能进行正选择分析了）。
alpha: 设置一个固定的alpha值或初始alpha值（推荐设置为0.5）。该值小于1，表示只有少数热点位置的替换率较高；该值越小，表示位点替换率在各位点上越不均匀；若设置fix_alpha = 1且alpha = 0则表示所有位点的替换率是恒定一致的。
Malpha: 当输入的多序列比对结果中有多基因时，设置这些基因间的alpha值是否相等：0，分别对每个基因单独计算alpha值；1，所有基因的alpha值保持一致。
ncatG: 序列中不同位点的变异速率服从GAMMA分布的，ncatG是其一个参数，一般设置为5，4，8或10，且序列条数越多，该值设置越大。
          对于密码子序列，当NSites设置为3时，ncatG设置为3；当NSites设置为4时，ncatG设置为5；当NSites值设置>=5时，ncatG值设置为10。
fix_kappa:设置是否给定一个Kappa值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的Kappa值。
kappa: 设置一个固定的Kappa值，或一个初始的Kappa值。
fix_omega: 设置是否给定一个omega值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的omega值。 
omega: 设置一个固定的omega值，或一个初始的omega值。
method: 设置评估枝长的ML迭代算法：0，使用PAML的老算法同时计算所有枝长，在clock = 0下有效；1，PAML新加入的算法，一次对一个枝长进行计算，该算法仅在clock参数值为1，2或3下工作。

/TJPROJ7/META_ASS/16s/yaoyuanyuan/personalscript/micro/software/PAML4.9/paml4.8/bin/codeml /TJPROJ7/META_ASS/16s/yaoyuanyuan/personalscript/micro/software/PAML4.9/paml4.8/test/codeml.ctl

部分分析结果展示: