基于基因树-物种树调和的物种树定根方法

引言

系统发育树描述了基因、物种或其他生物实体之间的进化关系，在组织和解释生物多样性方面发挥着核心作用。我们解释系统发育树的大多数方式都需要对进化树进行定根，这样才能知道进化的时间方向。例如，推断祖先-后代关系或确定顶端分类群之间的关联程度都需要有根的树。比较系统发育分析，如重建祖先状态，也需要一颗有根的树。这是系统发育学中许多困难的来源，因为所有最常用的替代模型都是时间可逆的。时间可逆模型是指状态（氨基酸或核苷酸）1和2在两个方向上的预期替换数相同的模型，假设两种状态的频率处于平衡状态（即，其中为状态的平衡频率，为从初始状态经历时间后变为状态的概率）。具有此属性的模型在计算上易于处理，但其似然值不受根位置的影响，因此推断出来的树本质上是无根的。请注意，即使是在树查看软件中打开时看起来有根的树，在大多数情况下都是任意根的，应该在无根模式下更准确地查看。

考虑到标准方法的局限性，最广泛使用的树定根方法是在分析中包含外群。外群是一组已知的分类群，根据先验信息，这些分类群在剩余的内群分类群的最后一个共同祖先之前已经分化。在推断出一棵无根树后，可以将根放在将外群与内群分开的分支上。这种方法在某些情况下效果很好，但当外群与内群只有较远的亲源关系时，可能会导致问题。这在深层微生物进化分析中经常出现的一种情况——在细菌和古菌的根的推断中考虑到外群的长分支——但也跨越生命树。在这些情况下，外群定根会在树中引入一条很难相对于内群放置的长分枝，并会诱发系统发育错误，例如长枝吸引，其中群内的长枝分类群被吸引而靠近外群。较远外群的另一个问题是，它们的分子序列可能很难与内群序列自信地对齐，或者它们可能共享较少的共同基因，从而限制了可用于推断内群分类群关系的数据量。在其他情况下，可能没有或不存在合适的外群。例如，不能使用外群来确定整个生命树的根。在这些情况下，必须找到另一种定根的策略。

也许最广泛使用的外群定根的替代方法是基于分子钟概念的方法，即分子进化速率（序列替换）随时间恒定的假设。这些方法使用沿树的序列分歧来推断根的位置。中点定根——在最长的两末端的中间位置定根——是一个简单的方法，它保证如果根据严格的分子时钟序列进化，返回正确的根。中点定根容易受到违反分子钟假设情况的影响，例如沿着单个末端分支的加速进化可能会导致推断错误的根。最近，一些如最小祖先偏差 (minimal ancestor deviation) 的方法，旨在通过搜索后代分支之间的最小速率变化的根，利用所有末端到末端距离的信息，可以提供概念上于中点定根相似但更稳健的结果。松弛分子钟 (relaxed molecular colocks) 可以明确地模拟谱系之间的进化速率变化，并结合化石校准的信息，也可以用于搜索最有可能产生序列数据的根，即最大化其边缘似然值。然而，与在给定固定根的情况下估计内部节点年龄的后验分布相比，计算这种边缘似然值需要大量计算，因此在实践中很少使用RMC方法来推断根。与RMC一样，非平稳和不可逆替代模型也可用于推断根，因为数据的似然值取决于根的位置，尽管在实践中它们尚未得到广泛应用。

另一方面，调和方法 (reconciliation methods) 的目的是从基因树和物种树之间的系统发育不一致中提取定根的信息，结合来自基因复制和转移的信号。基因树-物种树调和工作通过使用一系列基因和物种层的进化事件（如基因复制、转移和丢失）将基因树绘制到物种树中。基于简约思想的调和方法，在给定预先指定的事件成本的情况下，寻求意味着事件最少的调和，而概率方法依赖于沿着物种树的基因树进化的生成模型。概率方法有两个重要的优点：i）它允许估计模型的参数，尤其是复制、转移和丢失率；ii）使用序列进化和基因树-物种树调和的联合模型，以一致的方式考虑基因树的不确定性。后一点很重要，因为仅基于序列的基因树重建几乎总是涉及在统计上等价或弱区分的关系之间进行选择，这些关系可以更好地基于假定的物种树来解决。的确，基于联合似然值的基因树重建，也就是所谓的物种树感知的方法，比单独使用序列重建的要精确得多。因此，预期会导致更少的、会干扰定根的虚假复制、转移和丢失事件。

本文的重点介绍了一种基于物种树的调和，即混合似然估计 (ALE)，用于定根基因树和物种树的方法 (包含一批脚本程序)，该方法还可对有根物种系统发育中的祖先基因内容进行推断。ALE方法在微生物基因组分析方面有一些优势。调和模型包括基因转移，因此ALE可用于分析经历过一个或多个HGT的、不能用于串联或超树分析的基因家族。这意味着ALE可用于推断有根物种树，其使用的原核基因家族比例远高于可通过串联或超级树方法进行分析的小子集，因为基因家族树可能因水平转移以及家族特有的复制和丢失事件而彼此不同。从数据中可以推断出复制 (duplication)、转移 (transfer) 和丢失 (loss) 事件的发生率，该方法整合了单基因树中的过度不确定性，以及它们与物种树调和的方式。模拟结果表明，ALE可以准确地为物种树定根，并且可以从数据中准确估计复制、转移和丢失事件。

使用ALE估算物种树根的工作流程

数据的准备

ALE软件要求FASTA格式的基因组 (Prokka注释的faa文件) 中，序列的ID要用下划线_来分割基因组ID和序列编号（例如：‘GCA_005280195.1_1'）。其中，基因组ID也是基因组注释faa文件的文件名。

推断基因家族

使用OrthoFinder2进行基因家族推断。

基因家族序列的比对

大多数比对软件要求每个基因家族至少有四条序列。因此，在比对前需要需要筛选可比对的基因家族。

即使是小数据集，也可以产生数千、甚至上万个基因家族，导致线性计算需要相当长的时间；为了提高计算速度，可以GNU parallel等并行软件。进入包含所有可比对的基因家族序列文件的目录内，执行：

ls *.fasta | parallel mafft --auto {} ">" {}.aln

该命令以自动模式在当前目录中以fasta结尾的每个文件上运行应用程序mafft。它将一次运行多个作业，最多可运行计算机上的核心数。--auto自动模式将选择最佳比对模型进行比对。然后可以使用BMGE等软件修剪比对结果。下面的shell脚本并行执行BMGE，从比对中修剪对齐不好的位点。

ls *.aln | parallel java -jar BMGE.jar -i {} -t AA -m BLOSUM30 -of trimmed_{}

推断基因家族的自展树分布

然后对这些比对结果，用IQ-TREE来估计中每个基因家族的ML树，并使用超快自展树的分布来捕获基因树的不确定性。基因树的分布在ALE中用于估计条件分支概率，进而估计基因树概率。因此，也可以用基因树的MCMC样本（例如，使用PhyloBayes进行采样）作为ALE的输入，这其实也是ALE最自然的输入数据。在实践中，对于大型数据集来说，这可能在计算上非常耗时，所以超快自展树集是一种可接受的折衷方案。

输入的基因树质量对于后续分析非常重要，因此替代模型的选择就成为了关键。下面的命令通过添加两个混合模型的模型搜索来执行IQ-TREE，并创建超快自展树分布。

iqtree2 -s <Trimmed_gene_family_ID.aln> -m MFP -madd LGC20,LGC60 --score-diff ALL -B 10000 -wbtl

要将自展树分布转换为ALE对象，使用ALEobserve：

ls *.ufboot | parallel ALEobserve {}

推断无根物种树

推断无根树可用supermatrix和supertree两种方案。supermatrix方案将所有单拷贝直系同源基因分别比对后进行串联，视其为一个supermatrix然后用IQ-TREE等软件构建物种树。supertree方案，将所有单拷贝直系同源基因分别构建基因树后，用ASTRAL来整合所有基因树进而获得物种树。

为物种树确定候选根位置

为了评估每个根位置的可能性，无根物种树可以手动地在不同候选位置上定根。 手动定根可以用FigTree等可视化软件处理，也可以用脚本 root_tree_in_position.py处理。

python3 root_tree_in_position.py unrooted_st.newick species1 species2

执行该脚本可以在unrooted_st.newick树上，species1和species2两物种最近的祖先分枝上定根。

调和候选有根树

使用ALEml_undated，将每个ALE对象与候选有根物种树进行调和。ALEml_undated至少需要两个参数：候选有根物种树和ALE文件。有时可能需要调整调和模型的其他几个参数。例如，通过指出每个基因组估计缺失的部分，可以在分析中考虑基因组的不同完整性。fraction_missing缺失比例值可通过BUSCO或CheckM来估计 (假设基因组完整性的估计值为80%，那么fraction_missing就等于0.2)。要添加fraction_missing参数，只需添加：fraction_missing=fraction_missing.txt。fraction_missing.txt包含每个序列的相应缺失比例 (每行记录一个基因组的缺失比例：“Sulfolobus:0.20”)。

ls *.ale | parallel ALEml_undated root_1_st.newick {}

ALEml_undated使用的调和模型的参数可以用各种命令行参数进行修改（请参见ALE GitHub更多细节）。例如，可以使用delta、tau和lambda参数将D、T或L速率单独设置为零 (tau=0导致无任何基因转移的调和，delta=0导致无任何基因复制的调和)；否则，这些参数的值将通过最大似然法进行估计。

根据物种树候选根位置的数量，调和有根树的过程需要执行相应的次数。为了区分不同的运行过程与结果，ALEml_undated的执行需要在不同的目录内完成。

比较对候选根位置的支持

基因家族的似然值可以用来比较物种树不同根的支持度，就像传统系统发育学中的位点似然值可以用来比较不同基因树的支持度一样。CONSEL软件可以完成相关的统计检验，在给定的位点（或基因家族）似然值集合中选择树。CONSEL软件需要的似然值数据，记录在保存于为每个候选根位置作调和分析的目录内的.uml_rec中。为提取这些信息，可执行脚本程序：

python3 write_consel_file_p3.py root1 root2 root3 root4 > likelihoods_table.mt

命令行中输入的如root1是1号候选根位置的调和结果目录。似然值表likelihoods_table.mt中，记录的似然值顺序和命令行中输入的目录顺序一致。执行AU检验，需要完成以下步骤：

makermt likelihoods_table.mt consel likelihoods_table catpv likelihoods_table

所得结果根据相应检验的P值来确定哪些候选根的位置是可接受的 (P>0.05)，哪些候选根位置是被拒绝的 (P<0.05)。

至此通过ALE的基因树-物种树调和方法来为物种树定根的分析就完成了！

评估根信号的性质（鲁棒性检查）

在传统的系统发育学中，位点剥离 (site stripping)——去除快速进化或组分有偏差的位点——通常用于评估系统发育信号的稳健性。其基本原理是，这些位点中的信息可能会产生误导（或至少难以充分建模），例如，与真实的进化信号相比，更倾向于长分支 (或组分) 吸引。在整个树中，从快速进化或有组分偏差的位点获得支持的关系，被认为不如对移除这些位点具有鲁棒性的关系可靠。在系统基因组学水平的根或拓扑推断中也存在类似的情况：基因家族在大小、进化速率和DTL事件发生率方面存在差异，考虑到可用的方法，一些家族可能比其他家族更容易建模。与传统案例类似，从特定数据子集（例如，丢失率或转移率非常高的基因家族）获得支持的关系可能被认为不如在整个家族分布中获得支持的关系可靠。与传统的位点过滤一样，这一过程可能是探索数据的最佳方式。目前尚不清楚筛选“异常”位点或基因家族是否总能改善推断，也没有明确的阈值，例如，快速替代位点或频繁转移的基因在该阈值下可被确定为不可靠。

脚本DTL_ratio_analysis_ML_diff.py可以用来评估高复制、转移和丢失率（DTL）基因家族的影响。该脚本根据基因家族的DTL率从分析中依次删除基因家族，并重新计算每个根的总似然值，以及与最大似然根的似然值差异。因此，从这些分析中，我们可以可视化特定种类的基因家族（例如，DTL率特别高或特别低的基因家族）在整个分析中对给定根的支持程度。

DTL_ratio_analysis_ML_diff.py根据复制率、转移率或丢失率，或丢失/物种形成率、复制/物种形成率或转移/物种形成率，创建3个图。一次只能测试上述指标中的一个。第一个图显示了基于所选度量的家族按顺序移除后不同根的总似然值，第二个图显示了每个根与最大似然值根的似然值差异，第三个图显示的与第二个图相同，但仅适用于具有高物种代表性的基因家族。该脚本需要两个名为roots_to_test.txt和species_list.txt的文本文件。roots_to_test.txt文件应包含uml_rec文件的目录名。物种列表species_list.txt文件应包含数据集中的所有物种，并应与调和分析中使用的物种树的格式相匹配。该脚本还需要两个命令行参数：由AU检验确定的最大似然根目录的名称，以及要按D、T、L、LS、DS或TS删除基因家族的度量。该脚本还可仅使用高物种代表性调和文件重新绘图。在与write_consel_file_p3.py运行的相同目录中运行DTL_ratio_analysis_ML_diff.py脚本。

最可能有根物种树的基因含量进化

一旦确定了最可能的根位置，这种技术允许用户量化基因内容进化中复制、转移、丢失和起源的相对贡献。这些预测值可以针对物种树的每个分支进行计算，从而能够深入了解基因含量如何随时间演化。branchwise_number_of_events.py脚本，用于计算物种树所有分支的复制、转移、丢失和起源事件的预测总数，并估计每个内部树节点的祖先基因组中存在的基因数量。

python3 branchwise_number_of_events.py > dtloc.tsv

ALE的输出还提供了每个节点上存在的基因家族的估计值。因此，我们可以在内部节点模拟基因家族的存在和缺失，重建祖先基因组。 Ancestral_reconstruction_copy_number.py脚本可以用来重建内部节点上的所有基因，有三个命令行参数：第一个是基因家族在调和分析中出现在该分支上的百分比，即0.5意味着出现在50%的调和中 (如果想进行更严格的祖先重建，增加百分比)；第二个和第三个参数是内部节点的范围。

python3 Ancestral_reconstruction_copy_number.py 0.5 6 10

然后，可以分别使用eggNOG mapper和BlastKOALA，使用基因本体和KEGG数据库对感兴趣的序列进行功能注释。