关于 PAML 的一二三事

---

PAML 简介

PAML，全称 Phylogenetic Analysis by Maximum Likelihood，由国人大佬杨子恒教授开发，是一款用于进行各种系统发育分析的软件包。里面的诸多模型可以说是系统发育学文章里的常客，其所涵盖的功能包括但不限于：

• 重建系统发育树。
• 检测选择压力。
• 估计进化速率。
• 估计物种或基因的分歧时间。

想要最详细地了解这个软件，可直接下载官方 tutorial 进行啃读：

PAML Manual：http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf

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

也有其他博主写的精细介绍，详见：PAML-discussion-group from 简书

如果有些报错问题或者运行问题在网上找不到解释，也许可以看看 PAML 的谷歌论坛。

提示，本文将基于 linux 系统的操作进行讲解。

更新日志

2023.05.20 补充了支位点模型的局限性。

2023.06.21 补充了自由比模型的相关内容。

2023.06.24 补充了自由比模型的批量运行方法。

2023.07.03 补充了 mcmctree 的运行介绍，精简部分内容。

2023.07.19 补充了 mcmctree 的相关内容，修改了后记内容。

2023.08.19 补充了正选择检测（支位点模型）的内容。

正选择检测（支位点模型）

如果想要系统地了解正选择检测的各种方法及其局限性等相关内容，可以参见官方出版论文教程：初学者的 PAML 正选择检测指南。

该部分侧重于检测特定分支在特定位点上的正选择（因此被称作支位点模型 Branch-Site Model）。其他的还有：

• 支模型（Branch model），检测作用于特定谱系的正选择，m0 vs branch model。
• 位点模型（Site model），检测作用于基因中位点的正选择，可选 M0 vs M1a M1a vs M2a 及 M7 vs M8。

原理及运行

这里首先详细讲一讲它的原理，PAML tutorial 中说法：

The branch-site model A (see the section Codon substitution models above) is specified by using both variables model and NSsites.

Model A: model = 2, NSsites = 2, fix_omega = 0

This is the alternative model for the branch-site test of positive selection. The null model is also the branch-site model A but with ω2 = 1 fixed, specified by

Model A1: model = 2, NSsites = 2, fix_omega = 1, omega = 1

简单来说，正选择检测的原理就是拿出两个模型（一个认为有正选择，一个认为没有）的最大似然对数值，通过似然比检验计算 p 值。

具体说法（可以跳过，无伤大雅）：

正选择检测原理

支位点模型的基本原理是允许树的不同分支和序列的不同部位之间 ω 发生变化，该模型将树的分支分为前景分支（树文件中用特殊符号标记）和背景分支。

在前景分支上，某些位点可能发生正选择，而在背景分支上，这些位点可能没有受到正选择的影响。因此，支位点模型允许我们在指定的前景分支上检测正选择信号：

The branch test (Yang 1998) assigns two different ω parameters (ωF and ωB) to the foreground and background branches on the tree. The null hypothesis of the test is H0: ωF = 1, whereas the alternative is H1: ωF ≥ 1, with ωB to be a free parameter under both hypotheses.

在实际应用中，支位点模型通常包含两个子模型：一个所有分支都中性进化（ω = 1）的零模型（A1）和一个认为前景枝存在正选择（ω > 1）的备择模型（A2）。通过计算两个模型的最大似然值并进行似然比检验，评估正选择的统计显著性。

似然比检验的大致原理是：比较两个模型的最大似然值，从而确定哪个模型更能解释观测数据。步骤大致有：

1. 分别为零模型和备择模型计算最大似然值（给定模型下观测数据的最大概率）。
2. 计算似然比统计量（两模型最大似然值之差的两倍）。
3. 似然比统计量服从卡方分布，其自由度等于两个模型之间的参数差。根据卡方分布的临界值和实际计算得到的似然比统计量，确定拒绝还是接受零假设。

对一个基因进行正选择检测需要运行两次 PAML，并根据输出文件中的最大似然值和自由度进行卡方检验，判断其是否存在显著正选择：

Model A1：

seqfile = 这里要填你的序列文件路径
treefile = 这里要填你的树文件路径，要检测是否有正选择的分支记得标记
outfile = 这里要填你的输出文件，PAML会将结果输出到此处
noisy = 9   
verbose = 0   
runmode = 0  
seqtype = 1   
CodonFreq = 2   
clock = 0
aaDist = 0
model = 2
NSsites = 2   
icode = 0   
Mgene = 0
fix_kappa = 0   
kappa = 2   
fix_omega = 1   
omega = 1   
fix_alpha = 1   
alpha = .0  
Malpha = 0   
ncatG = 3   
getSE = 0   
RateAncestor = 0 
fix_blength = 0  
method = 0   
Small_Diff = .45e-6
cleandata = 1

这些参数的含义（可以跳过，无伤大雅）：

各参数具体含义

noisy ：信息打印级别。

verbose ：输出信息详细程度。

runmode ：设置为 0，PAML 将使用你提供的进化树。

seqtype ：设置为 1，指定核苷酸序列。

CodonFreq ：密码子频率模型，设置为 2 时可以捕获第三个密码子位置的碱基偏好，而这通常与选择压力较低的密码子部位有关。

clock ：分子钟模型，支位点模型不接受其他选择。

aaDist ： 氨基酸替换模型，该处不使用。

model 和 NSsites，设置为 2 代表支位点模型。

icode ：遗传密码表，0 是标准密码表，需根据自身情况修改。

Mgene ：多基因模型，0 表示使用单一模型。针对所有基因使用相同的模型。

fix_kappa：是否固定kappa值，0 表示不固定。这允许程序自动估计转换/颠换比率。

kappa：初始 kappa 值。程序会自动调整 kappa 值以找到最佳估计。

fix_omega：是否固定 omega 值，1 表示固定。

omega：初始 omega 值。

fix_alpha：是否固定 alpha 值，1 表示固定。意味着不考虑不同位点之间的替换速率变异。

alpha：初始 alpha 值。表示不使用 Gamma 分布模型

Malpha：混合 Gamma 分布，0 表示不使用。

ncatG：Gamma 分布类别数。没有使用 Gamma 分布模型，所以类别数设置为 1 即可。

getSE：是否计算参数的标准误差，0 表示不计算。

RateAncestor：是否计算祖先节点的替换速率，0 不计算，如果想顺便重建祖先序列的话也可以调 1，就是慢了一些。注：调 1 方便后续比较祖先序列确定平行替换。

fix_blength：是否固定分支长度，0 表示不固定。这允许程序自动估计分支长度。

method：计算最大似然值的方法，表示使用一阶梯度法。这是 codeml 中常用的方法，适用于大多数情况。

Small_Diff：用于收敛判断的阈值。EasyCodeml 同款配置。

cleandata ：是否删除缺失数据和停止密码子，1 表示要删除，具体如何选择看研究需求。

Model A2：

seqfile = 这里要填你的序列文件路径
treefile = 这里要填你的树文件路径，要检测是否有正选择的分支记得标记
outfile = 这里要填你的输出文件，PAML 会将结果输出到此处
noisy = 9   
verbose = 0   
runmode = 0  
seqtype = 1   
CodonFreq = 2   
clock = 0
aaDist = 0
model = 2
NSsites = 2   
icode = 0   
Mgene = 0
fix_kappa = 0   
kappa = 2   
fix_omega = 0 
omega = 1.5   
fix_alpha = 1   
alpha = .0  
Malpha = 0   
ncatG = 3   
getSE = 0   
RateAncestor = 0 
fix_blength = 0  
method = 0   
Small_Diff = .45e-6
cleandata = 1

相较于 A1，它的变动有：

• fix_omega 变成 0，这允许 ω 自由估计，这时 omega 的值将作为优化过程的起点。
• omega 变成 1.5，表示在优化过程开始时，ω 的初始值为 1.5。

在不同的 paper 中或许会看到不同的 omega 值（常是 1.5 或 2），但这一般无关紧要，因为在后续的优化过程中，该值会变化以最大化似然。这个值只是一个起始点，并不代表最终 ω。选择一个合适的初始值只是有助于加速优化过程并提高收敛的稳定性。但大多数情况下，不同的初始值通常会收敛到相似的结果。

运行（上面两个模型分别操作，控制文件中就是上面两个模型的具体参数）：

$ codeml xxx(你的控制文件)

不同模型跑出的结果（outfile）的似然值一般都是这种样子：

lnL(ntime:  7  np: 11):  -1225.318164

查找 lnl 就能找到，两个文件里都有，记录下两个文件中的 np 和 lnL （上例中为 11 和 -1225.318164）。计算 D 和自由度：

• D = 2 | lnL1 - lnL0 |
• df = np1 - np0

关于此处有个问题，从数学角度出发，使用似然比检验时更高的似然值意味着模型更能适应或描述数据，因此正选择模型（Model A2）的似然值显著高于中性模型（Model A1）时才应该能够成为数据中存在正选择的证据。

但是不管是网上的各个教程还是 EasyCodeml 中，对于似然值之差都取绝对值进行似然比检验，这导致有些时候正选择模型的似然值就算更低也能得到统计显著性（这种情况虽然少但确实会发生），这又给我带来了两个疑问：

• Model A2 包含更多的参数，这提供了更多的自由度来捕捉数据中的模式，因此它的似然值理应更高（或至少一样）。
• 为什么两个模型存在显著差异（不管哪个拟合的更好）就能说明存在显著正选择？

目前我并不确定具体结论，现已在 PAML discussion group 上发贴求助，等得到具体回答以后会对文章内容进行修改。

然后使用卡方检验，这里可以用 python scipy 包中的 scipy.stats.chi2 进行。

以上就是整个正选择检测的流程，接下来放一个脚本，该脚本可以用于批量进行正选择检测，该脚本原作者是一个简书博主，这里对他的脚本进行了一定的更改以适应合理的批量处理：

# callCodeml.py
import re
import sys
import os
from scipy.stats import chi2
import time
from typing import List, Tuple

def help_info() -> None:
    print("This script is used to call codeML.\nUsage: python3 callCodeml.py Dir treeFile cpuNum\n")

def get_args() -> Tuple[str, List[str], str, int]:
    args = sys.argv[1:]
    if len(args) < 3:
        help_info()
        sys.exit(1)
    
    dir_path = os.path.abspath(args[0])
    if not os.path.isdir(dir_path):
        help_info()
        sys.exit(1)
    
    tree_file = os.path.abspath(args[1])
    try:
        cpu_num = int(args[2])
    except ValueError:
        help_info()
        sys.exit(1)
    
    seqs = [x for x in os.listdir(dir_path) if not x.startswith('.')]
    return dir_path, seqs, tree_file, cpu_num

def create_dir(base_dir: str, name: str) -> None:
    dirs = [
        f'{base_dir}/',
        f'{base_dir}/{name}',
        f'{base_dir}/{name}/null',
        f'{base_dir}/{name}/alte'
    ]
    for dir_path in dirs:
        os.makedirs(dir_path, exist_ok=True)

def create_ctl(base_dir: str, dir_path: str, name: str, tree_file: str) -> None:
    ctl_null = f'''
    seqfile = {dir_path}/{name}
    treefile = {tree_file}
    outfile = {base_dir}/{name}/null/{name}_null.res
    noisy = 9
    verbose = 0
    runmode = 0
    seqtype = 1
    CodonFreq = 2
    clock = 0
    aaDist = 0
    model = 2
    NSsites = 2
    icode = 0
    Mgene = 0
    fix_kappa = 0
    kappa = 2
    fix_omega = 1
    omega = 1
    fix_alpha = 1
    alpha = .0
    Malpha = 0
    ncatG = 3
    getSE = 0
    RateAncestor = 0
    fix_blength = 0
    method = 0
    Small_Diff = .45e-6
    cleandata = 1
    '''
    ctl_alte = ctl_null.replace(f'{base_dir}/{name}/null/{name}_null.res', f'{base_dir}/{name}/alte/{name}_alte.res')\
                       .replace('fix_omega = 1', 'fix_omega = 0')\
                       .replace('omega = 1', 'omega = 1.5')
    
    with open(f'{base_dir}/{name}/null/null_profile.ctl', 'w') as f:
        f.write(ctl_null)
    with open(f'{base_dir}/{name}/alte/alte_profile.ctl', 'w') as f:
        f.write(ctl_alte)

def run_create(base_dir: str, dir_path: str, seqs: List[str], tree_file: str) -> None:
    for name in seqs:
        create_dir(base_dir, name)
        create_ctl(base_dir, dir_path, name, tree_file)
    
    with open(f'{base_dir}/codeml.sh', 'w') as f:
        f.write('''
#!/bin/bash
list=$(ls -1)
count=$(ls -1 | wc -l)
echo "Total number: $count"
dir=$(pwd)
echo "Start Now..."
for i in $list
do
    {
    echo "cd $dir/$i/null; codeml ./null_profile.ctl > log.txt; cd $dir/$i/alte; codeml ./alte_profile.ctl > log.txt" >> paml_command
    }
done
wait
''')
    print(f"{time.strftime('%Y-%m-%d %H:%M:%S')} Files created.\n")

def run_codeml(base_dir: str, num: int) -> None:
    os.chdir(base_dir)
    os.system('bash codeml.sh')
    os.system(f'cat paml_command | parallel --no-notice -j {num}')
    print(f"\n{time.strftime('%Y-%m-%d %H:%M:%S')} Codeml Done...\n")

def get_res(path: str) -> Tuple[float, int, float]:
    with open(path, "r") as f:
        t = f.read()
        kappa = float(re.findall(r"kappa\s*=\s*([\d\.]+)", t)[0])
        lnL = float(re.findall(r'lnL\s*\(.*\)\s*=\s*([-\d\.]+)', t)[0])
        np = int(re.findall(r'lnL\s*\(.*\)\s*=\s*[-\d\.]+\s*\(\s*(\d+)', t)[0])
        return lnL, np, kappa

def run_stat(base_dir: str, name: str) -> None:
    try:
        path_alte = f'{base_dir}/{name}/alte/{name}_alte.res'
        path_null = f'{base_dir}/{name}/null/{name}_null.res'
        
        lnL0, np0, _ = get_res(path_null)
        lnL1, np1, kappa = get_res(path_alte)
        
        lnl = abs(lnL0 - lnL1) * 2
        np = abs(np0 - np1)
        pvalue = 1 - chi2.cdf(lnl, np)
        
        with open(f'{base_dir}/result.tsv', 'a') as f:
            f.write(f"{name}\t{lnL0}\t{lnL1}\t{np0}\t{np1}\t{kappa}\t{pvalue}\n")
    except Exception as e:
        print(f"{time.strftime('%Y-%m-%d %H:%M:%S')} {name} failed to stats: {e}")

def main() -> None:
    start_time = time.strftime('%Y%m%d-%H%M%S')
    base_dir = f"{os.getcwd()}/WorkingDir_{start_time}"
    
    dir_path, seqs, tree_file, cpu_num = get_args()
    run_create(base_dir, dir_path, seqs, tree_file)
    run_codeml(base_dir, cpu_num)
    
    for name in seqs:
        run_stat(base_dir, name)
    
    print(f'Total time used: {time.time() - start_time:.2f} seconds\n')

if __name__ == "__main__":
    t0 = time.time()
    main()
    print(f'Total time used: {time.time() - t0:.2f} seconds\n')

相较于原来的版本，变动有：

• 使用了 parallel 来控制并行命令的数量，因为一开始会对所有的序列文件同时开始处理，但如果用到了很多（比如几千个）序列的话这样做就会把 CPU 挤爆，所以加上一些控制的手段保证能几十个几十个来运行。
• 把文件名识别去除了（之前貌似有限制以特定字符开头或结尾啥的），过程中可能会有部分报错但不会影响结果。

运行方式：

$ python callCodeml.py Dir treeFile numofcpu

• Dir 指定存放 PAML 格式序列文件的文件夹。
• treefile 指定物种树。
• numofcpu 指定预计使用的线程数。

所有结果文件会输出在当前目录下以 WorkingDir 开头的文件夹中，所有结果的统计将出现在该文件夹下的 result.tsv 中。

运行前可能需要注意的事项：

• 输入的树文件要为无根树，输入有根树会导致 PAML 估计错误。
• 输入的序列文件要是 PAML 需要的格式（非 fasta）。
• 确保前景枝已被标记，具体的标记方法可以看之前推荐的教程。

你可能会感兴趣的问题：

• 正选择检测时，应当输入每个基因的基因树作为树文件还是输入物种树作为树文件这一个问题，事实上在很多论文里都是有分歧的，两种方法都有人做，不过也有文章已经指出过，当基因树和物种树拓扑结构不一致时，将物种树作为输入树将会错误地估计枝长，导致更高的假阳性率。

Gene tree discordance causes apparent substitution rate variation.

from Systematic Biology

Discussion group 里也有人提过这个问题，杨子恒教授的说法：

The best tree to use would be one that accurately describes the true evolutionary history of the aligned sequences.

If the gene and species trees disagree, you’ll need to decide which is a better approximation of the true history.

– a gene tree may be a poor indicator of true history simply because of limited sample size, or because of unusual patterns of sequence evolution (say, if there was strong convergent evolution)

– a species tree may be a poor indicator of true history if gene duplications occurred or if speciation events were tightly clustered (say, if there was incomplete lineage sorting)

Regardless of which tree you select, reviewers may ask that you conduct some/all of the analyses again on the other tree to ensure the results are robust.

• 经过卡方检验发现部分基因存在显著正选择，但是并没有正选择位点。这可能是因为整个基因表现出了正选择但是单独位点的信号太弱。

Can we say a gene is under positve selcetion if the LRT tests results is significant but BEB detects no positively selective sites?

• 是否应当去除外群进行分析？杨子恒教授认为，如果外群离得不是很远，那么将其纳入正选择检测可能是有益的。

Question on branch labeling with two species (multiple strains)

• 不同版本的 PAML 跑出来的结果可能会不同，请做好版本控制工作。
• 是否要删除比对中的 gap ？杨子恒团队建议删除主要是 gap 或难以对齐的区域之后，使用 cleandata = 0 保留数据中的信息。

还有很多很好的问题有博主整理出来了，可见 PAML-discussion-group from 简书。这里我就只放一些我在分析时遇到的几个问题还有对应的解释出来。

局限性

虽然说支位点模型是最广为应用的正选择检测方法之一，但是它的很多假设也在一定程度上削弱了它的可靠性。

• 支位点模型假设所有分支所有位点的进化是独立的，在很多情况下这并不成立（你可能想知道 遗传连锁 或 协同进化）。
• 支位点模型先验地将所有分支分为了存在正选择的前景枝和负选择或中性选择的背景枝两类，当背景枝的情况不符合这个假设时会导致非常高的假阳性率或假阴性率（Sergei et al.）。
• 支位点模型假设所有的分支在所有的位点中共享相同的 dS，在很多情况下这并不成立。

此外，也要考虑到一些可能存在的客观不足，比如说：

• 样本量过小（小样本偏差）。如果分析中涉及到的物种太少，那么物种的特异性就会被过度放大（考虑少数物种时以为存在正选择信号，但是加入更多物种后发现实际上并没有），而且参数估计可能会有很大的不确定性，导致结果偏向某一方向。
• 某些正选择信号可能是遗传漂变造成的（轻微有害突变被漂变固定），这在有效种群大小较低的物种中更有可能发生。

所以，如果想要得到更准确的正选择结果，可以考虑使用 Hyphy 等其他生信软件（或模型）对结果进行补充，同时注意采取严格的假阳性控制手段（例如 BH 矫正等）。

Hyphy 的模型有时候会表现得更好，现在越来越多的相关文章都是同时使用 PAML 和 Hyphy 进行分析，然后将两者的结果一起展示。

dN / dS 计算（自由比模型）

原理及运行

这一部分也可以叫做选择压力分析，使用的是 PAML Branch Model 中的 free-ratio model。

与之前讲到的支位点模型不同，自由比模型允许 ω(dN/dS) 在所有分支中变化，它对所有分支的 ω 进行独立统计（包括末端分支和祖先分支）。

PAML 的控制文件中，使用 model = 1 来指定自由比模型。

具体如下：

seqfile = 这里要填你的序列文件路径
treefile = 这里要填你的树文件路径，要检测是否有正选择的分支记得标记
outfile = 这里要填你的输出文件，PAML 会将结果输出到此处

noisy = 9
verbose = 1
runmode = 0

seqtype = 1
CodonFreq = 2
model = 1
NSsites = 0
icode = 0
Mgene = 0

fix_kappa = 0
kappa = 2

fix_omega = 0
omega = 1

fix_alpha = 1
ncatG = 8

getSE = 0
RateAncestor = 0
Small_Diff = .5e-6
cleandata = 1

部分参数并没有详细指定，PAML 会选择默认选项（默认情况可以参见 pamlDOC.pdf）

然后直接运行即可：

$ codeml xxx(你的控制文件)

得到的结果可在 outfile 中查看，一般而言各个分支的 ω 值都在最后一行（以井号 # 作为标识），就像这样：

(A #0.1 , B #0.2 , C #0.3)

和支位点模型相同，这里列出批量运行自由比模型的 Python 脚本：

# -*- coding: utf-8 -*-
# free-ratio-calcu
# Author: Juse

import sys
import os

args = sys.argv[1:]
if len(args) == 0:
	print("Usage: python free-ratio_calcu.py treefile seqdir outputDir")
	sys.exit()
try:
	treefile = os.path.abspath(args[0])
	seqdir = os.path.abspath(args[1])
	outputDir = os.path.abspath(args[2])
except:
	print("Please check your command.")
	sys.exit()

paml_seq = []
for paml in os.listdir(seqdir):
	paml_seq.append(f"{seqdir}/{paml}")

for seqid in paml_seq:
	outid = f"{seqid.split('/')[-1]}.mlc"
	with open("codeml.ctl", "w") as f:
		f.write(f'''
seqfile = {seqid}
treefile = {treefile}
outfile = {outputDir}/{outid}

noisy = 9
verbose = 1
runmode = 0

seqtype = 1
CodonFreq = 2
model = 1
NSsites = 0
icode = 0
Mgene = 0

fix_kappa = 0
kappa = 2

fix_omega = 0
omega = 1

fix_alpha = 1
ncatG = 8

getSE = 0
RateAncestor = 0
Small_Diff = .5e-6
cleandata = 1
			''')
	os.system("codeml")

具体脚本已经上传到 GitHub 上（链接），运行：

$ python free-ratio_calcu.py treefile seqdir outputDir

• treefile 指定物种树文件。
• seqdir 指定序列文件夹，应将所有要运行自由比模型的序列放置在其中并且不包括任何其他文件，序列应为 PAML 格式。
• outputDir 指定输出文件夹，所有的结果文件将在该处产生，建议指定一个空文件夹以衔接后续批量提取。

在这里我并没有像支位点模型那个脚本一样设置多线程操作，所以如果数据较多的友友可以尝试自行修改脚本以提高效率。此外序列的格式要为 PAML 所需格式而非传统 fasta 格式，需要进行序列格式转换的可以使用 Phylosuite 或者 JuseKit。

批量提取结果

一般来说我们关注的都是末端分支 ω 值，然后取出进行统计比较，如果想要快捷提取出不同物种的 ω 值可以考虑使用 shell 的 grep 或者 Python 的 re。

在这里给出我的 Python 脚本实例（下载链接）：

# free-ratio-omega
# Author: Juse

import sys
import re
import os

spefile = sys.argv[1]
mlcfiledir = sys.argv[2]


species = []
with open(spefile, 'r') as f:
    for line in f:
        if line.strip() != '':
             species.append(line.strip())
species[:] = sorted(species)

with open('free-ratio.result', 'w') as op:

    head = '\t'.join(species)
    op.write(f'File\t{head}\n')

    for mlc in os.listdir(mlcfiledir):
        op.write(mlc+'\t')
        mlcpath = f'{mlcfiledir}/{mlc}'
        with open(mlcpath, 'r') as m:
            file_content = m.read()
            for spe in species:
                search = fr'{spe} #(\d+(?:\.\d+)?)'
                omega = re.search(search, file_content).group(1)
                op.write(omega + '\t')
        op.write('\n')

要求的输入有两个：

• 物种名文件，格式如下（一行一个物种名，填需要统计 ω 值的物种）：

SpeA
SpeB
SpeC

• PAML 运行结果文件（也就是上面的 outfile）。

脚本运行方式：

$ python free-ratio-omega.py spename_file path_of_outfile

把 spename_file path_of_outfile 两项换成对应文件和路径即可。

统计结果会出现在运行命令的路径中（free-ratio.result），结果文件实例：

File	A	B	C
1	0.00733407	0.0001	0.00574998
2	0.119151	0.0917734	0.0648513
3	0.0653332	0.0288158	0.0296085

需要注意的点：

• 在一些基因中，可能会存在有些分支几乎没有变化的情况，这时很容易出现极端值（999 或 0.0001），因此在分析之前最好排除这些基因。关于这一点，部分文献中也有讲到一定的过滤方法，例如根据 N*dN S*dS 的大小进行过滤。
• 本文章提供的脚本仅支持获取末端分支的 ω 值，如有需要提取某一祖先分支 ω 值，可以结合 Python 的 ete3 包完善脚本进行。

然后就可以使用 R 或其他软件根据得到的结果进行绘图：

估计物种分歧时间

这一功能通过 PAML 的 mcmctree 实现，虽然该方法较本文中已经描述的几种方法要更复杂得多（涉及到许多参数的设置），但是由于其拥有专门的 tutorial 教程，因此这里我只挑出其中一些我个人认为比较重要的地方进行展示。

首先，请过目 mcmctree 的 tutorial，这对于理解程序的运行和各个参数而言还是非常重要的，其中还展示了相关示例和实践，可以跟着做一遍。此外，也可以看看简书博主总结的教程：mcmctree估算物种分歧时间。

如果不嫌看英文麻烦，也可以看看这篇教程：Example of mcmctree analysis

如果想系统地了解这部分的理论背景，可以看这篇综述文章：https://www.nature.com/articles/nrg.2015.8

这里讲一些比较关键的点：

首先是估计物种分歧时间所需要的条件：一个有根树（无枝长信息）以及相关的化石矫准。也可以用其他文献中估计的分歧时间作为校准点（这种情况被称为二次校准），不过相对来说没那么严谨；或者说以某个地质历史事件发生时间作为校准点，例如某个岛屿分裂导致地理隔离并随之引起两物种的分化。在该分析中，校准信息越准确，分析结果越可靠。

参数设置的相关注意事项：

• 在松弛钟模型（clock = 2 or 3）下，如果树根部没有校准点，那么就需要通过 RootAge 指定一个分歧时间上限（所有物种的最近共同祖先形成时间）。

Juse 注：松弛钟模型允许不同分支进化速率不同，这与严格钟模型相反。

• 选择时间单位时应使节点的分歧时间大致在 0.01-10 的范围内（超过 10 应该也没有关系但不要低于 0.01），比方说如果分歧时间在大约 100-1000MY，那么可以将 100MY 视作一个时间单位。此外，一些参数的设置需根据你对时间尺度的选择进行调整。

Juse 注：具体的例子在 mcmctree tutorial 的 11 页可见。

在 clock = 2 时，只有 rgene_gamma 要随之变动。而 clock = 3 时，sigma2_gamma 也要变化。

• 指定先验替换率，具体怎么设置取决于对替换率有多大的把握，如果并无太大把握就将 αμ 设置成 1 1.5 或 2，然后调整 βμ 至其合理。

Juse 注：有些生物类群的替换率可能能从文献中得到，例如：例子。

虽然 PAML 现在使用的 Dirichlet 先验能一定程度上减小错误替代率带来的影响，但关于替换率的参数设定仍需要谨慎，因为该值会直接影响分歧时间的后验估计。

如果没有相关信息，可以通过指定 clock=1 运行 baseml 来获得大致的替换率估计并用于先验设定（详见之前提供的英文教程或后文补充部分）。

• 最好进行多次分析，确保结果相似（要使用不同的 seed 或直接将 seed 设置为 -1）。对于分析是否已经收敛，可以用 tutorial 中提到的方法去判断，也可以看一看这篇文章中阐述的方法（基于 R 实现）。若结果不相似，说明分析可能还未收敛，此时可能需要对 burnin sampfreq nsample 三个参数进行一定的修改。
• 此处涉及了一个新的参数为 ndata，指定的是序列的分区数量，tutorial 中所涉及到的序列就为分区数据（分三个部分，分别是密码子的第一、第二和第三个碱基），如果没有分区的话指定 ndata = 1 即可。
• 关于选择软边界还是硬边界，文章指出使用软边界能够更加灵活并且减弱错误化石校准的影响。
• 关于出生率、死亡率和抽样分数的先验设置，文章发现其对于后验时间估计只有很小的影响，但依然建议尝试改变这些值以观察后验估计的鲁棒性。

Bayesian Estimation of Species Divergence Times Under a Molecular Clock Using Multiple Fossil Calibrations with Soft Bounds

from Molecular Biology and Evolution

• 关于怎么选择模型，你可能会对以下信息感兴趣：

不同 Model 指定的具体模型：

0:JC69，1:K80，2:F81，3:F84，4:HKY85

ChatGPT3.5 的回答

在这些模型中，JC69（Jukes-Cantor 69）是最简单的模型，它假设所有碱基或氨基酸的变异率相等，并且没有考虑碱基或氨基酸之间的相似性。因此，JC69模型通常只用于较近的亲缘关系分析。

K80（Kimura 2参数）模型是比JC69更复杂的模型，它考虑了不同碱基或氨基酸之间的相似性，并引入了转换和转换/颠换的不同变异率。

F81（Felsenstein 81）模型是另一个比JC69更复杂的模型，它考虑了碱基或氨基酸之间的不同变异率，并假设平稳状态的分布是相等的。

F84（Felsenstein 84）模型是F81模型的改进版本，它引入了不同的碱基或氨基酸转换和转换/颠换的变异率。

而HKY85（Hasegawa-Kishino-Yano 85）模型是这些模型中最复杂的一个。它在F84模型的基础上引入了碱基或氨基酸的酶位点偏好，并考虑了碱基或氨基酸的相对频率。因此，HKY85模型在描述分子演化过程中更准确和灵活，具有更高的泛用性。

ChatGPT4.0 的回答

1. **JC69 (Jukes-Cantor 1969)**：是最简单的模型之一。它假定四种碱基（A, T, C, G）替换的概率是均等的，即任何两种碱基之间的替换率是相等的。
2. **K80 (Kimura 2-parameter)**：比JC69稍微复杂一些。它引入了两个参数，一个是转换的替换率（purine之间或pyrimidine之间的替换），另一个是转变的替换率（purine和pyrimidine之间的替换）。
3. **F81 (Felsenstein 1981)**：在JC69的基础上增加了考虑站点频率不均的特点，即每种碱基出现的概率可以不同。
4. **HKY85 (Hasegawa, Kishino and Yano 1985)**：这个模型进一步扩展了K80，允许碱基频率不等，并且区分了转换和转变的替换率。
5. **F84 (Felsenstein 1984)**：类似于HKY85，但在处理转换和转变的替换率上略有不同。

综合来看，HKY85和F84是你所列举的模型中较为复杂的，它们在处理碱基替换时考虑了更多的参数。由于HKY85和F84在处理替换率和频率方面都较为灵活，因此它们通常具有较高的泛用性。在许多情况下，HKY85模型是一种非常受欢迎的选择。

在文章的具体示例中，分别使用了 JC 模型和 HKY+G 模型进行对比，JC 在部分情况下会显得不可靠。

所以如果不缺计算资源和时间的话， Model = 4 或许会是一个万金油选择（同时设置 alpha = 0.5 & ncatG = 5 以使用 HKY85+G5 Model）。

实战演练

由于该部分 tutorial 已经提供了相关的示例，因此这里为了我将直接示范其通过近似似然节省运行时间的方法。

首先，为什么要使用近似似然 —— 因为近似似然能大幅节约时间，这是 tutorial 的原文：

For large alignments, calculation of the likelihood function during the MCMC is computationally expensive, and estimation of divergence times is very slow. Thorne et al.[8] suggested using an approximate method to calculate the likelihood that improves the speed of the MCMC dramatically.

这是论文片段：

Because Bayesian molecular-clock dating relies on expensive MCMC sampling, the time savings obtained by using the approximation can be quite dramatic, up to 1000× depending on the dataset.

本次示例使用的文件均在 PAML 根目录的 examples/DatingSoftBound 下，且均在 tutorial 中可见。

首先跳转到该文件夹，该文件夹中有一控制文件为 mcmctree.ctl，将其内容改为：

        seed = -1
     seqfile = mtCDNApri123.txt
    treefile = mtCDNApri.trees
    mcmcfile = mcmc.txt
     outfile = out.txt

       ndata = 3
     seqtype = 0    * 0: nucleotides; 1:codons; 2:AAs
     usedata = 3    * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV)
       clock = 2    * 1: global clock; 2: independent rates; 3: correlated rates
     RootAge = '<1.0'  * safe constraint on root age, used if no fossil for root.

       model = 4    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
       alpha = 0.5    * alpha for gamma rates at sites
       ncatG = 5    * No. categories in discrete gamma

   cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

     BDparas = 1 1 0.1  * birth, death, sampling
 kappa_gamma = 6 2      * gamma prior for kappa
 alpha_gamma = 1 1      * gamma prior for alpha

 rgene_gamma = 2 20 1   * gammaDir prior for rate for genes
sigma2_gamma = 1 10 1   * gammaDir prior for sigma^2     (for clock=2 or 3)

    finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0 or 1): times, musigma2, rates, mixing, paras, FossilErr

       print = 1   * 0: no mcmc sample; 1: everything except branch rates 2: everything
      burnin = 2000
    sampfreq = 10
     nsample = 20000

$ mcmctree mcmctree.ctl

可以在相关文件中查看详细的序列和物种树校准信息，运行上述命令后，文件夹中会出现一个名为 out.BV 的文件，将其重命名为 in.BV。

$ mv out.BV in.BV

随后修改控制文件，将 usedata = 3 换为 usedata = 2 ，再次运行并得到结果。

比对一下结果，使用近似似然时（用时 10 s）：

t_n8           0.1959 ( 0.1627,  0.2477) ( 0.1588,  0.2391)  0.0803  (Jnode 12)
t_n9           0.1571 ( 0.1464,  0.1629) ( 0.1479,  0.1639)  0.0160  (Jnode 11)
t_n10          0.0928 ( 0.0814,  0.1074) ( 0.0804,  0.1062)  0.0258  (Jnode 10)
t_n11          0.0640 ( 0.0593,  0.0739) ( 0.0587,  0.0723)  0.0136  (Jnode  9)
t_n12          0.0256 ( 0.0191,  0.0336) ( 0.0188,  0.0331)  0.0143  (Jnode  8)
t_n13          0.0484 ( 0.0378,  0.0607) ( 0.0374,  0.0602)  0.0228  (Jnode  7)
mu1            0.4790 ( 0.3951,  0.5705) ( 0.3933,  0.5675)  0.1742
mu2            0.1643 ( 0.1270,  0.2092) ( 0.1262,  0.2078)  0.0816
mu3            2.9379 ( 2.0577,  3.7933) ( 2.0485,  3.7801)  1.7316
sigma2_1       0.0436 ( 0.0022,  0.1503) ( 0.0001,  0.1206)  0.1205
sigma2_2       0.0941 ( 0.0103,  0.2748) ( 0.0001,  0.2292)  0.2291
sigma2_3       0.2195 ( 0.0713,  0.5266) ( 0.0517,  0.4588)  0.4070
lnL          -16.1493 (-24.6090, -9.3610) (-23.9400, -8.9710) 14.9690

不使用近似似然时（usedata = 1，用时 24 min 38 s）：

t_n8           0.1931 ( 0.1637,  0.2304) ( 0.1609,  0.2260)  0.0651  (Jnode 12)
t_n9           0.1572 ( 0.1470,  0.1630) ( 0.1483,  0.1637)  0.0154  (Jnode 11)
t_n10          0.0912 ( 0.0815,  0.1036) ( 0.0809,  0.1025)  0.0217  (Jnode 10)
t_n11          0.0631 ( 0.0592,  0.0713) ( 0.0585,  0.0699)  0.0114  (Jnode  9)
t_n12          0.0243 ( 0.0191,  0.0311) ( 0.0188,  0.0306)  0.0119  (Jnode  8)
t_n13          0.0464 ( 0.0375,  0.0572) ( 0.0370,  0.0566)  0.0195  (Jnode  7)
mu1            0.4824 ( 0.4064,  0.5676) ( 0.4033,  0.5633)  0.1600
mu2            0.1677 ( 0.1323,  0.2090) ( 0.1300,  0.2061)  0.0761
mu3            2.7935 ( 2.2004,  3.3122) ( 2.2439,  3.3458)  1.1019
sigma2_1       0.0340 ( 0.0022,  0.1169) ( 0.0001,  0.0933)  0.0933
sigma2_2       0.0712 ( 0.0049,  0.2187) ( 0.0001,  0.1797)  0.1796
sigma2_3       0.0717 ( 0.0109,  0.2105) ( 0.0015,  0.1746)  0.1731
kappa_1       12.3972 (10.4875, 14.6043) (10.3315, 14.4080)  4.0765
kappa_2        9.3277 ( 7.2435, 11.8446) ( 7.1484, 11.6949)  4.5464
kappa_3       33.8883 (30.1020, 38.0932) (30.0334, 37.9816)  7.9483
alpha_1        0.2192 ( 0.1746,  0.2733) ( 0.1717,  0.2689)  0.0972
alpha_2        0.0436 ( 0.0019,  0.1100) ( 0.0000,  0.0987)  0.0987
alpha_3        3.9802 ( 2.7839,  5.7941) ( 2.6174,  5.5002)  2.8828
lnL        -29997.5335 (-30011.0870, -29985.7710) (-30010.5710, -29985.3530) 25.2180

对比不使用近似似然和使用近似似然时的结果可以发现，两次分析所得到的各节点分歧时间都是非常相近的，但是不使用近似似然分析所用的时间却花费了 140 倍以上。

此外，对于蛋白和密码子序列数据，若要使用近似似然方法节约时间可能需要进行额外的操作，详情可见 Tutorial 4: Approximate likelihood with protein data。

要补充的一些话

本文中关于 mcmctree 的部分有很多相关的结果信息并没有详细说明，这主要是因为这一部分与之前讲的两个部分有个根本的不同 —— 它具有一个单独的 tutorial 并且 PAML DOC 中也有专门讲解，这些内容基本上涵盖了所有我未提及到的其他地方。

俗话说得好，好记性不如烂笔头。跟着看跟着学跟着做才是最好的学习方法，做出来了一遍还怕后面做不成吗？

估计替换率的方法

详尽版：A Step-by-Step Tutorial: Divergence Time Estimation with Approximate Likelihood Calculation Using MCMCTREE in PAML

创建一个 baseml.ctl 文件，写入以下内容：

seqfile = xxx
treefile = xxx
outfile = mlb

noisy = 3
verbose = 1
runmode = 0

model = 7
Mgene = 0

fix_kappa = 0
kappa = 2

fix_alpha = 0
alpha = 0.5
Malpha = 0
ncatG = 5

fix_rho = 1
rho = 0.
nparK = 0

clock = 1
nhomo = 1
getSE = 1
RateAncestor = 0
cleandata = 0

其中 seq 为序列文件（PAML 或 Phy 格式），treefile 为根处有校准的有根树文件，例子：

7  1
((((human, (chimpanzee, bonobo)), gorilla), (orangutan, sumatran)), gibbon)'@1.0';

在准备好相关文件后，运行以下命令：

$ baseml baseml.ctl

输出的 mlc 文件中存在以下内容：

Substitution rate is per time unit

    0.090798 +- 0.006418

根据该值，假设 m = s = 0.09，则：

• a = (0.09/0.09)^2 = 1
• b = 0.09/0.09^2 = 11.1

因此可在运行 mcmctree 时将 rgene_gamma 设置为：rgene_gamma = 1 11.1

---

除了已介绍的模型，还有一些模型例如 Two-ratio model 和 Fixed ratio model 等没有讲到，不过这些模型的运行只需改几个参数，包括它们之间的比较（比如说 Two-ratio vs. Fixed）同样是通过似然比检验，因此后面不会再特别补充，以后的变动可能主要是对现有的内容进行完善（不过也说不定）。