利用 OrthoFinder、IQtree、Notung、iTOL 绘制基因树

本文为以 4 个兰科物种中 NB-ARC 结构域的序列信息作为输入，构建 有根基因树。步骤参考文献 Genome- Wide Analysis of the Nucleotide Binding Site Leucine-Rich Repeat Genes of Four Orchids Revealed Extremely Low Numbers of Disease Resistance Genes，Front. Genet.，08 January 2020。

我们知道构建基因家族的有根基因树与基因家族的扩张和收缩分析 需要物种树作为支持（详情参见：OrthoFinder 2.0 原理及所涉及的相关概念 ）。本文通过 OrthoFinder 推断物种树 。

OrthoFinder 虽然可以生成基因树，但是组成基因树的基因必须来自于 OrthoFinder 通过 MCL 算法聚类出的 Orthogroup。如果待研究的基因家族基因无法被聚在 1 个 Orthogroup 中，则无法通过 Orthogroup 将所有基因家族成员绘制在 1 个基因树中。为了提高分析的灵活性，本文使用 IQtree 推断无根基因树 ，然后利用 Notung 调和无根基因树与有根物种树，得到有根基因树并完成基因家族的扩张和收缩分析。最后利用 iTOL 对基因树进行可视化呈现。

一、使用数据

1. 构建物种树使用的是 4 个兰科物种的全部蛋白质序列信息。
2. 构建基因树使用的是 4 个兰科物种蛋白组中 NB-ARC 结构域序列信息（结构域长度 ∈ [ 75% mean-length，125% mean-length ]）。
   
   
   

二、分析流程

1. 使用 OrthoFinder 推断有根物种树

• 有关 OrthoFinder 的介绍详见 OrthoFinder 2.0 原理及所涉及的相关概念

使用 orthofinder 根据 4 个兰科物种的蛋白序列，构建物种树。

$ orthofinder -f protein_fasta/ -t 50 -a 50
-t number_of_threads：多序列比对及建树等所调用的 CPU 数
-a number_of_orthofinder_threads ：多序列比对及建树外其他操作所调用的 CPU 数 

OrthoFinder 将文件夹 protein_fasta 内所有以 .fa、.faa、.fasta、.fas、.pep 结尾的文件作为输入，读取文件内的蛋白序列。有根物种树结果如下图（图 1）所示（通过 megax 查看）

与参考文献中物种树结构一致（图 2），说明了 OrthoFinder 结果的准确性 。由于 OrthoFinder 无需外群即可生根，为了判断无外群生根的准确性，本实验以 Amborellla trichopoda 作为外群进行建树，命令如下：

$ orthofinder -b protein_fasta/OrthoFinder/Results_Jan11/WorkingDirectory/ -f Atrichopoda/

• -b <dir1> -f <dir2> 中 <dir1> 是指原先分析输出中的 WorkingDirectory 文件夹路径，<dir2> 是指新添物种的蛋白序列文件夹路径。命令表示在原先分析 <dir1> 的基础上添加 <dir2> 中物种的蛋白序列进行分析。这种做法相比从头分析省去了原先蛋白序列间比对的时间，OrthoFinder 将只进行新添序列间及新添序列与原序列间的比对。物种树结果如下：

可以发现，不加外群和加外群物种树的拓扑结构是一致的，验证了 OrthoFinder 不加外群构建有根物种树的准确性 。

• Amborella trichopoda，无油樟，双子叶植物，其被认为是被子植物中最早分化出来的一支，和其他几乎所有的被子植物都没有明显的亲缘关系，常被用作基因分析的外群。Amborella 的基因组数据下载自 JGI phytozome
  https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Atrichopoda
  https://data.jgi.doe.gov/refine-download/phytozome?organism=Atrichopoda&expanded=291。

由于 OrthoFinder 不支持输入用户定义的 Orthogroup，而 4 个物种的 NB-ARC domain 会被 OrthoFinder 分割为多个 Orthogroup 并分别构建基因树，无法将 4 个物种的 NB-ARC 构建在 1 棵基因树中。这里使用 IQtree 来推断基因树结构。

2. 使用 IQtree 推断无根基因树

• 有关 IQtree、MUSCLE 的介绍详见 MUSCLE、IQtree 软件及使用简介

使用 MUSCLE 对 NB-ARC domain 蛋白序列进行多序列比对。

muscle -in NB-ARC.domain.fasta -out NB-ARC.domain.afa -maxiters 2
-in 输入文件 
-out 输出文件名，输出文件默认为 Fasta 格式
-maxiters 最大迭代次数

将多序列比对的结果作为 IQtree 的输入，推断无根基因树结构。

$ iqtree -s NB-ARC.domain.afa -nt 50 -m JTT+F+R10 -n 2 -alrt 1000
-s  输入文件
-nt 线程数
-m  指定使用模型：JTT+F+R10
-n  指定进化树迭代次数：2
-alrt 使用SH-aLRT检验并指定重复抽样次数：1000

252 个 NB-ARC domian 蛋白序列信息中 6 个蛋白没有通过 IQtree 的 chi2 检验，最终构建出包含 246 个蛋白的无根基因树。

3. 使用 Notung 根据有根物种树为无根基因树生根且推断基因复制、转移、丢失事件

Notung 输入：有根物种树、无根基因树

• 物种树必须包含基因树中所有基因的物种，额外的物种将被 Notung 忽略。
• 输入的基因树中的每个基因需要包含物种信息（格式：gene_species）。Notung 以 “_” 作为分隔符，基因名中最后一个 “_” 后内容作为物种名。如 Notung 认为 XP_020599319.1_Phalaenopsis_equestris 的基因名为 XP_020599319.1_Phalaenopsis，物种名为 equestris。所以如果物种名中包含 “_” 需替换为 “-” 或其他字符。

在 Notung 中导入基因树、物种树后，通过 Rooting Mode 将无根基因树转化为有根基因树。Rooting Mode 根据有根物种树为每个边计算 DTL 分数（DTL 分数越小的越适合做根），Notung 会高亮（红色）最小值及附近 [(max-min) × 5% ] 的边。用户通过鼠标点击确定选择那个边作为根。将生根后的基因树以 NEWICK 格式导出。

• Notung 的 Reconciliation Mode 功能的输入是 有根基因树，通过比较基因树和物种树推断基因复制、转移、丢失事件。
• Notung 的 Rooting Mode 功能的输入是 无根基因树，通过比较基因树和物种树推断最可能的生根边，根据用户选择的根推断基因复制、转移、丢失事件。本文使用的是 Rooting Mode 功能。

将有根物种树、无根基因树导入后发现 Notung 会计算出许多适合的生根位点（下图红线）。在没有外群的情况下，难以确定适合的生根位点。

为此，本文原先的数据集中加入外群（Amborella trichopoda）后重新进行 MUSCLE、IQtree 分析。Notung 计算出适合的生根位点如下图所示。可以发现，凭借外群能方便寻找适合的生根位点。从 Notung 的底边栏可以看到 Notung 推断 NB-ARC 发生了 254 次基因复制事件，76 次基因丢失事件。将生根后的基因树以 NEWICK 格式导出。

4. 利用 iTOL 在线工具绘制基因树图

iTOL：https://itol.embl.de/upload.cgi

将 Notung 生成的有根基因树文件作为输入，提交至 iTOL。图中

• 土黄：Amborella trichopoda（外群）
• 红色：Dendrobium catenatum
• 紫色：Phalaenopsis equestris
• 青色：Apostasia shenzhenica
• 绿色：Gastrodia_elata