主要

系统发育树描绘了生物实体之间的进化关系。这些实体可以是物种,就像生命之树一样1,2,3,4.它们也可以是肿瘤进展树中的癌细胞5或发育谱系树6、传染病暴发中的病毒和细菌菌株7、细胞或树木中的基因,用于在模型和非模型物种之间传播分子功能注释8,9.由于这种普遍性,推断系统发育树的方法在所有生命科学中是最常用和引用最多的软件工具之一。

在物种树推断的背景下,全基因组测序的可用性使得考虑每个分类单元的标记基因与基因组提供的标记基因一样多成为例行公事。这种“系统发育”方法解决了真核生物生命树的许多关键方面,例如被子植物深部分支之间的关系10,海鞘在脊索动物中的位置11,蜕皮动物分支12,Lophotrochozoa分支13以及主要多足类分支之间的关系14,等等。

然而,尽管测序的质量和成本迅速提高15,16,推断系统发育树所需的数据分析仍然非常费力且计算密集17.系统发育研究需要多个昂贵的步骤,每个步骤都可能是主要的研究工作(图1)。1):原始读取的管理,从头组装通常包括多轮纠错和使用一种或多种技术搭建脚手架18、重要基因的注释和表征、直系同源基因的鉴定和比较,以及从直系同源标记进行树推断。目前的最佳实践通过长读长和短读长测序等技术组合以及跨多个管道的多轮参数优化来优化这一过程。

图1:策略和管道说明。
图1

a,Read2Tree 旨在在使用许多物种时绕过许多耗时且昂贵的管道步骤来获得系统发育树,因此从读取到树。b,Read2Tree 管道概述。

目前的趋势是对更多的物种和样本进行测序。地球生物基因组计划于2018年11月启动,旨在在未来十年内对“地球上所有150万种已知的动物、植物、原生动物和真菌物种”进行测序19.组成联盟在简化和优化测序和注释过程方面取得了进展,但正字推理和树推理步骤仍然极具挑战性。与此同时,各个实验室正在进行大量的基因组测序活动,每项研究的样本量为数百到数千个基因组变得很普遍16.然而,根据感兴趣的物种,往往缺乏高质量的参考基因组,而且单个实验室往往缺乏计算基础设施或专业知识,无法在各个分析步骤中充分利用数据。这在主要联盟主导的研究中得到了体现,这些研究需要数年时间和数百万美元来阐明某些感兴趣物种的进化,或者最近使用各种管道来评估变异并报告严重急性呼吸系统综合症冠状病毒 2 (SARS-CoV-2) 的组装。因此,一个主要的瓶颈是对这些大规模数据集的统一分析,以避免某些偏差或伪影。

在本文中,我们介绍了 Read2Tree,这是一种推断物种树的方法,它的工作原理是将原始测序读数直接处理成相应的基因组,绕过基因组组装、注释或全序列比较。与当前已建立的管道相比,Read2Tree能够在很短的时间内提供数百个样本的完整系统发育比较。至关重要的是,在不影响生成的树的准确性的情况下实现了加速。此外,Read2Tree还能够仅使用低覆盖率(0.1×)数据集以及RNA与基因组测序来提供准确的树木和物种比较,并可进行长读长或短读长。这使得 Read2Tree 成为一种高度通用的方法,可以从单个样本中获取关键见解,可扩展到数千个样本。为了建立这种方法,我们评估了其在跨越不同王国的基因组和转录组数据集、发散时间和测序技术上的性能。随后,我们应用 Read2Tree 构建一棵大型酵母生命树,并将其用于比较 SARS-CoV-2 样本——从而突出了 Read2Tree 的准确性(例如,与国家生物技术信息中心 (NCBI) 分类相比)和速度。

结果

最先进的系统发育管道需要许多步骤,这些步骤既耗时又容易出错(图 1)。使用Read2Tree,我们直接处理原始测序读数,并重建传统树推理方法的序列比对(图1)。图1b和补充图1b。1). 我们首先将原始读数与来自全基因组参考直系同源组(OG)的核苷酸序列进行比对;我们使用Mafft20默认值)(图 1)。1b,1)。在每个OG中,我们从与参考序列对齐的reads中重建蛋白质序列(图1)。1b,2)。重要的是,参考OG中的这些序列并不局限于单拷贝标记基因,例如线粒体细胞色素c氧化酶I基因或BUSCO基因21;它们还包括多个旁系同源基因以及非通用基因。这是通过利用从我们实验室开发的直系同源基质 (OMA) 资源中分析的 2,500 个不同基因组计算的 OG 来实现的22,23.接下来,我们保留最佳参考引导的重建序列,使用重建的核苷酸碱基的数量作为标准(图1)。图1b、3和补充图1b。随后,将选定的共识添加到OG的多序列比对(MSA)中(图1)。1b,4)。最后,可以使用常规方法进行假定的 OG 选择和树推理(我们使用 IQTREE24默认情况下;无花果。1b,5)。有关各个步骤的更多详细信息,请参阅方法

这样一来,Read2Tree就能够绕过基因组组装、注释、同源性和正字推断,在很短的时间内报告推定OG的关键信息。此外,由于每个样本都是独立处理的,因此Read2Tree可以并行处理输入基因组,并根据输入基因组的数量进行线性缩放。

覆盖范围和参考距离对精度的影响

我们在各种条件下测试了Read2Tree,包括两种序列(DNA与RNA)、三种靶种(拟南芥酿酒酵母和Mus musculus)、三种测序技术(Illumina、PacBio和Oxford Nanopore Technologies(ONT))、六级测序覆盖率(范围从0.2×到20×)和六组不同的参考物种(与跨越10亿年进化的目标越来越远) (图。2a). 对于序列重建精度(图 2a).2b),我们测量了重建序列的正确性(“精确度”)和重建序列的完整性(“召回率”)。对于树木重建的精度(图1)。图2c和补充图2c。6),我们将重建的树与已知的物种系统发育进行比较,并报告了重建树的精确度和召回率,以至少90%支持的分支而言。

图2:使用三种不同数据集、六种不同覆盖水平和三种测序技术的Read2Tree基准。
图2

a, 参考数据集的系统发育树。深紫色(底部)是用于制图的物种。颜色表示物种移除,以评估参考数据集中最近邻的依赖性。时间点是从 timetree.org 获得的59.b,Read2Tree序列更相似(百分比同一性)和更完整,覆盖率增加,与更密切相关物种的距离减小。Illumina数据的最佳序列鉴定性。颜色传达了与最近的参考物种(参考文献)的进化距离的增加。c,在将树枝折叠到90%以下后,使用Read2Tree重建的树木的精确度和召回率。

总的来说,Read2Tree在序列重建方面能够保持高精度(图1)。2b)和树木重建(图2b)。2c) 在所有数据集中,根据数据集难度的不同,召回级别也不同。首先,我们评估了单个数据集的覆盖范围为0.2×至20×的影响。我们观察到,增加测序覆盖率对精密度影响不大,主要是降低召回率:在大多数配置中,即使覆盖率低至0.2×,Read2Tree也可以在序列水平上保持90-95%的精确度(图1)。小鼠的转录组短读长数据上获得了最好的低覆盖率结果,在0.2×覆盖率下,精密度达到98.5%。为了评估Read2Tree的多功能性,我们在DNA和RNA数据集中对其进行了基准测试。这通常没有太大的影响,但转录组RNA结果(在小鼠数据集中)受平均覆盖率差异的影响略小,这可能是由于这些数据中不均匀的基因表达水平导致的覆盖率差异很大(图1)。2b,c)。接下来,我们评估了Read2Tree是否能够利用当前的测序技术。为此,我们将其应用于传统的短读长、Oxford Nanopore 和 PacBio 长读长读长。为了实现这一点,Read2Tree 内置了略有不同的映射策略,用于长读取和短读取(方法)。如图。图2b,c显示,Read2Tree在每种测序技术中都保持了高准确度,但我们观察到与传统的短读长相比,其准确度最高。我们尚未评估可能改变这一结果的最新测序技术,如PacBio HiFi或Illumina infinity。

最后,我们评估了Read2Tree相对于手头样本与参考集中最近亲戚之间的进化距离的鲁棒性。这通常很关键,因为人们可能不知道组装的最接近的祖先,或者它不可用25.因此,我们在从700万年前到11亿年前的广泛进化距离范围内测试了Read2Tree。虽然这些肯定是极端情况,但总体而言,Read2Tree 能够成功应对它们。图2b,c显示参考集的选择主要影响回忆,更近的参考基因组导致更多的重建位置。值得注意的是,即使在参考文献非常遥远的数据集中,Read2Tree也能够保持高精度,例如,在参考文献集中没有任何脊椎动物基因组的情况下处理小鼠RNA测序(RNA-seq)数据。

我们还在模拟数据上测试了 Read2Tree,覆盖率在 0.1× 到 10× 之间,到最近参考的距离在 2 到 150 个点接受突变 (PAM) 单位之间变化——其中 100 PAM 对应于平均每个位点一个替换。除了最极端的情况(PAM >120或覆盖率<0.5×;补充图。7).

鉴于物种、覆盖率、测序技术、测定(DNA 和 RNA)和模拟数据的广泛基准,我们观察到 Read2Tree 确实是一种高度通用且准确的工具,可以直接从原始读数中重建系统发育。

比基于装配的树更快,而且通常更准确

接下来,我们将 Read2Tree 的性能与传统的装配流水线进行了比较。为此,我们使用 Canu 在与上一节相同的数据集中生成了从头组装和蛋白质预测26用于 PacBio 和 ONT 数据以及 Megahit27与 SoapDeNovo 一起28Illumina读数(方法)。常规组装使用OMA独立处理,包括相同的导出参考基因组,因为OMA独立之前被证明可以识别最准确的系统发育标记基因29.为了在用于树推断的串联对齐中包含直系同源标记,我们要求通常设置的最小阈值为 80% 的分类单元存在。如上所述,我们通过删除参考树上的物种来改变数据集中最接近的剩余物种(图 1)。2a). 使用不同的覆盖率和参考集,我们获得了每个物种的 42 个数据点。对于这些数据点中的每一个,我们分别执行了正字推理并记录了其计算时间。放置在各个OG中的序列比例显示出高水平的变异(补充图1)。对于蛋白质组的每个组装和变异,我们计算了组装或Read2Tree所得树与使用拟南芥酿酒酵母的高质量基因组组装获得的树之间的拓扑距离。

3 显示了总体结果,突出显示了 Read2Tree 的性能。也许不出所料,我们观察到覆盖率水平对基于组装的方法的性能产生了深远的影响,使它们无法处理低于 5-10× 的覆盖率。因此,对于这些数据集,我们只报告 Read2Tree 结果。

图 3:Read2Tree 与具有汇编、正字预测和 MSA 计算的常规流水线的比较。
图3

a,使用参考树与 Read2Tree 的树或来自组装方法的树之间的差异来比较树。对于深蓝色,我们只有 Read2Tree 树,因为无法获得这些低覆盖率的程序集。低于零(深蓝色或浅蓝色)时,Read2Tree 更准确,而高于零(红色)时,组装方法更准确,灰色表示方法之间没有差异。b,从读取到级联 MSA 可用性所需的墙时间的比较,显示了可用最近剩余参考和覆盖率的依赖性。

在可以比较这两种方法的情况下,传统的从头组装方法优于Read2Tree的唯一情况是覆盖率高,并且与最接近的参考物种非常遥远(>500 Mya)(图1)。3a,每个图形的右上角区域)。在所有其他情况下,Read2Tree 的准确性优于传统方法。具体来说,在更高覆盖率的酵母数据集上,组装和Read2Tree总体上表现良好——我们从未观察到在获得的树和参考树之间超过两个不同的分支。在至少 10× 的覆盖率和遥远的参考物种的情况下,传统的组装方法优于 Read2Tree(图 1)。图3a和补充图3a。4).

相比之下,在更复杂的拟南芥肌肉分枝杆菌数据集上,Read2Tree优于组装方法,与参考的差异较小(Read2Tree最多两个不同的分支,而传统方法最多四个)。在ONT数据上,Read2Tree的特点是读取时间更长,但错误率更高,在两个数据集上都优于传统方法。

最后,就计算时间而言,Read2Tree通常比传统方法快得多,在较大的基因组上快100倍(图100)。图3b和补充图3b。8b)。

总而言之,这些结果表明,Read2Tree在所有条件下都更快,并且在低覆盖率数据集和其他传统方法完全失败的数据集(长读长转录组学)中生成可靠的树。在更高的覆盖率水平下,Read2Tree推断的树木在质量上可以与那些从具有完整管道的组装参考物种中获得的树木相媲美,特别是当应用于更复杂的基因组时,除非最近的参考物种非常遥远(>5亿年)。

我们还将 Read2Tree 与 Mash 进行了比较,这是一种基于 k-mer 的快速方法30通常用于细菌基因组。虽然 Mash 的无对齐方法甚至比 Read2Tree 快得多,但生成的树远不如 Read2Tree 或基于汇编的方法准确(补充图 1)。说明了为什么像 Mash 这样的无比对方法虽然对快速近似非常有用,但通常不适合重建高质量的系统发育树。

精确重建 435 种酵母生命树

为了评估Read2Tree的潜在大规模应用,我们将其应用于从原始读段中重建大型酵母系统发育。得益于Read2Tree处理低覆盖率数据集的能力,我们可以将分析扩展到NCBI序列读取档案(SRA)数据库(2018年11月,404个物种)中可用于出芽酵母的所有Illumina单端和双端、ONT、PacBio和454个测序读取数据集,以及从OMA数据库(2018年发布,3,063个OG)获得的31个参考物种。使用自动检索和绘图方法,我们能够获得 404 个物种的直接序列(补充文件 1)。Read2Tree可以在大约一个月的计算中处理这些数据集(按顺序添加每个物种,并在30个中央处理器(CPU)上并行执行映射 - 每个参考一个CPU),由于其“令人尴尬的并行”架构,每个样本都独立处理,直到系统发育推断(10× Illumina:~20分钟,使用四个线程)。

这些数据集中有很大一部分最近被用于构建 363 种萌芽酵母物种的系统发育31.这包括一个包含 196 个新组件及其注释的数据集31.这项巨大的努力将酵母生命树划分为13个主要分支,并强调了水平基因转移对酵母物种进化的影响31.由于最先进管道的复杂性,它还消耗了数百万个 CPU 小时和多年的工作。此外,传统的基于组装的方法无法将低覆盖率样品纳入其分析中。我们能够使用Read2Tree使用一小部分资源来扩展这项工作。

使用 Read2Tree,我们能够在 435 个样本(包括 31 个物种作为参考)中计算和生成这种大型系统发育。一些样本由于覆盖率太低(约3.1)而失败×假设平均基因组大小为12 Mbp。尽管如此,使用Read2Tree,即使在覆盖率低于5×的水平下,我们也能够包括多个样本,据报道,OG中放置了2,500多个序列(补充图1)。14). Read2Tree能够重建系统发育,并报告了每个样品组装的系统发育相关基因,总体上显示出与参考数据相似的GC水平(补充图1)。15). 我们没有观察到每个物种放入OG的序列数量与其个体覆盖率之间的相关性,这也证明了这一点。14,相关系数0.2)。

考虑到共同物种的子集,我们的结果与Shen等人的结果高度一致。31(图。4 和补充图。1213):两棵树都表现出与 NCBI 分类树相似的距离——我们的树有 297 次分裂,而 Shen 等人有 291 次分裂。在直接比较中,Shen et al. 和 Read2Tree 彼此之间更相似,只有 128 个不同的分裂(分支的差异为 20%),而不是 NCBI 分类法。在支撑位低于 90 的树枝倒塌后,保守的 NCBI 树和我们的分裂次数差异为 29 次,我们和 Shen 等人之间的分裂次数差异为 25 次。其中 24 个分裂是 Read2Tree 和 Shen 等人共同的。为了更深入地了解这些差异的性质,我们评估了与NCBI分类法的一致性,以达到两个不同分辨率的水平:科和属。在较粗的科水平上,Read2Tree与6个科的NCBI分类法更一致,而Shen等人在一个科中更一致(补充图1)。10). 在更精细的属水平上,Read2Tree对4个属的分类学更符合NCBI的分类法,而Shen等人则为10个属(补充图10)。11).

图4:Read2Tree与最先进的系统发育管道之间的高度一致性31.
图4

最上面一行显示完整的树和用于计算外圆树的对齐矩阵。红点表示引导程序低于 100 的节点。物种Naumovozyma dairenensis,以前被错误分类32,33,以红色突出显示。最底行显示修剪成重叠叶子集的树木。

尽管如此,Read2Tree 和 NCBI 分类法之间仍然存在某些差异。虽然解决大多数此类情况本身就构成了整个后续研究,但我们能够解释一个明显的分歧:Naumovozyma dairenensis 被归入 CUG-Ser1 分类,而根据 NCBI 分类学,它应该是囊菌酵母在 Saccharomyceces sensu lato 组中 酵母 Saccharomyces sensu lato 组在 Saccharomycetaceae 家族中。但是,这是文献中报告的错误元数据的情况32,33.

鉴于这种系统发育,我们现在可以在几分钟内使用 Read2Tree 轻松更新和扩展它,并生成其他序列。这样可以深入研究酵母的比较基因组学,并进一步探索菌株之间的差异及其对生命、食品生产等的影响。对于其他生物体来说,这也很容易重现,因为Read2Tree能够跨越相对于参考树的较大进化距离。

用于人畜共患病监测和人类流行病学的Read2Tree

为了进一步说明Read2Tree的多功能性,我们用它来重建一个系统发育,包括来自OMA冠状病毒数据库的各种冠状病毒,以及存放在SRA的215个原始冠状病毒测序样本。 除了推定的SARS-CoV-2序列外,我们还纳入了两个来自蝙蝠(SRR11085797(参考文献)的样本。34)和SRR11085736(参考文献)。35)) 和一个来自水貂36(SRX9605666)。

重建的系统发育与从UniProt参考蛋白质组获得的谱系分类完全一致。特别是,该树不仅恢复了主要的冠状病毒属(Alpha-、Beta-、Gamma-和Delta冠状病毒),而且还恢复了所有完全一致的亚属(图1)。5).

图5:Read2Tree正确地对最近的SARS-CoV-2序列进行了分类,并概括了各个变体的演变。
图5

所有属(整个树中的灰色框)和亚属(彩色标签)都已正确划定。插图的重点是树中包含 215 个 SARS-CoV-2 样本的部分,并且关注的变体(彩色标签)一致地聚集在树上,表明 Read2Tree 可用于对样本进行分类。

第一个蝙蝠样本对应于 RaTG13 的读数,RaTG13 是迄今为止发现的 SARS-CoV-2 的近亲34.事实上,在我们的树中,它正好落在 SARS-CoV-2 分支之外。另一个蝙蝠样本也可以被确认为甲型冠状病毒即鼻病毒亚属35.同样,我们可以确认水貂样本的分类,作者将其鉴定为甲冠状病毒亚属 Minacovirus36.

SARS-CoV-2序列在冠状病毒生命之树中的位置也与我们先前对它们的了解一致。参考基因组,2020年1月初报道的Wuhan-胡-1序列(参考文献)。37),位于子树的底部。在它之前分支出来的唯一三个序列是SRR11092056-8,它们是在大流行开始时从严重肺炎患者那里获得的34.最后,我们注意到,分析中包含的关注变体清楚地显示为树上的不同分支。

为了实证测试我们方法的可扩展性,我们还使用 Read2Tree 处理了 10,283 个 SARS-CoV-2 样本。重建的树根据疾病控制中心关注的变体分类对序列进行聚类,进一步证明该工具可用于快速可靠地对 SARS-CoV-2 变体进行分类(补充图 1)。17). 对于其他对照组,同样的观察结果也成立——仅使用编码基因标记运行 Read2Tree(补充图 17)。16),并使用 FastTree38作为树推理方法(补充图。18).

总体而言,Read2Tree在不同冠状病毒序列中的应用说明了该工具能够同时处理该病毒家族的大量系统发育呼吸39以及对单个 SARS-CoV-2 关注变体进行分类所需的深度。这使得Read2Tree适用于人畜共患监测和人类流行病学40.

讨论

我们介绍了Read2Tree,这是一种扩展和简化比较基因组学费力过程的方法:组装,注释和系统发育比较。这些步骤的计算成本高昂且容易出错,并且需要专业知识。使用Read2Tree,我们可以直接从原始读数中重建系统发育相关基因,从而能够以最低的计算和覆盖要求对手头的物种进行放置和比较。与最先进的管道相比,该方法的效率使得使用一致的方法并行处理大量样品成为可能,并且不会影响准确性。

目前大规模比较基因组学或一般比较基因组学项目的固有问题最近从获得准确的组装转向对这些组装的注释和管理。这在一定程度上是由于测序技术的进步导致了长读长16,18,也是由于汇编算法的创新41,42.这些步骤仍然需要很高的DNA质量,而且通常更昂贵,但可以进行大型项目,如脊椎动物基因组计划43,人类泛基因组44和端粒到端粒45项目。然而,在所有这些情况下,基因组的注释以及连续性和准确性方面的改进仍然是主要瓶颈。此外,我们还表明,Read2Tree能够在很短的时间内对所有三种测序技术(Illumina、ONT和PacBio)进行准确分析。此外,大型联盟也可以从运行 Read2Tree 中受益,尽管拥有高覆盖率的数据集,以独立质量控制 (QC) 他们的组装和树构建方法。

一个主要优点是,尽管从头组装是侧步,但Read2Tree可以在没有密切参考基因组的情况下运行;事实上,我们展示了准确的树木重建,包括从数亿年差异分开的物种中测序读数。尽管我们也达到了这种鲁棒性的一些限制,但当Read2Tree同时受到非常高的发散性和低测序覆盖率的影响时,应该注意的是,随着越来越多的物种在生命之树上进行测序,进化距离将趋于减少。

此外,虽然基因组资源的大多数作者将注释集与组装的序列放在一起,但并非所有作者都这样做。直接从原始读段处理基因组的能力不仅规避了这一限制,而且还可以减少因过度依赖特定参考基因组而产生的偏差。已经有一些初步的努力来“非人化”非人类类人猿基因组46,但许多其他分支仍然受到类似偏差的影响,通过处理原始读取可以大大减少这种偏差。

我们在大规模酵母数据集上展示了 Read2Tree 的速度和准确性。在这里,Read2Tree 能够从直接从 SRA 检索的原始读取样本中重建高质量的树。 尽管覆盖率水平存在差异和其他可能的技术偏差,但还是实现了这一目标。

在第二个说明性应用中,我们根据原始冠状病毒测序数据重建了一棵树,其中包括来自正在进行的 SARS-CoV-2 大流行的 10,000 个样本。在这里,Read2Tree再次能够正确地分类和放置所有样本,无论是在冠状病毒科的整个范围内,还是在SARS-CoV-2样本之间的微小变异深度上,其中系统发育标记基因的最佳选择通常取决于序列差异的水平47.

我们还将 Read2Tree 与超快、无对齐方法 (Mash) 进行了比较,其中 Read2Tree 实现了更高的精度(补充图 1)。5). 在目前的形式下,Read2Tree具有与Kraken2等宏基因组分类器不同的功能。48)或离心机49.事实上,虽然这些工具试图利用已知的特征序列进行读取级分类,但Read2Tree旨在通过推断用于系统发育树推理工具的大型多位点输入数据矩阵来有效地提取全基因组(或转录组范围)的系统发育信号,这一步骤已被证明对于解决困难的系统发育至关重要17,29,50,51,52.尽管如此,Read2Tree可以进一步开发,以处理宏基因组样本 - 通过将其与基因组分箱预处理步骤相结合。近年来,已经提出了许多不同的基因组分箱方法,无论是通过“差异覆盖”方法,该方法利用样本之间的相关丰度来识别来自同一物种的读数53,54,55,使用Hi-C协议,这使得在物理上接近地识别DNA的某些部分成为可能56,57或单细胞技术58.

总体而言,Read2Tree是一种重建系统发育重要基因和表征手头样本或整个样本集合的方法,无需预处理,计算资源很少,生物信息学专业知识也很少,因此可以研究大量基因及其进化。这将使系统发育重建工作更快、更全面——从微小的病毒基因组到大型真核生物基因组,还有细胞谱系、癌症树和生物学和医学中的其他类型的系统发育。

方法

Read2Tree 方法的说明

Read2Tree 将各种公开可用的工具用于其某些步骤(MAFFT20、NextGenMap60和 Samtools61),并以结构化的方式使用它们,从读取和引用 OG 到直接输入树推理工具(默认为 IQTREE)的串联对齐24.为此,它需要两组输入数据:(1)一组可以直接从OMA数据库获得的参考OG,以及(2)来自单个物种的要映射的reads。Read2Tree 管道按以下方式工作。首先,它使用REST-API从OMA浏览器中从选定的参考OG中检索DNA序列,然后将这些序列分类到每个物种的一个文件中。同时,它使用 AA 序列和 MAFFT 计算比对20,然后使用密码子信息生成 DNA 比对。计算完成后,将所有读段映射到DNA参考物种,并构建共识序列(局部组装)。由于我们的本地组件是参考引导的,因此它们永远不会比 OG 的最长或最短序列部分更长或更短。然后,使用最佳选定参照的坐标将局部组件放置到对齐中。因此,不需要新的比对,我们可以确保正确的AA / DNA被放置在比对中的正确位置。然后将每个 OG 的结果对齐连接起来,并计算一棵树。有关 Read2Tree 内部工作的更多详细信息,请参阅补充图。1.

Read2Tree 可以在映射步骤中使用多个实例并行化。建议先计算参考集。然后,可以拆分映射步骤,以便每个映射都可以在高性能集群上作为单个作业提交执行。

OG选择

OG选自OMA62使用标记基因导出功能 (https://omabrowser.org/oma/export_markers/)。对于所有物种,覆盖物种的最大数量设置为0.8,最大标记数量设置为-1(无限制)。所选物种显示在图中。1一个

拟南芥酿酒酵母的全基因组测序读数来自SRA数据库,用于PacBio、Illumina和Oxford Nanopore技术。还从SRA数据库中获得了所有三种技术的信使RNA-seq读长。读取的子采样是在 Python 中执行的(参考文献)。63).对于 PacBio 和 ONT 读数,对子采样进行了优化,使累积碱基数符合预期覆盖率。对于覆盖率测试,假设小鼠的累积基因长度(转录组)为 38 Mbp,而 thale cress 和 12 Mbp 酵母基因组长度为 12 Mbp,则对读段进行子采样。对读数进行采样以获得 20×、10×、5×、1×、0.5× 和 0.2× 的覆盖率。大酵母树的读数是从SRA数据库(补充文件1)中获得的。冠状病毒的读数是从SRA数据库(补充文件1)中获得的。所有 SRA 编号均可在补充文件 1 中找到。

参考树构造

使用图中定义的物种计算了三个评估物种的参考树。2a. 物种选自OMA62如 OG 选择中所述。使用 MAFFT 对单个 OG(基因标记)进行比对20版本 7.310 (–maxiter 1000–local),并使用 IQTREE 推断树24版本 1.6.9 (-m LG -nt 4 -mem 4 G -seed 12345 -bb 1000)。对于用于测试对参考数据集的依赖性的参考树,从现有对齐中删除特定物种,并使用 IQTREE 计算树,如前所述。所有参考树均可在补充文件 2 中找到。为了突出进化的年份,我们使用时间树收集了时间59(2022 年 4 月)。

Read2Tree 运行

对于单一物种运行(图。23),Read2Tree 使用默认参数运行。对于大酵母树(图1)。4),Read2Tree分多个步骤运行。首先,获得参考数据集(使用 –reference 选项)。然后,并行化映射,以便对于每个物种单独执行针对单个参考的映射(使用 –single_mapping 选项)。这意味着,对于每个物种,执行了 31 次并行映射。此外,假设基因组长度为 12 Mbp,则将读取覆盖率超过 20× 的物种采样至 20× 覆盖率。最后,将所有映射的物种合并在一起并连接起来,以提供多序列输出(使用-merge_all_mappings)。

准确性评估

我们通过获取放置在 OG 中的每个 Read2Tree 重建序列(针对每个物种、覆盖率、技术和去除水平)来评估序列重建的准确性,并对其原始 OG 进行 blastp(ncbi-blast,2.8.1 版)搜索,其中包含来自目标物种的高质量组装的原始序列。准确度以原始细胞百分比身份来衡量,召回率为所有OGs的串联MSA中获得的氨基酸总数。此外,我们评估了 Read2Tree 重建序列的顶部命中是否与其组装的同物种对应序列最相似,或者用作重建参考的序列或该特定 OG 的任何其他随机序列部分(补充图 1)。3).

组件

对于这三个物种,全基因组数据是按照最佳实践或默认参数,使用单独的测序技术特定组装程序组装的。对于Illumina,我们首先使用了megahit27(版本 1.2.9)具有用于组装重叠群的默认参数。随后,SOAPdenovo28(版本 2.04-r241)用于脚手架:首先,SOAPdenovo-fusion -D -K 41 -c megahit.contigs.fa -g scaffold_prefix -p 20,然后是 SOAPdenovo-63mer map 和 scaff,并在配置文件上提供推荐参数。对于 ONT 读取,我们使用 Canu 组装读取26(版本 2.0)具有指定的基因组大小 (genomeSize) gnuplotTested = true -nanopore-raw 和 useGrid = false 参数,以仅在集群上的一个节点上本地运行它。最后,对于 PacBio 连续长读长数据,我们还使用了具有类似参数的 Canu(2.0 版),但指定了 -pacbio-raw 参数。所有运行时间均使用 linux 时间测量,并记录墙和 CPU 时间。RNA-seq数据的组装方式与全基因组不同。对于Illumina RNA-seq,我们使用了Trinity64(版本 2.8.5),并具有以下参数:–seqType fq–max_memory 50 G–左 reads1.fq.gz–右 reads2.fq.gz–CPU 6–trimmomatic–full_cleanup–输出前缀。它们会自动执行 Trimmomatic 并遵循 Trinity 的建议。

组装基因组的正字学预测

对于每个程序集(物种、技术和覆盖级别),我们使用 Slurm 调度器在 UNIL HPC 集群上运行 OMA 独立(版本 2.3.3)。为此,我们收集了图中描述的所有物种。2 使用 OMA All versus All 导出功能。然后,我们根据图去除相关物种。2、每次在集合中添加小鼠、酵母或水芹的组装,并使用标准参数运行正字预测(OMA 版本 2.2.1)。因此,例如,对于Illumina M. musculus 10×组装,我们对所有参考数据集运行了七次OMA,并且与其最近亲戚的距离越来越近。我们总共运行了 126 个不同的 OMA 运行,其中包含 7 种参考蛋白质组变体和 3 种技术变体,对拟南芥S. cervisiea 具有三个覆盖级别。此外,我们还对M. musculus 5×、10×和20×Illumina组装体进行了21次OMA。all-versus-all 部分在 1,000 个节点上并行化,最后一部分在具有 40 G 内存的单个节点上运行。为了获得用于树推断的 OG,我们应用了 0.8 的分类占用阈值,如前所述。根据 Shen 等人的程序过滤 OG(见下文)。OG使用MAFFT单独对齐20版本 7.310 (–maxiter 1000–local) 并连接,并使用 IQTREE 推断树24版本 1.6.9 (-m LG -nt 4 -mem 4G -seed 12345 -bb 1000)。

树与树的比较

使用多个树距离度量将每个 Read2Tree 树与拟合参考进行比较。为了实现拓扑相似性,我们使用了两种方法,一种是使用 Robinson-Foulds 距离并计算两棵树之间不同拆分的数量,另一种是将每个节点折叠起来,并支持引导程序,然后计算重叠拆分的数量。然后,我们将重叠拆分数除以引用中的拆分数,将精度定义为重叠拆分数除以 Read2Tree 树中的拆分数。

大酵母树

对于大型酵母树,我们于 2018 年 11 月从 SRA 中提取了所有可用的酵母数据集(406 个物种,补充文件 1),并将 Read2Tree(标准参数)应用于从 OMA 数据库(2018 年 11 月)中提取的 31 个酵母物种使用标记导出功能,最小物种覆盖率为 0.8(3,082 个 OG)。所选物种可在补充文件3中找到。使用 Read2Tree 根据其排序方法映射 SRA 数据库中的读取。为了将我们的分析与Shen等人进行比较,我们的目标是拥有尽可能多的共同物种。为此,我们用测序读数来补充我们的树,这些测序读数是我们从15个物种的组装基因组中模拟出来的,这些物种存在于Shen等人的树中,但在我们的数据集中缺失(补充文件1)。使用InSilicoSeq(版本1.3.0 https://github.com/HadrienG/InSilicoSeq,–model hiseq -n 600000)进行模拟。从Shen等人的树中绘制物种图。31在我们的树上,我们利用 ete3 的 NCBI 界面获得了物种/菌株的分类单元鉴定(参考文献)。65).对于无法自动定位的物种,我们使用NCBI分类界面(https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi)获得了分类单元鉴定。

过滤OGs酵母,如Shen等人。

鉴于重建的序列放置在各自的OG中并添加到它们的比对中,我们决定按照ref的协议计算一棵树。31.简而言之,从 3,082 个比对中,我们选择了包含超过 171 个物种的比对,从而产生了 1,829 个 OG。然后,我们使用 phyutils 2.2.6 (seqs -aa -clean 0.01) 来清理对齐。由于我们的方法没有将来自同一物种的多个序列放入一个OG中,因此我们跳过了推定的旁系同源物的删除。在对齐方式中,我们将所有“X”更改为间隙字符“-”。然后,我们应用了 trimAl 版本 1.4.rev15 (-gappyout)。接下来,我们删除了长度短于它们所属的每个 OG 的修剪 MSA 长度的 50% 的蛋白质序列。我们还删除了修剪后的MSA总长度为<167个氨基酸位点的OG。这导致了 926 次对齐。通过这些对齐,我们使用具有自动模型选择功能的 IQTREE(版本 1.6.9)来计算树。然后,我们在基因树中鉴定出分支长度超过所有分支长度中位数20倍的物种。我们从各自的排列中删除了这些物种,再次控制了超过 171 个物种。然后,我们使用 IQTREE (-seed 12345, -m LG + G4, -bb 1000, -nt 20) 计算树。

大酵母树比较

使用所有分类单元鉴定,我们检索了当前的NCBI参考分类法和每个物种的分类。然后,我们比较了三棵树(NCBI,Read2Tree和Shen等人。31) 使用重叠叶集上的 Robinson-Foulds 距离。此外,我们在树上叠加了 Shen 等人的分类。最后,我们使用从NCBI分类信息中提取的特定分组(目、科和门)的祖先节点比较了包含最多单系物种的祖先节点。所有比较均使用自定义 Python Jupyter 笔记本进行。此外,我们还收集了有关GC内容和输入覆盖率与映射率的数据。用 ete3 (ref.65).七巧板图是使用 dendextend R 库生成的66.使用 phylo.io 获得并排拓扑比较67.

冠状病毒科(Coronaviridae)树重建

标记基因从至少四个物种的 https://corona.omabrowser.org/ 输出。这些基因的DNA序列是从同一资源中获得的。使用自定义脚本添加了四个具有 SARS-CoV-2 参考基因组基因间区域的额外组。我们从参考基因组MN908947组装中提取了至少 30 bp 的连续块,这些块既没有被任何 CDS 区域覆盖,也没有被不属于 https://corona.omabrowser.org 资源中任何 OMA 组的蛋白质(即 ORF8 和 ORF10)覆盖。这导致了四个区域(1..265;26473..26522;27760..27893;29675..29903),我们将其视为额外的组。SARS-CoV-2 样本是从 Nextstrain open (https://data.nextstrain.org/files/ncov/open/global/metadata.tsv.xz)7.使用自定义 Python 脚本(包含在下面的链接存储库中)获得了具有跨越所有不同分支的 SRA 种质的不同样本。SRA读取种质以及Nextstrain 的分支注释可在补充文件 1 中找到。读取是从 SRA 数据库下载并修剪的。将Read2Tree应用于该数据集,所有获得的reads都映射到标记基因。Read2Tree 使用标准参数运行。通过去除间隙超过 70% 的色谱柱来过滤所得的超基质比对。这删除了 30,969 列,产生了大小为 295 × 42,669 的超矩阵。最后,使用 IQTree2 (ref.24) (版本 2.2.0-beta) 参数 -m GTR -ninit 2 -me 0.05。作为其他控件,我们使用 FastTree 计算了树38版本 2.1.11 而不是 IQTREE2,并且没有额外的四个额外组。所有树木均可在补充文件 1 中找到。

对于包含 10,283 个样本的放大实验,我们使用了相同的协议,除了读取注释的来源。在这里,我们使用了 https://harvestvariants.info/ 的分支注释(种质和注释可在补充文件 1 中找到)。

模拟系统发育分析

模拟的系统发育包括使用 ALF 包的 15 个物种树的固定拓扑68(版本 0.99)。我们将导致其中一个物种(感兴趣的物种)的分支长度更改为 2 PAM 和 150 PAM 之间。对于每次运行,我们都会推断 OMA 组(不包括感兴趣的物种)。然后,使用 art_illumina69,我们生成了长度为 100 和 150 bp、覆盖率为 0.1 至 10 的 DNA 测序读段(配对端)。接下来,对于每种情况,我们运行 Read2Tree 来推断系统发育。最后,我们根据ALF的输出计算了推断物种树和真实物种树之间的Robinson-Foulds度量。

与 Mash 的比较

我们以已建立的程序集作为我们从 NCBI 下载的参考。随后,我们使用了 Mash(2.3 版)草图30大小为 10 m(默认为 k = 21),然后用 Mash 距离获得基因组之间的距离,并根据该参考集分析读数。最后,我们应用了 RapidNJ70(版本 2.3.2)在从 Mash 获得的距离矩阵上推断物种树。我们在提供的参考文献中针对不同的距离执行此操作,始终将来自拟南芥的读数与组件进行比较。

报告摘要

有关研究设计的更多信息,请参阅本文链接的《自然投资组合报告摘要》。