防杠叠甲：下述所有内容均为本人理解，欢迎友好的学术讨论。觉得看不进去的，写的不好的请windows用户移步右上角×标记，mac用户移步左上角，手机用户请关机。
持续更新中
将依据[Beginner’s Guide on the Use of PAML to Detect Positive Selection]介绍PAML的正选择有关内容
一 不同的密码子替换模型检测适应性进化
- branch model
branch model 假设系统发育树上的不同枝有不同的 ω 值。它们可以用来检测作用于特定谱系的正选择，而不需要对整个系统发生树的 ω 比率进行平均。例如，在检测基因加倍时间后的正选择是useful的，在这种情况下，其中一份拷贝可能获得了新功能，因此有一个更快的进化速率。 -
site model
site model 将基因中任何位点(密码子)的 ω 比率视为来自统计分布的随机变量，因此允许 ω 在密码子之间变化。ω > 1的密码子被定义为正选择。
Likelihood ratio test(LRT)被用来比较零假设（ω < 1）和备择假设（ω > 1）。
可以使用Naïve Empirical Bayes (NEB) method 和Bayes empirical Bayes (BEB)检测哪个密码子受到正选择。
图1 - Branch-site models
Branch-site models 旨在检测预先设置的谱系上的一少部分（only a few sites）的正选择。受到正选择检测的值叫前景枝（foreground branches，ω2 > 1）, 其他的枝叫背景枝（background branches）。背景枝有两类位点，保守的位点（0 < ω0 < 1）和中性的位点（ω1 = 1）。卡方检测鉴定模型选择，BEB检测受选择位点。
下面将会介绍CODEML的四种分析模式：
1）所有位点、所有枝系具有相同的ω
2）检测正选择进化下一部分位点的正选择（site model）
3）检测系统发育树种一个特定的枝或一些特定的枝的正选择（branch model）
4）特定枝一些位点的正选择（branch site model）
Protocol
The Control File
control file包含了COMEML运行需要的参数，所有的可以参考PAML的说明，这里只提及分析有关的参数。
control file的一些规则：
1）一行一行读，如果有相同的两行，第二行的内容会覆盖第一行。
2）一些内容对不同的model都适用。
3）空行会被程序忽略，*之后的内容会被认为是注释。
The Input/Output Files
control file的第一块指定了输入文件的路径：
1)序列比对文件(seqfile)，PHYLIP格式，每一个基因的比对文件都有一个开头(header, 例如"12 1989")，如果一个文件有多个基因比对，需要多个header。
2)系统发育树(treefile)，Newick格式，需要删除bootstrap值，枝长信息最好删除。可能需要一个header（自身体会好像不是刚需），如果有多个树，header是必要的。
3)主要结果储存的输出文件(outfile)
Screen Output
输出到屏幕和输出文件的内容总量是由noisy和verbose控制。
Input Sequence Data
用来推断的数据类型（e.g., nucleo-tides, amino acids, or codons that are to be translated into amino acids by CODEML）是由seqtype控制。
ndata: 基因或者alignments的数量
icode: genetic code
cleandata: if sites with ambiguity data and alignment gaps should be kept (cleandata = 0) or removed (cleandata = 1)
在使用CODEML之前，序列是要对齐的，内含子，非编码区以及终止子要去除。为了保持阅读框架（reading frame），一个有用的策略是首先对齐蛋白质序列，然后相应地构建密码子对齐序列。推荐删除主要是gap的列或者区域，然后使用cleandata=0
Substitution Model
几个参数用于指定分析中使用的模型。第一，ω随着谱系变化（model）,随着位点变化（NSsites）,或者两者都有（model and NSsites）。第二，密码子频率（CodonFreq）是可以相同的（1/61）或者是不同的来定义密码子bias。FmutSel模型（CodonFreq=7）定义每个密码子的fitness都是60（=61-1）。FmutSel0模型（CodonFreq=6）是FmutSel的特殊情况，对于同义密码子的fitness是一样的，因此只使用19（=20-1）个氨基酸fitenss参数。该模型假设氨基酸频率是由蛋白质的功能需求决定的，但同义密码子的相对频率完全由突变偏好参数决定。在这些突变选择模型下，estFreq使用从数据中观察到的密码子频率或最大似然估计密码子频率。取决于用这些参数定义的模型，ω值可以通过fix_omega和omega来指定的或由软件计算。clock指定分子钟(i.e., rate constancy among lineages）。因为我们在这里讨论的模型都是时间可逆的，因此我们不假设分子钟（clock=0），并且系统发育树是无根树（unrooted）。Note，几乎所有的系统发育树构建软件例如RAxML, IQ-TREE和MrBayes都是无根树。
One Ratio With Homogeneous ω Across Lineages and Sites
CODEML中最简单的模型就是M0（one-ratio）模型,其假定多有谱系的所有位点只有一个ω值（图3A）。大多数情况下，这种假设是不现实的，因为大多数位点都被期望在ω < 1的约束下。因此，尝试在M0的假设下寻找正选择的证据是会失败的。事实上，当我们想鉴定正选择时，最好的方式是在后续的描述中使用site或branch-site模型计算ω。然而，M0模型给出了一个零假设，可以作为与其他模型比较的参考，以确定更复杂的模型是否能更好地拟合数据。M0模型下计算出的ω也给出了在位点和谱系上平均的选择性约束的总体度量。以 codeml codeml-M0.ctl运行后，一些内容会打印在屏幕上，一些内容会保存到输出文件（out_M0.txt）内，输出文件可以分为5个部分。
1. 输入比对文件的位点模式总结Summary of the Site Patterns in the Input Alignment
输入的比对序列和相对应的简略比对在输出文件的顶部被输出，每一个位点模式只被呈现一次，相对应的频率（pattren counts）在下方的空白显示。如果使用cleandata=1，删除gap和ambiguous位点前后的比对情况都会被打印。注意，如果gap被删除后，位点会重新编号，这会影响site和branch-site模型下的输出。
2. 输入比对和选择模型总结Summary of the Input Alignment and the Model Selected
CODEML的版本和比对文件的名字会一并输出。后面会是control file里指定模型的具体信息和物种个数以及密码子个数。
3. 核酸和密码子频率的总结Summary of Nucleotide and Codon Frequencies
每个序列中观察到的核苷酸和密码子频率以及它们在所有序列中的平均值被打印在输出文件中，然后是模型下的密码子频率，可以使用诸如evolver之类的软件进行模拟。
4. 树得分和模型参数值计算总结Summary of Tree Scores and Estimated Model Parameter Values
某一模型下的log-likelihood，free parameters的数量，树的枝长和其他的参数在前四行输出。 Free parameters包括枝长，equilibrium frequencies，转换和颠换的比率（κ）和omega分布参数。有n个物种的无根树有2×n-3的枝长，3个mutation-bias parameters （four frequencies with the sum to be 1）。
因为指定的是M0模型(model = 0 and NSsites = 0)，每一个枝有相同的ω 值(0.3357) 。The tree length is dN = 1.1386 at the nonsynonymous sites and dS = 3.3921 at the synonymous sites, with ω = dN/dS = 0.3357，意味着该基因总体受到净化选择，即非同义突变的固定概率(即ω = dN/dS = 0.3357)平均为同义突变的三分之一。
在下一个模块中，我们聚焦三个鉴定正选择的模型，site，branch，branch-site模型
位点ω不一致的Site模型
在这个板块中，在允许ω沿密码子变化的site模型下运行CODEML（图3B）。在这个例子中，运行CODEML基于M0模型（类似之前）和下面的位点不一致模型：M1a,M2a,M7和M8（table1）。编辑配置文件，并重命名为codeml-sites.ctl，并运行CODEML。
输出文件包含控制文件中定义的五个模型中的每个模型的一个部分，其中提供了log-likelihood 和总参数的数量，后者用于确定每个模型比较的自由度。
为了比较nested位点模型，使用LRT来进行检验，其定义为零假设和备择假设的log-likelihood差值的2倍，2Δℓ = 2(ℓ1 − ℓ0), where ℓ0 is the log-likelihood score for the null model, whereas ℓ1 is the log-likelihood under the alternative model。LRT统计量与χ2分布进行比较，其自由度等于两个模型之间自由参数(在每个模型的输出文件中给出)数量的差异。结果总结在表3中。
M0（one-ratio）和M1a（nearly netural）模型可以使用LRT相比。这是对氨基酸位点之间的选择压力变异性的测试，而不是对正选择的测试。相比于M0，M1a更适合示例数据，with 2Δℓ = 559.26，表明ω所反映的选择压力在不同位点之间差异很大。
与M1a相比，M2a增加了一些受到正选择（ ω2 > 1）的位点。这并没有提高模型的适合度因为2Δℓ = 0（table3）。
同时，做了另一个正选择测试通过比较M7(beta, null model)和M8 (beta&ω, alternative model)。M8更适合数据，with 2Δℓ = 12.54，意味着存在带ω >1的正选择位点。阳性选择位点的测试与M1a-M2a和M7-M8比较给出了相互矛盾的结果。当正选择的证据存在但不是很强时，注意到M1a-M2a检验更为严格，因为弱正选择的位点往往被集中到ω1 = 1的位点类别中。
site模型下的参数的最大似然估计（MLE）列在table4。可以使用bash脚本从输出文件中提取这些值。或者，我们在GitHub存储库中的分步教程中包含代码片段。
https://github.com/abacus-gene/paml-tutorial/tree/main/positive-selection/00_ data
尽管M0假定在一个基因中所有的密码子有相同的ω值（model = 0 and NSsites = 0），位点模型假定了几种位点类型（图3A、B）。例如，在M8模型下，94.2%的位点（p0）具有ω的beta分布，然而5.8%的位点的ω = 1.841，表明存在一小部分氨基酸残基在正选择下。这些信息能在输出文件中找到，其位置在MLE开头的下一行。再一次运行中如果指定了多个位点模型，每一个模型在输出文件中有自己的输出信息块，在这里也会显示估计的参数。
在M2a和M8模型下，使用BEB（Bayes Empirical Bayes）方法计算来自不同位点类别的每个位点的后验概率。在正选择类别下具有高后验概率的位点更可能受到正选择。M8模型下的输出如下所示：
在每一行，第一列显示位点的位置（例如，10,25,108，123），后面是氨基酸的类型。第三列（Pr (w > 1)）显示该位点来自正选择类的后验概率。最后一列显示了该位点ω分布的后验均值和标准差。注意，BEB计算仅在正选择模型(即M2a和M8)下进行，而不是在零模型(即M1a或M7)下进行。
在该例子中，14个位点具有大于50%的概率有正选择（ ω > 1），上述只列出4个。例如，第一个序列的位点10是丝氨酸，是正选择类的概率是0.54，ω后验概率的分布的平均值是1.468，标准差是0.634。在这里，BEB的计算只是提供改位点受到正选择的微弱证据，因为ω > 1的后验概率很低并且ω的分布在每个位点是弥漫分布。同时，M1a-M2a的LRT比较不显著并且M7-M8的比较只在5%的水平下显著。总的来说，这些结果显示在这个基因上一些位点受到了正选择，但强度不够强。
沿分枝的不同ω的Branch模型
在这个模块，将展示使用CODEML运行branch模型（图3C），其假定ω沿树的不同分枝变化。分支模型是通过在树文件中使用标记（tags）标记分支来指定的：#0（默认），#1等。至多允许8个不同分枝的ω。在这里我们只考虑两种类型。前景枝（foreground）用#1标记（ ω1），其他的枝作为背景枝且标记为#0（默认，ω0）。注意，#0不用在树文件中标记。可以在Newick树中手动添加前景枝标记，或者使用GitHub中的脚本(https://github.com/abacus-gene/paml-tutorial/tree/main/positive-selection/00_data)，或者使用PhyloTree和EasyCodelML。
这里，我们介绍四种测验，设置了不同的枝作为前景枝：（1）chicken，（2）duck，（3）chicken和duck，（4）chicken和duck以及他们的共同祖先。零假设在所有的测验里都是M0，其假设所有枝ω相同。除了测验4，其他测验均使用无根树。
在第一个测验里，我们想知道chicken是否受到正选择。树文件如下所示：
然后，在其他三种不同的检验下运行相同的分析。配置文件和之前相同，除了以下的部分：
通过使用model=2，我们指定与背景枝相比，前景枝可能在不同的ω下进化。运行其他检验时，注意更改相应的输出文件名。
输出文件格式与之前相同，唯一的区别是输出文件中有一个额外的块，其中包含输入树的每个分支的dN/dS比率。
以chicken为例，输出结果如下：
在此树中，每个分支的估计ω比显示在“#”之后。这里，所有背景枝的ω0 = 0.293，而前景枝的ω1 = 1.171。
注意到以duck为前景枝时，ω1 = 999。999是数值的上限意味着无穷，如果在有关分支上缺乏同义替换，则可能发生这种极端估计。在这种情况下，虽然不能准确的估计ω1，但LRT依然有效。检验3时，ω0 = 0.287，ω1 = 0.711，与检验4一致。结果显示，与其他枝相比chicken和duck有更高的ω值，暗示了可能的正选择。
可以使用LRT检验假设的显著性。根据LRT，我们得出结论，分支模型比M0模型更适合所有测试假设的数据。换言之，不同检验得到的ω值与背景枝的ω比有显著不同。
在这里，我们展示了(1)如何在分枝模型下执行CODEML，对树上不同的进化谱系或分枝假设不同的选择压力(即不同的ω比)，以及(2)如何使用M0作为零模型进行影响预先指定分支的正选择的LRT。我们进行了四个测试，部分是为了说明，如果模型为根周围的两个分支指定了不同的进化过程，则可以使用有根树。注意，测试中前景分支的识别取决于被测试的生物学假设，这应该是先验的。如果在没有任何先验假设的情况下，将树的每个分支依次指定为前景进行分支测试，则需要对多重测试进行校正。
沿分枝和位点的不同ω的Branch-site模型
在本节中，我们在Branch-site模型下运行CODEML，其假设ω在谱系和位点之间都是变化的（图3D）。该模型可用于检测沿预先指定的前景枝影响特定氨基酸位点的正选择。
树文件和branch模型的格式相同，前景枝用#1做标签。我们使用之前branch模型的树文件和按如下修改branch模型的配置文件：
输出结果包括了一个模块，其中列出了branch-site模型A中假设的k = 4个位点类的比例估计值以及背景和前景分支的ω值（table 6）。site class 0 从0-1（0 < ω0 < 1），对于前景枝和背景枝，这一类的位点经历纯化选择。在site class 1，ω值固定为1。这些位点受到中性选择。在site class 2a 或 2b，前景枝经历正选择（ω2 ≥ 1），背景枝受到纯化选择（ 0 < ω0 < 1，site class 2a）或者经历中性进化（ ω1 = 1，site class 2b）。
计算出的ω显示（table 6）在鸡和鸭谱系中存在正选择的位点。事实上，正选择似乎影响了鸟类进化中的所有三个分支。
