PAML软件中的codeml命令可以使用Maxmum Likelihood方法对密码子或氨基酸的多序列比对结果进行分析和进化参数计算 。这两种数据共用一个命令，程序运行时的参数绝大部分是公用的，但仍需要分清两种数据各自的参数意义。PAML软件中进行正选择分析的示例数据位于examples/lysozyme/目录下。
1.输入文件（1/3）：含有拓扑结构的树文件
输入 一个树文件input.trees。该文件可以是有根树，也可以是无根树，也可以带有枝长信息。codeml命令用于分析各分枝上的omega值，若是有根树的话，对于分析root节点所在的分枝，其结果可能是不对的。因此，作者推荐使用无根树进行正选择分析。无根树是3分叉的树，最少需要3个物种；而3个物种的无根数有3个分枝（2n – 3），对其进行分析，仅能得到某个物种所对应分枝的omega值，这种情况下没太大意义；而使用4个物种的无根树进行分析，可以得到两个物种共同祖先分枝的omega值。故作者推荐使用至少4~5个物种进行分析。一个树文件内容示例：
7 1 ((Hsa_Human, Hla_gibbon) #1,((Cgu/Can_colobus, Pne_langur) #1, Mmu_rhesus), (Ssc_squirrelM, Cja_marmoset));
文件内容分两行：第一行表述树中有5个物种，共计1个树，两个数值之间用任意长度的空分割。此外，Newick格式的树尾部一定要有分号，没有的话程序可能不能正常运行。
2. 输入文件（2/3）：密码子序列比对结果文件 input.txt
7 390 Hsa_Human AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT... Hla_gibbon AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT... Cgu/Can_colobus AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT... Pne_langur AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT... Mmu_rhesus AAGATCTTTGAAAGGTGTGAGTTGGCCAGAACT... Ssc_squirrelM AAGGTCTTCGAAAGGTGTGAGTTGGCCAGAACT... Cja_marmoset AAGGTCTTTGAAAGGTGTGAGTTGGCCAGAACT...
PAML要求输入的Phylip格式，其物种名和后面的序列之间至少间隔两个空格（是为了允许物种名的属名和种名之间有一个空格）。
3. 输入文件（3/3）：codeml命令的配置文件 codeml.ctl
软件难点在于配置文件的理解，常用示例如下：
*输入输出参数：
seqfile = input.txt *设置输入的多序列比对文件路径。
treefile = input.trees *设置输入的树文件路径。
outfile = mlc *设置输出文件路径。
noisy = 9 *设置输出到屏幕上的信息量等级：0，1，2，3，9。
verbose = 0 *设置输出到结果文件中的信息量等级：0，精简模式结果（推荐）；1，输出详细信息，包含碱基序列；2，输出更多信息。
getSE = 0 *设置是否计算并获得各参数的标准误：0，不需要；1，需要。
RateAncestor = 0 *设置是否计算序列中每个位点的替换率：0，不需要；1，需要。当设置成1时，会在结果文件rates中给出各个位点的碱基替换率；同时也会进行祖先序列的重构，在结果文件rst中有体现。
*数据使用说明参数：
runmode = 0 *设置程序获取进化树拓扑结构的方法：0，直接从输入文件中获取进化树拓扑结构（推荐）；1，程序以输入文件中的多分叉树作为起始树，并使用星状分解法搜索最佳树；2，程序直接以星状树作为起始树，并使用星状分解法搜索最佳树；3，使用逐步添加法搜索最佳树；4，使用简约法获取起始树，再进行邻近分枝交换寻找最优树；5，以输入文件中的树作为起始树，再进行邻近分枝交换寻找最优树；-2，对密码子序列进行两两比较并使用ML方法计算DnDs，或对蛋白序列进行两两比较计算ML距离，而不进行其它参数（枝长和omega等）的计算。
* fix_blength = 1 *设置程序处理输入树文件中枝长数据的方法：0，忽略输入树文件中的枝长信息；-1，使用一个随机起始进行进行计算；1，以输入的枝长信息作为初始值进行ML迭代分析，此时需要注意输入的枝长信息是碱基替换率，而不是碱基替换数；2，不需要使用ML迭代计算枝长，直接使用输入的枝长信息，需要注意，若枝长信息和序列信息不吻合可能导致程序崩溃；3，让ML计算出的枝长和输入的枝长呈正比。
seqtype = 1 *设置输入的多序列比对数据的类型：1，密码子数据；2，氨基酸数据；3，输入数据虽然为密码子序列，但先转换为氨基酸序列后再进行分析。
CodonFreq = 2 *设置密码子频率的计算算法：0，表示除三种终止密码子频率为0，其余密码值频率全部为1/61；1，程序分别计算三个密码子位点的四种碱基频率，再算四种碱基频率在三个位点上的算术平均值，计算61种密码子频率时，使用该均值进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；2，程序分别计算三个密码子位点的四种碱基频率，计算61种密码子频率时，使用相应位点的碱基频率进行计算，这61种密码子频率之和不等于1，然后对其数据标准化，使其和为1，得到所有密码子的频率；3，直接使用观测到的各密码子的总的频数/所有密码值的总数，得到所有密码子的频率。一般选择第三种方法，设置CodonFre的值为2。第一种方法让所有密码子频率相等，不符合实际轻卡；第二种方法没有考虑到三种位点各碱基频率的差异，得到的密码子频率不准确；第三种方法能比较准确计算出各密码子的频率；第四种方法由输入数据得到真实的各密码子频率，但有时候有些非终止密码子的频率为0，可能不利于后续的计算。
cleandata = 1 *设置是否移除不明确字符（N、？、W、R和Y等）或含以后gap的列后再进行数据分析：0，不需要，但在序列两两比较的时候，还是会去除后进行比较；1，需要。
* ndata = 1 *设置输入的多序列比对的数据个数。
clock = 0 *设置进化树中各分支上的变异速率是否一致，服从分子钟理论：0，变异速率不一致，不服从分子钟理论（当输入数据中物种相差较远时，各分支变异速率不一致）；1，所有分枝具有相同的变异速率，全局上服从分子钟理论；2，进化树局部符和分子钟理论，程序认为除了指定的分支具有不同的进化速率，其它分枝都具有相同的进化速率，这要求输入的进化树信息中使用#来指定分支；3，对多基因数据进行联合分析。
Mgene = 0 *设置是否有多个基因的多序列比对信息输入，以及多各基因之间的参数是否一致：0，输入的多序列比对文件中仅包含一个基因时，或多个基因具有相同的Kappa和Pi参数；1，输入文件包含多个基因，这些基因之间是相互独立的（这些基因之间具有不同的Kappa和Pi值，且其进化树的枝长也不相关）；2，输入文件包含多个基因，这些基因具有相同的Kappa值，不同的Pi值；3，输入文件包含多个基因，这些基因具有相同的Pi值，不同的Kappa值；4，输入文件包含多个基因，这些基因具有不同的Kappa值和不同的Pi值。当值是2、3或4时，多个基因的进化树枝长虽然长度不一样，但是呈正比关系。这点和参数值等于1是不一样的。
icode = 0 *设置遗传密码。其值1-10和NCBI的1-11遗传密码规则对应：0，表示通用的遗传密码。
Small_Diff = .5e-6 *设置一个很小的值，一般位于1e-8到1e-5之间。推荐检测设置不同的值在比较其结果，需要该参数值对结果没什么影响。
*位点替换模型参数：
model = 1 *若输入数据是密码子序列，该参数用于设置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 = 0 *设置氨基酸之间的距离。
NSsites = 0 *输入数据时密码子序列时生效，用于设置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 = 1 *序列中不同的位点具有不同的碱基替换率，服从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 = 0 *设置一个固定的alpha值或初始alpha值（推荐设置为0.5）。该值小于1，表示只有少数热点位置的替换率较高；该值越小，表示位点替换率在各位点上越不均匀；若设置fix_alpha = 1且alpha = 0则表示所有位点的替换率是恒定一致的。
Malpha = 0 *当输入的多序列比对结果中有多基因时，设置这些基因间的alpha值是否相等：0，分别对每个基因单独计算alpha值；1，所有基因的alpha值保持一致。
ncatG = 5 *序列中不同位点的变异速率服从GAMMA分布的，ncatG是其一个参数，一般设置为5，4，8或10，且序列条数越多，该值设置越大。
*对于密码子序列，当NSites设置为3时，ncatG设置为3；当NSites设置为4时，ncatG设置为5；当NSites值设置>=5时，ncatG值设置为10。
fix_kappa = 0 *设置是否给定一个Kappa值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的Kappa值。
kappa = 2 *设置一个固定的Kappa值，或一个初始的Kappa值。
fix_omega = 0 *设置是否给定一个omega值：0，通过ML迭代来估算Kappa值；1，使用下一个参数设置一个固定的omega值。
omega = .4 *设置一个固定的omega值，或一个初始的omega值。
method = 0 *设置评估枝长的ML迭代算法：0，使用PAML的老算法同时计算所有枝长，在clock = 0下有效；1，PAML新加入的算法，一次对一个枝长进行计算，该算法仅在clock参数值为1，2或3下工作。
4. 运行codeml进行分析
运行程序的命令：
codeml codeml.ctl
5. 正选择基因分析
对基因进行正选择分析，可以按如下步骤进行：
(1) 收集该基因在多个物种中的CDS和protein序列（可以使用orthoMCL分析结果）。 (2) 对该基因的蛋白序列进行多序列比对，再根据CDS序列转换成Codon序列比对结果。这个步骤需要自己编写一些程序来实现。 (3) 根据Codon序列的比对结果构建系统发育树。推荐使用RAxML软件使用ML算法对所有的单拷贝同源基因进行物种树构建。然后，输入物种树信息、Codon多序列比对信息，使用codeml进行omega计算和选择压模型检验。 (4) 快速过滤非正选择基因：设置runmode = -2运行codeml命令对两两序列比较，使用ML方法计算dNdS，若存在omega值大于1，则认为该基因可能属于正选择基因，进行后续分析。 (5) 正选择基因的初步鉴定方法：设置model = 0、NSsites = 1 2运行codeml命令（site models）分别对正选择模型(M2a)和近中性进化模型(M1a)进行检测，若两种模型的似然值相差较大，通过自由度为2的卡方检验，可以确定该基因是否为正选择基因。从程序结果中（例如：lnL(ntime: 11 np: 16): -870.867266 +0.000000）可以找到M2a模型的似然值l1和M1a模型的似然值l0，再计算2Δl = 2(l1-l0)，通过命令“chi2 2 2Δl”根据卡方检验算出p值。此外，也可以设置NSsites = 7 8进行计算，则表示分别对M8(beta & omega)和M7(beta)模型进行检测，和前者类似，若两种模型的似然值相差较大，通过自由度为2的卡方检验，可以确定该基因是否为正选择基因。根据PAML作者所说，M2a和M1a的比较，比M7-M8的比较更严格。即若想得到更多的正选择基因，可以使用M7_VS_M8分析。此外，根据M2和M8模型的BEB（Bayes empirical Bayes）分析结果，还可以得到在多序列比对结果中第一条序列上的正选择位点。根据PAML作者所说，这种site models之间的比较，能用于检测是否在序列上不同的位点具有不同的omega值，而不是用于正选择检测（We suggest that The M0-M3 comparison should be used as a test of variable ω among sites rather than a test of positive selection）。 (6) 指定分化枝上的正选择基因鉴定方法：经过上一步初步鉴定后，设置model = 2、NSsites = 2、fix_omega = 0、omega = 2.0运行codeml命令（branch-site model A）；再设置model = 2、NSsites = 2、fix_omega = 1、omega = 1运行codeml命令（modified branch-site model A / null model）。进行这两种模型分析时，要求输入的树文件中对目标分化枝进行标注。对这两种模型进行LRT分析，计算2Δl = 2(l1-l0)，注意是前者的似然值减后者（null model）的似然值；再使用自由度为1的卡方检验，通过命令“chi2 1 2Δl”计算出的值再除以2，即得到p值。 (7) 目标分化枝上的快速进化基因（Rapidly evolving gene）鉴定方法：设置model = 0、NSsites = 0运行codeml命令，再设置model = 2、NSsites = 0运行codeml命令（这时要求输入的树文件中对目标分化枝进行标注）。然后比较两种运行模式下的似然值，通过自由度为1的卡方检验，可以鉴定该基因在目标分化枝的omega值和背景差异显著。