基于扩增子数据的系统发育树的构建和展示   

仪器设备

普通个人电脑 (Windows10系统64位版、CPU ≥ 双核、内存 ≥ 4 G、硬盘 ≥ 20 GB)

软件和数据库

1. gitforwidnows 2.23.0 (http://gitforwindows.org)
2. R 4.0.3 (https://www.r-project.org)
3. Rstudio 1.2.5019 (https://www.rstudio.com/products/rstudio/download)
4. USEARCH v10.0.240 (https://www.drive5.com/usearch/download.html)
5. MUSCLE (http://www.drive5.com/muscle/)
6. IQ-TREE v2.0.3 (http://www.iqtree.org)
7. trimAL (http://trimal.cgenomics.org/downloads)
   

软件的安装和使用

一、 首先在C盘根目录新建名为bin的目录，其具体位置为C:\bin。

二、 USEARCH软件的下载和安装
USEARCH软件 (http://www.drive5.com/usearch/download.html) 是Robert C. Edgar开发的一款超快的扩增子数据分析软件，该软件在序列比对、原始数据质控、OTU聚类、多样性分析等多领域广泛应用 (Edgar, 2013)。在使用前需选择合适的USEARCH软件版本进行下载，其下载界面为 (如图1)。


图1. USEARCH10软件的下载

注：32位usearch为免费版，但限制使用内存4GB，64位为收费版本，没有内存使用限制。选择接受许可协议，版本必须选择v10.0，选择填写邮箱，提交收到链接，下载后改名为usearch10并将软件放在Windows10系统中C:\bin目录中。

三、 R语言的下载安装和使用
R语言是目前生物学、经济学等领域最流行的统计分析语言，下载最新版R语言 (下载页面：https://cran.r-project.org)；点击Download R for Windows完成R-4.0.3.win.exe安装程序的下载(如图2)；双击安装程序，建议语言选择英文安装。


图2. R语言的下载

四、 IQ-TREE的下载安装和使用
IQ-TREE软件是2015年发表的一款快速准确进行最大似然法 (Maximum Likelihood, ML) 构建系统发育树的软件 (Nguyen et al., 2015)，目前已经更新到v2.2.2版本。IQ-TREE软件应用一种快速、有效的随机算法，在近似的计算时间内具有比RAxML软件和PhyML软件更高的精确度。此外IQ-TREE软件的模型选择速度比jModelTest快10-100倍，其自展支持率估算比RAxML软件快10-40倍并且支持宏基因组等大数据计算 (Minh et al., 2020)。IQ-TREE软件的下载界面： (http://www.iqtree.org)，选择下载最新版64位的IQ-TREE (v2.2.2) 软件 (如图3)，解压后将其放在C:\bin目录中。


 图3. IQ-TREE软件的下载

五、 MUSCLE的下载安装和使用
MUSCLE软件是一款快速多重序列比对软件，MUSCLE软件具有比CLUSTALW等软件更快的比对速度以及精确度，它能在数分钟内完成数百条序列的比对。迄今，MUSCLE软件已经被引用了超过37,000次，是生物学领域中最为广泛使用的软件之一 (Edgar, 2004)。MUSCLE软件的下载界面 (如图4)：(http://drive5.com/muscle)，选择下载最新版的MUSCLE软件，解压后将其放在C:\bin目录中。


 图4. MUSCLE软件的下载

六、 trimAL 的下载安装和使用
完成精确的多序列比对后，需过滤掉一些低质量以及高变异度的序列区域。trimAL软件能快速，精确切除和过滤低质量以及高变异度的序列，仅保留进化保守的区域用于后续分析。trimAL软件的下载界面 (如图5)：(http://trimal.cgenomics.org/downloads)，选择下载最新版的trimAL软件，解压后将其放在C:\bin目录中。


 图5. trimAL软件的下载

七、 Rstudio的下载安装和使用
Rstudio的下载页面：(https://rstudio.com/products/rstudio/download)。从网页中选择下载最新版的Rstudio（如RStudio Desktop 1.3.1093），双击安装程序进行默认安装。Rstudio安装完成后，按如图所示步骤调出“Terminal”界面，然后在Terminal窗口中输入ls (LS的小写)，按回车进行测试，如果出现“command not found”错误，请按照下图重新进行操作和设置 (如图6)。


 图6. Rstudio软件的设置及Terminal的调用

八、 “git for windows”的下载安装和使用
“git for windows”软件 (v2.28.0) 是一款能在Windows系统下运行的命令行工具，能在Windows下运行部分Linux代码 (下载页面：https://gitforwindows.org/)，按照默认参数右键管理员安装Git-2.28.0-64-bit.exe即可 (如图7)。


图7. “git for windows”软件的下载

九、 添加程序位置至Windows系统中的环境变量
我的电脑-右键属性-按右侧截图操作，测试是否安装成功：在RStudio的Terminal下输入：usearch10，按回车，如有出现USEARCH10的版本信息，则表明安装成功。若不成功，可检查环境变量配置，按下图进行操作 (如图8)：


图8. Windows10系统的环境变量的设置

分析步骤

一、 数据获得和操作流程概况
准备输入数据
本文中所有的测试数据集、所需软件、R脚本以及生成的结果均已上传到百度网盘，如需测试和使用，可点击如下链接进行下载和使用 (百度网盘链接： https://pan.baidu.com/s/15MIJzV0_5kNV3VF_o0rKZw 提取码：5ru9)。具体技术流程概况以及各个软件的功能与彼此之间的联系如图9所示。


图9. 系统发育树的构建流程图

二、 数据处理相关的命令行操作均在Rstudio “Terminal”界面进行

1. 首先在C盘新建名为Tree_data目录
   #切换到Tree_data目录中
   $ cd C:\Tree_data
   $ mkdir -p result/tree
   $ cd result/tree
   
   
2. 研究者可以根据实际情况0.001 到 0.01的阈值筛选高丰度OTU。
   #统计OTU表中OTU数量，代码如下：
   $ tail -n+2 ../otutab.txt | wc -l
   
   #按相对丰度0.2 %筛选高丰度OTU，代码如下：
   $ usearch10 -otutab_trim ../otutab.txt -min_otu_freq 0.002 -output otutab1.txt
   
   #统计筛选OTUs表特征数量，代码如下：
   $ tail -n+2 otutab1.txt | wc -l
   
   #提取ID用于提取序列，代码如下：
   $ cut -f 1 otutab1.txt | sed '1 s/#OTU ID/OTUID/' > otutab_high.id
   
   
3. 在进行完OTU筛选后要根据OTUs的ID提取每个OTUs对应的fasta格式的代表性序列，手动整理物种注释和分组信息表annotation.txt，如图10所示。
   #筛选高丰度菌/指定差异菌对应OTUs的代表性序列，代码如下：
   $ usearch10 -fastx_getseqs ../otus.fa -labels otutab_high.id -fastaout high_otus.fa
   $ head -n 10 annotation.txt
   
   
   图10. annotation.txt文件中包含的内容
   
   

三、 序列对齐及系统发育树的构建
#构建进化树，实现高丰度菌的进化树的分组信息展示与美化。
#起始文件为 result/tree目录中high_otus.fa (序列)、annotation.txt (物种和相对丰度)文件
# Muscle软件进行序列比对和对齐，代码如下：
$ cd Tree_data/result/tree
$ muscle -in high_otus.fa -out otus_aligned.fa

#trimAL软件进行低质量以及高变异度的序列的过滤和修剪，代码如下：
$ trimal -in otus_aligned.fa -out otus_aligned_trimed.fa -gt 0.95

#利用IQ-TREE软件进行ML系统发育树的构建，代码如下：
$ mkdir -p iqtree
$ iqtree -s otus_aligned_trimed.fa -bb 1000 -redo -alrt 1000 -m MFP -nt AUTO -pre iqtree/training_otus

#参数简介：
-m 参数：指定模型选项，MFP表示ModelFinder Plus (自动默认)。
-redo参数：之前运行成功后生成了相应的文件，指定redo会重新跑一遍覆盖之前的文件。
-pre参数：将结果输入到iqtree文件夹中，且生成文件的前缀为training_otus。
-alrt参数：是否启用 SH-aLRT检验。

四、 iTOL网站进行系统发育树的编辑和展示
在运用iTOL在线工具进行系统发育树的美化之前，首先要使用“table2itol.R”这个R包 (https://github.com/mgoeker/table2itol) 生成用于系统发育树编辑和美化的注释文件。然后访问并登陆iTOL网站 (http://itol.embl.de/)，上传otus.nwk，再拖拽以下命令行生成的不同注释文件 (分别位于“plan1”、“plan2”和“plan3”三个文件夹) 于iTOL主界面的系统发育树图上即完成系统发育树的美化。
# plan1生成外圈颜色、形状分类和丰度文件，代码如下：
$ cd Tree_data/result/tree
$ Rscript ../../script/table2itol.R -a -c double -D plan1 -i OTUID -l Genus -t %s -w 0.5 annotation.txt

# plan2生成丰度柱形图注释文件，代码如下：
$ Rscript ../../script/table2itol.R -a -d -c none -D plan2 -b Phylum -i OTUID -l Genus -t %s -w 0.5 annotation.txt

# plan3生成热图注释文件，代码如下：
$ Rscript ../../script/table2itol.R -c keep -D plan3 -i OTUID -t %s otutab.txt

#参数简介
-a：找不到输入列将终止运行 (默认不执行)
-c：将整数列转换为factor或具有小数点的数字
-t：偏离提示标签时转换ID列
-w：颜色带，区域宽度等
-D：输出目录
-i：OTUs列名
-l：OTUs显示名称如种/属/科名
注：当需要标注的颜色过多时，R脚本会采用形状+颜色的方式对类别进行区分。

结果与分析

1. IQ-TREE运行完成后会在iqtree文件下生成多个文件，主要包括程序运行日志training_otus.log、ML树文件 (含有UFBoot或BP/SH-aLRT评估分支置信度) 和系统发育树文件training_otus.contree，本文测试数据中生成的系统发育树文件名称为training_otus.contree (如图11)。
   
   
   图11. IQ-TREE软件建树生成的结果文件
   
   
2. 首先进入iTOL在线网站，点击右上角注册 (仅限新用户) 和登陆，(如图12)：
   
   
   图12. iTOL在线网站主界面
   
   
3. 完成账号登录后，点击My Tree按钮，进入个人主界面，然后点击“Upload tree files”上传文件夹中的training_otus.contree文件，获得如下树形图 (如图13)。导入树文件之后，可以在iTOL在线网站的右上角选择下图红框中的“Basic”和“Advanced”进行系统发育树的编辑 (比如树形的变换、自展支持率的显示、字体大小和颜色的调整、分支的位置变换等)。iTOL在线网站具有强大的系统发育树编辑和美化功能，研究者可根据需求进行各种个性化的调整。关于编辑的具体使用方法可以参考iTOL的官方帮助文档 (https://itol.embl.de/help.cgi) 和iTOL官方视频资料 (https://itol.embl.de/video_tutorial.cgi)。
   
   
   图13. iTOL在线网站的系统发育树展示和编辑
   
   
4. 系统发育树的编辑：按住鼠标左键，将“plan1”文件夹中的“iTOL_labels-Genus.txt”文件拖到iTOL网页的当前主界面上，iTOL在线网站会自动将所有OTUs替换成其物种注释对应的属名 (如图14)。
   
   
   图14. 系统发育树中分支节点名称替换
   
   
5. 按微生物门分类水平进行系统发育树着色：按住鼠标左键，将“plan2”文件夹中的“iTOL_treecolors-Phylum.txt”文件拖到iTOL网页的当前主界面上，iTOL在线网站会自动按门水平对系统发育树进行着色 (如图15)。同时注意，此处着色可分别进行标签、分支和全树着色。
   
   
   图15. 系统发育树图按门水平进行着色
   
   
6. 添加其它系统发育树分组注释：使用三个不同分组的高丰度OTUs做进化树，可以把三个分组的高丰度OTUs的相对丰度，用柱状图形式进行展示。按住鼠标左键，将“plan2”文件夹中的“iTOL_simplebar-A.txt”、“iTOL_simplebar-B.txt”和“iTOL_simplebar-C.txt”文件拖到iTOL网页的当前主界面上，iTOL在线网站会自动将“KO”、“OE”和“WT”三个分组的高丰度OTUs用柱形图形式进行展现 (如图16)。
   
   
   图16. 添加高丰度OTUs不同分组柱形图的系统发育树图
   
   
7. 在“plan1”、“plan2”和“plan3”文件夹中有很多其它系统发育树美化和编辑的文件。例如，“plan3”文件夹中热图的添加等，研究者可以根据需求进行灵活添加，最终达到系统发育树的构建和展示目的。最后，研究者可以点击iTOL在线网站右上角的“Export”按钮，导出编辑完成的系统发育树图。本文最终生成的系统发育树图 如图17所示。
   
   
   图17.最终生成的系统发育树图
   

小结

本文简要介绍了微生物扩增子数据中高丰度OTUs数据的筛选，代表性序列及对应物种注释的获取，以及系统发育树的构建方法等。该流程展示了一套完整的扩增子数据构建系统发育树操作步骤，以帮助研究者学习和掌握接近发表质量要求的系统发育树图的构建方法。研究者在建树过程中可以根据研究领域参考同行文献对系统发育树进行一些细节参数的调整，并使用AI (Adobe Illustrator) 对特殊图形、字符和树进行进一步修改和美化。

致谢

感谢“git for windows”软件开发者提供的Git 软件 (https://github.com/git-for-windows)及mgoeker在GitHub网站上开发和公开分享的R语言包“table2itol” (https://github.com/mgoeker/table2itol)。本文分析方法已应用于待发表文章“Distribution and variations of mycotoxin producing fungal community in major rice production areas of China”。

参考文献

1. Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res 32(5):1792–1797. 
2. Edgar, R. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods 10(10): 996–998.
3. Minh, B. Q., Schmidt, H. A., Chernomor, O., Schrempf, D., Woodhams, M. D., von Haeseler, A. and Lanfear, R. (2020). IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol 37(5): 1530–1534.
4. Nguyen, L. T., Schmidt, H. A., von Haeseler, A. and Minh, B. Q. (2015). IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol 32(1): 268–274.