EasyCodeML在适应性进化分析中的应用

研究背景

在DNA水平上检测正选择，并找出对应的氨基酸位点是适应性进化分析中不可或缺的一项重要内容。在基于密码子 (Codon-based) 进化模型的生物软件和资源中，PAML包中的CodeML 程序目前使用最为广泛 (Yang, 2007)。根据序列中不同位点或系统发育树中不同分支上ω 值 (非同义替代与同义替换的比值，即dN/dS或Ka/Ks) 不同，CodeML中的密码子模型可以分为枝模型 (Branch model)、位点模型 (Site model)、枝位点模型 (Branch-site model) 和进化枝模型 (Clade model) 四大类(Yang and Nielsen, 1998 and 2002; Forsberg and Christiansen, 2003; Bielawski and Yang, 2004; Zhang et al., 2005)。
    位点模型 (Site model) 常用于检测核苷酸序列中位点的正选择 ( ω >1) 作用，该模型假设数据集中不同密码子的选择压不同，而支系 (Lineage) 间的选择压力相同。位点模型中共有8个不同假设的模型 (表1)，其中M0和M3、M1a和M2a、M7和M8以及M8a和M8为四对巢式模型 (Nested model)，M0、M1a、M7和M8a为零假设模型 (null model)，M3、M2a和M8为备选假设模型 (alternative model)。除M0 vs. M3 (不是用于检测正选择) 外，如果备选模型显著优于零假设模型(Likelihood ratio test, p < 0.05)，表明存在正选择作用 (Yang et al., 2000)。检测正选择作用的巢式模型中，M7 vs. M8 是检测正选择最为常用的模型对，但由于M8a vs. M8得到的假阳性结果更低，现在也渐渐被使用。当数据中检测到正选择作用后，可以应用贝叶斯经验贝叶斯法 (Bayes empirical Bayes, BEB) 进一步鉴定可能受正选择作用的氨基酸位点 (Yang et al., 2005)。

表1. 位点模型 (Site model) 中不同模型的背景假设

软件运行环境及信息

1. EasyCodeML程序为Java编写，支持跨平台使用。运行前需确保操作系统中已安装了JAVA运行环境 (JRE 1.6或更高版本)。JRE最新版本可从JAVA官网下载，下载链接为：https://www.java.com/en/download/manual.jsp
2. EasyCodeML v1.41 https://github.com/BioEasy/EasyCodeML
3. 本文示例数据可以访问作者永硕网盘获取（http://raindy.ys168.com/）；EasyCodeML软件子目录Example中的Example2.pml和Example2.tre文件与示例数据相同。
   

实验步骤

1. 分析前数据处理
   
   1.1

   比对序列文件
   用于选择压力的序列文件必须是比对后的PAML格式，扩展名为PML。非PAML格式的序列文件，可以通过EasyCodeML工具栏下的"Sequence Format Convertor"进行转换，通过拖曳操作将不同的序列格式转化为PAML格式 (图1)。格式转换前，需要检测编码蛋白的核苷酸序列长度是否为3的倍数。除外，序列中的终止密码子也必须删除，否则容易产生假阳性结果。
   
   
   图1. 比对序列格式转换
   
   
   1.2

   系统发育树文件
   用于适应性进化分析的树文件必须是不带任何数值的 Newick 格式。如果树文件中的数值未被清除，可以使用工具栏下的"Tree Cleaner"工具进行过滤 (图2)。
   
   
   图2. 树文件过滤器
   
   
2. 数据加载及分析参数配置
   

   EasyCodeML分析所需要的数据文件准备就绪后，可在Preset Mode (预置模式) 中"Setup"标签下加载对应的数据文件，并配置主要参数 (图3)。
   
   
   图3. EasyCodeML参数配置界面
   
   

   2.1

   选择一个本地文件夹作为工作目录，可以直接通过拖曳进行操作。
   2.2

   选择Site Model作为分析模型 (默认)。
   2.3

   载入已转为PAML格式的比对序列文件。
   2.4

   载入已过滤掉数值的 Newick 格式树文件。
   2.5

   "nt"为线程数，如果计算机资源支持多线程，可以将"nt"选项后的数值调至系统支持的最大值，但一般不超过 7 (Site Model中总共有7个模型)。
   2.6

   物种对应的密码子可以通过"icode"选项下拉菜单进行选择，示例数据为标准密码子，故使用默认密码子 (0: standard)。
   2.7

   通过"Check"检查树文件中类别名称与比对序列的名称是否一致，当弹出"SeqFile And TreeFile is Ready"时，说明两者名称一致，可以进行下一步的CodeML分析。
3. CodeML 分析
   启动CodeML分析前，需要点击"Save Current Profile"以使前面配置的参数信息生效。点击"Run CodeML"即可开始CodeML分析，其运行时间主要取决于数据中序列长度和数量。当CodeML运行完毕，程序会自动弹出"Running Finished"提示，在LRT标签下将对数似然法 (log‐likelihood value, lnL) 和参数个数值 (the number of parameter, np) 等关键参数显示在界面中(图4)。
   
   
   图4. 似然率检验结果
   
   
4. 结果分析及绘图
   LRT结果显示，M1a vs. M2a、M7 vs. M8和M8a vs. M8三对巢式模型的p值均小0.05，说明备选模型显著优于零假设模型，表明存在正选择作用。通过"Export"模块可将主要结果及参数导出，并自动生成为准发表级的Excel表格(图5)。
   
   
   图5. EasyCodeML 自动生成的准发表级数据表
   
       如图5所示，表中最后一列"Positive Sites"是由BEB法鉴定可能受正选择作用的氨基酸位点。虽然M2a中也检测了正选择作用，但其灵敏度未能鉴定受正选择压力的氨基酸点。相比之下，M8模型更为灵敏，示例数据中共找到了9个显著的正选择位点(后验概率 ≥ 0.95)。
       在位点模型中，模型对M0 vs. M3虽然不是用于检测正选择，但M3模型中生成的的rst文件，可用于dN/dS的滑窗分析(Sliding windows analysis) 以直观显示序列中不同位点的选择压力变化。如图6所示，可以提取倒数第二列数据进行绘图，横坐标为比对序列Codon位置，纵坐标为dN/dS比值，预览图如Gao et al. (2017) 一文的Figure 3所示，本文不再赘述。
   
   
   图6. 不同位点的dN/dS值
   

常见问题及解决方法

在分析过程中，EasyCodeML用户可以会遇到一些如下问题：

1. 为什么EasyCodeML 无法运行，解压缩包死循环？
   EasyCodeML.jar，是主程序文件，正常情况下，直接双击即可运行。但如果*.jar文件被识别为压缩包，则程序会被系统关联的解压缩软件进行解压缩处理。
   解决方法：是去除*.jar与解压缩软件的关联，比如WinRAR之类的，在参数设置中去除 jar 的关联，重新关联为JAVA程序即可。
2. 为什么找到的正选择位点在原始比对序列中找不到？
   在EasyCodeML中，默认启用了Clean data 功能(Clean data = 1)，即：比对序列中如带有gap，则gap所对应的一列序列被忽略不进行分析，从而导致位置发生偏移。
   解决方法：(1) 需要重新运行CodeML，分析前将 Clean data选项前的勾选取消，即Clean data = 0，此时数据含有gap运行CodeML分析，其运算量也相应增大；(2) 分析前将带gap的同一列碱基全部手工删除，找到正选择位点的氨基酸位置后，再还原对应到原始比对序列上。
   

小结与建议

EasyCodeML以CodeML为内核，可以实现CodeML分析界面化、自动化操作，生成准发表级的数据表格。应用EasyCodeML进行适应性进化分析，简单易用。但受内核所限，该软件仅适用于中、小型数据集。对于多基因的CodeML分析，推荐使用LMAP软件 (Maldonado et al., 2016)；对于基因组级别的数据，则推荐使用软件PosiGene (Sahm et al., 2017) 和GWideCodeML (Macías et al., 2020)。

致谢

EasyCodeML自2015年首次推出测试版，时至目前最新的1.4版，期间历尽无数次的反复调试。在此，特别感谢华南农业大学陈程杰博士的辛勤付出和广大EasyCodeML用户的反馈意见使得程序日臻完善。

竞争性利益声明

作者声明没有利益冲突

参考文献

1. Bielawski, J. P. and Yang, Z. (2004). A maximum likelihood method for detecting functional divergence at individual codon sites, with application to gene family evolution. J Mol Evol 59(1): 121-132.
2. Forsberg, R. and Christiansen, F. B. (2003). A codon-based model of host-specific selection in parasites, with an application to the influenza A virus. Mol Biol Evol 20(8): 1252-1259.
3. Gao, F., Zou, W., Xie, L. and Zhan, J. (2017). Adaptive evolution and demographic history contribute to the divergent population genetic structure of potato virus Y between China and Japan. Evol Appl 10(4): 379-390.
4. Macías, L. G., Barrio, E. and Toft, C. (2020). GWideCodeML: a python package for testing evolutionary hypotheses at the genome-wide level. G3: Genes Genom Genet 10(12): 4369-4372.
5. Maldonado, E., Almeida, D., Escalona, T., Khan, I., Vasconcelos, V. and Antunes, A. (2016). LMAP: lightweight multigene analyses in PAML. BMC bioinformatics 17: 354.
6. Sahm, A., Bens, M., Platzer, M. and Szafranski, K. (2017). PosiGene: automated and easy-to-use pipeline for genome-wide detection of positively selected genes. Nucleic Acids Res 45(11): e100.
7. Yang, Z. (2007). PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol 24(8): 1586-1591.
8. Yang, Z. and Nielsen, R. (1998). Synonymous and nonsynonymous rate variation in nuclear genes of mammals. J Mol Evol 46(4): 409-418.
9. Yang, Z. and Nielsen, R. (2002). Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol Biol Evol 19(6): 908-917.
10. Yang, Z., Nielsen, R., Goldman, N. and Pedersen, A. M. (2000). Codon-substitution models for heterogeneous selection pressure at amino acid sites. Genetics 155(1): 431-449.
11. Yang, Z., Wong, W. S. and Nielsen, R. (2005). Bayes empirical Bayes inference of amino acid sites under positive selection. Mol Biol Evol 22(4): 1107-1118.
12. Zhang, J., Nielsen, R. and Yang, Z. (2005). Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol Biol Evol 22(12): 2472-2479.