摘要:序列比对是许多生物信息学分析的基础。许多比对工具通过将序列分割为连续的固定长度片段(称为k-mers)开始。使用较长且独特的seed比对速度更快,而使用较短且避免突变的seed则更准确。在此,我们介绍了X-Mapper,旨在通过包含缺口的动态长度seed(称为
研究论文
● 期刊:Genome Biology(IF:10.1)
● DOI:
https://doi.org/10.1186/s13059-024-03473-7●原文链接: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-024-03473-7
● 第一作者:Jeffry Gaston
● 通讯作者:Eric Alm (ejalm@MIT.EDU), An-Ni Zhang (章安妮) (anni.zhang@ntu.edu.sg)
● 发表日期:2025-1-22
● 主要单位:谷歌、麻省理工学院、南洋理工大学
摘要Abstract
序列比对是许多生物信息学分析的基础。许多比对工具通过将序列分割为连续的固定长度片段(称为k-mers)开始。使用较长且独特的seed比对速度更快,而使用较短且避免突变的seed则更准确。在此,我们介绍了X-Mapper,旨在通过包含缺口的动态长度seed(称为gapped x-mers)实现高速和高精度比对。分析人类参考序列时,我们观察到次优比对减少了11–24倍,而在细菌参考序列中不一致性降低了3–579倍,相较于其他比对工具,在非目标菌株和物种上的比对分别改进了53%和30%。其他基于seed的分析算法也可能受益于gapped x-mers。
引言Background
短读长测序应用广泛,例如微生物组的宏基因组学研究,数据分析要求研究人员处理海量数据集并应对多样化的基因组。在处理短序列的所有步骤中,序列比对是计算最密集的任务之一,比对的质量直接影响下游分析的效率和准确性。因此,基于测序的生物学研究中的这一数据挑战需要一种具有高准确性、灵活性和速度的序列比对工具。
大多数序列比对工具基于将序列分割为固定长度为k的seed子序列(k-mers)的搜索算法,并在参考数据库中寻找这些k-mers的现有出现位置,例如用于序列比对的Minimap、相似性搜索的Blast和宏基因组分类注释的Kraken。k-mer在序列比对中的使用使得处理大规模测序数据集成为可能,并推动了包括现代分子生物学和现代进化生物学在内的基于测序的生物学领域的发展。
然而,k-mer算法的已知挑战是需要选择一个合适的k-mer长度,而k-mer长度的选择对序列比对的效果有显著影响——固定的k-mer长度可能同时过长,无法匹配突变密集的基因组区域;同时又可能过短,无法区分重复基因组区域。例如,在一个突变密集区域(如每20个碱基对有1个突变)中,较长的k-mer(如30bp)可能无法找到匹配位置;而对于长度超过50bp的重复区域,较短的k-mer(如10bp)可能会找到许多匹配位置。不同物种的基因组多样性表明,基于固定大小k-mers的算法在宏基因组研究中的灵活性较差。
这一问题如图1a所示。在查询序列中的突变区域(常见于变异密集区域),k-mers(图1b中的k-mers 1-3)无法在参考基因组中找到匹配位置;这些k-mers过长,无法避免错配。相比之下,来自查询序列中重复区域的k-mers(如图1b中的k-mers 4和5)在参考基因组中会匹配到多个位置。因此,仅有12个匹配中的2个(蓝色高亮部分)对最佳比对有贡献(图1a),且较短的k-mer长度(图1b中的灰色匹配)效率较低。因此,对于高突变区域,较短的k-mer长度更有助于避免错配,而对于重复区域,较长的k-mer长度则更适合。然而,由于一个基因组中既包含突变区域又包含重复区域,固定大小的k-mers无法同时对所有基因组区域实现最优效果。尽管现有比对工具允许自定义k-mer大小以解决不同基因组的问题,但即使是最优化的k-mer长度也无法解决这一问题,因为所有基于k-mer的比对工具在一次运行中使用固定大小。
图1 | 一个示例展示了基于x-mer的算法(c)如何比基于k-mer的算法(b)更有效地处理低复杂度重复区域和高复杂度变异密集区域,以及使用gapped x-mer的算法(d)如何进一步改进变异密集区域的处理。
(a) 待比对序列比对到示例参考中的最优比对结果以及一些候选次优比对结果。参考序列被标记为变异密集区域和重复区域。
(b–d) 用蓝色高亮的k-mers/x-mers/gapped x-mers表示与最优比对一致的匹配。用灰色高亮的k-mers/x-mers/gapped x-mers表示不一致但支持次优比对的匹配。算法随后会识别出所有唯一位置(查询在参考中的位置),这些位置是任一k-mer/x-mer/gapped x-mer的映射结果。对于每个偏移位置,算法可能尝试将k-mer/x-mer/gapped x-mer的匹配扩展到邻近的碱基对,以在该位置为查询生成比对结果。连接k-mers/x-mers/gapped x-mers的每条虚线表示一个可能发生这种扩展的唯一查询偏移位置。
(b) k-mers 4和5在参考序列中匹配了多个偏移位置,而k-mers 1–3没有匹配。
(c) x-mers 1–3在参考序列中各匹配了1–2个偏移位置。
(d) gapped x-mers 1–2在参考序列中匹配了一个偏移位置。每个gapped x-mer中的“_”表示一个1 bp的缺口。
一种解决这些限制的方法是引入具有不同大小的k-mers,我们称之为x-mers。x-mer方法在查询序列中的突变区域缩短x-mer大小,允许绕过点突变(如图1c中的x-mers 1-2),避免了k-mers 1-3(图1b)遇到的问题。对于查询序列中的重复区域,我们延长x-mer大小以区分重复区域(如图1c中的x-mer 3)。在这一示例中,所有x-mers在参考基因组中匹配到1-2个位置,而4个匹配中的3个(蓝色高亮部分)对最佳比对有贡献(图1a),相比于使用k-mers找到的12个匹配中的2个(图1b)。此外,在最佳偏移位置,x-mers可以覆盖突变密集区域中的8bp中的6bp,而k-mers则无法覆盖突变密集区域。这些结果表明,x-mers可以更精确地描述适当的比对并利用查询中的更多信息。然而,适应密集点突变的短x-mers(如x-mers 1和2)更可能对多个次优偏移位置进行比对,从而降低其特异性。
虽然可变长度seed(x-mers)已被LAST用于局部比对,并且间隔seed首先由PatternHunter引入,后来被许多比对工具广泛采用,如SHRiMP2,我们引入了动态整合动态大小和间隔的k-mers的概念,我们称之为带缺口的x-mers。我们认为,带缺口的x-mer最有趣的贡献在于能够同时提供可变长度和间隔。这可以通过仔细选择存储哪些seed来实现,而不需要存储每个可能的seed,这种方法类似于minimizers。该算法在突变密集区域无需缩短k-mer大小,而是生成可以跳过不同突变组合的间隔(如图1d中的带缺口的x-mer 1),从而保持更高的特异性。因此,在这一示例中,每个带缺口的x-mer在参考基因组中恰好匹配一个位置,从而在一次搜索中识别出最佳比对。这种在精确性、特异性和覆盖范围方面的改进使带缺口的x-mers能够更高效、更快速地进行比对。
我们将带缺口的x-mer概念应用于短读比对并开发了一种短读比对工具,X-Mapper(https://github.com/mathjeff/Mapper)。总体而言,X-Mapper是一个高精度的序列比对工具,与其他比对工具相比,它减少了次优比对的数量。X-Mapper已在多种微生物组测序样本中进行了测试,包括全基因组测序数据和宏基因组数据,因为我们预计基于带缺口的x-mer的算法能更好地适应不同物种中包含突变区域和重复区域的多样化组合,而不是基于k-mer的算法。X-Mapper的高准确性、灵活性和速度可以为依赖短读长测序的多种生物学研究提供帮助。经过进一步优化后,X-Mapper还可能对人类基因组等更大规模的基因组研究有所帮助。我们设想,带缺口的x-mer算法的新概念可以激发新的seed生成方法,并有助于改善其他基于k-mer和x-mer的生物信息学应用,例如相似性搜索(BLAST,Diamond)、宏基因组中的分类赋值(Kraken)、多序列比对(MUSCLE,MAFFT)和基因组组装(SPAdes,IDBA)。
结果Results
X-Mapper算法
在本研究中,我们设计了一种基于gapped x-mer的算法,利用各种可能大小的gapped x-mer,并基于该算法开发了一种新的短读比对工具——X-Mapper。X-Mapper首先从1个碱基对开始构建x-mer的金字塔,并根据需要扩展到整个序列的长度(如图2,详细信息见“Methods”)。构建x-mer后,X-Mapper生成gapped x-mer,其长度基于参考序列长度进行调整,通过添加额外的x个碱基对和0到2个碱基对(具体为hashcode模3的伪随机数)以避免因某些碱基对数跳过gapped x-mer的生成。这些大约x个碱基对均匀分配为gap和extension碱基对,形成约x/2个碱基对的gap和x/2个碱基对的extension。每个x-mer根据自身内容选择gap的方向。随后,X-Mapper为每个gapped x-mer生成一个hashcode。参考序列中的gapped x-mer随后被保存到hashtable中。这些gapped x-mer由相同的算法针对参考基因组和每个查询序列生成,以便在hashtable中找到匹配项。
图2 | X-Mapper 的算法。X-Mapper 通过构建 x-mer 的金字塔开始。
X-Mapper在不同复杂性样本中的比对准确性
为了研究 X-Mapper 的比对准确性,我们首先将 X-Mapper 的比对penalty(alignment score)与基于 k-mer 的算法(如 Strobealign 和 Minimap2)以及基于 x-mer 的方法(如 Bowtie2、BWA(此研究中指 BWA MEM)和 LAST)进行了比较。由于每个比对工具在处理间隔的penalty公式上有所不同,跨工具比较双端比对结果的penalty并不明确。因此,在本次评估中,我们测试了单端测序数据的比对,比较了不同比对工具对相同read的比对结果。具体而言,我们检查了每个比对工具在相同penalty设置下将相同read比对到相同参考基因组的表现(图 3a)。我们测试了两个场景:一个是人类肠道微生物组的宏基因组比对到其自身组装,另一个是人类转录组数据比对到完整的人类基因组参考。测试中使用了四种penalty设置,包括 X-Mapper、Bowtie2、Minimap2 和 BWA 的默认设置。我们比较了这四种比对工具报告的相同read的比对结果,并将penalty最低(alignment score最高)的比对视为最佳比对。在这里,我们使用独立脚本重新计算了比对penalty,该脚本基于每个比对工具报告的 CIGAR 字符串和位置确定匹配数、不匹配数、模糊匹配数、间隙开启和间隙扩展。penalty较高(alignment score较低)的比对被定义为次优比对。随后,我们计算了每个比对工具报告的次优比对的read百分比。
图3 | 评估X-Mapper在不同复杂性测序样本中的比对准确性和一致性。
显示的penalty值是单个read的示例penalty值。
a , 比对准确性通过比较相同read在不同比对工具下比对到相同参考基因组的结果,在相同penalty(alignment score)设置下进行评估。具有最低penalty值(较高alignment score)的比对被认为是最优的,而所有其他具有更高penalty值(较低alignment score)的比对被认为是次优的。
b , 比对一致性通过比较同一比对工具在一个参考(Assembly A1)和一个复杂参考(多个基因组)上的相同read比对结果来评估。如果满足以下条件之一,则认为比对结果一致:(1) 在简单参考中的最优比对结果在复杂参考中也被报告;(2) 复杂参考比对结果优于简单参考(具有更低的penalty值或更高的alignment score),主要是由于组装错误造成的。然后对比对一致性在不同比对工具之间进行比较。
我们的结果表明,X-Mapper 的比对准确性高于所有其他比对工具,这通过较低的次优比对率得以体现(图 4)。具体而言,在比对一个人类肠道微生物组宏基因组时(图 4a,“all alignments”),X-Mapper 在 930 万read中显示了 0.05% 的次优比对率,这比 Strobealign (1.63%)、LAST (0.29%)、Bowtie2 (0.53%)、Minimap2 (0.44%) 和 BWA (0.29%) 在penalty设置 #1 下的比对率低 6–34 倍。在排除 Minimap2 报告中间软剪切的比对(图 4a,“alignments without middle soft clips”)后,X-Mapper 的次优比对率仍保持在 0.05%,比 Strobealign (1.55%)、LAST (0.26%)、Bowtie2 (0.48%)、Minimap2 (0.37%) 和 BWA (0.23%) 低 5–33 倍。为 Bowtie2 配置允许 seed 错配(“-N 1”)将其准确性从 0.48% 提高至 0.27%,这一趋势在所有penalty设置中均有所观察。因此,在后续分析中,我们使用了这一允许 seed 错配的 Bowtie2 配置。对于 Strobealign,我们在默认penalty和相对于不同alignment score调整后的penalty下进行了测试(详见“Methods”)。然而,其次优比对率仍然相对较高,最低为 0.96%。
在将人类转录组数据比对到完整的人类基因组参考时,我们发现 X-Mapper 的次优比对率为 0.48%,比 Strobealign (10.2%)、LAST (5.7%)、允许 seed 错配的 Bowtie2 (5.3%)、Minimap2 (6.6%) 和 BWA (5.7%) 在penalty设置 #1 下的比对率低 11–24 倍(图 4b,“all alignments”)。在移除 Minimap2 报告中间软剪切的比对后,X-Mapper 的次优比对率为 0.25%,比 Strobealign (4.6%)、LAST (2.4%)、允许 seed 错配的 Bowtie2 (1.58%)、Minimap2 (1.0%) 和 BWA (2.1%) 低 4–18 倍(图 4b,“alignments without middle soft clips”)。
图4 | X-Mapper 在应用于(a)细菌宏基因组和(b)人类转录组时,与其他比对工具相比,表现出更高的比对准确性和更低的次优比对率(报告次优比对的reads百分比)。最佳比对被定义为所有比对工具中报告的penalty最低(alignment score最高)的比对。每个样本在四种不同的penalty设置下进行了测试,包括:(1) Bowtie2 的默认设置,(2) BWA(本研究指BWA MEM),(3) X-Mapper,以及(4) Minimap2(x轴)。对于如Minimap2这样的局部比对工具,报告在contigs中部出现soft clips的reads被从下游分析中移除(“alignments without middle soft clips”)。
X-Mapper 在不同复杂性样本中的比对一致性
除了低比对 penalty,确保 aligner 在处理简单样品和复杂样品(多物种宏基因组)时,能够报告一致的比对结果也至关重要(图 3b)。为了测试比对一致性,我们将一个 B. fragilis 的 WGS 数据集(320 万成对端read,150 bp)比对到其自身的组装(表示简单参考)。我们还将该 WGS 数据集比对到包含 88 个人类肠道微生物组菌株(总计 441.4 Mbp)的基因组集合(表示复杂参考)。此集合包括:(1) B. fragilis 的 WGS 组装(Assembly1);(2) 5 个其他 B. fragilis 菌株参考;(3) 55 个非 B. fragilis 的 Bacteroides 物种;(4) 27 个非 Bacteroides 物种。在比对过程中,这 88 个基因组为 WGS read竞争。
我们发现 Bowtie2、Minimap2、BWA 和 Strobealign 的默认设置是报告一条比对结果。因此,我们还在这些 aligner 上启用了报告多条(例如 100 条)或所有比对结果的设置。结果是,启用这种“all”设置后通常会报告次优比对。虽然报告次优比对在技术上没有错,但对用户来说效率较低,因为需要在下游分析前筛选比对结果。此外,这些次级比对与最佳比对相比,通常多出 4–5 个点突变,即使最佳比对是完美匹配的。报告这些次优比对似乎没有帮助。因此,我们还计算了比对到复杂参考时报告次优比对的read百分比。
我们发现,与其他 aligner 相比,X-Mapper 的不一致性和次优比对率最低。例如,在 penalty 设置#1下,X-Mapper 的不一致率为 0.09%(在 3,208,556 个比对read中),次优比对率为 1.2e − 06(图 6)。Bowtie2(“all”模式报告所有比对)、Strobealign(“all”模式)和 LAST 的不一致率分别为 0.31%(在 3,209,652 个比对read中)、2.40%(在 3,210,658 个比对read中)和 3.25%(在 3,211,829 个比对read中)。然而,这些模式伴随着高比例的次优比对read,分别为 73.5%、53.4% 和 6.5%。其他 aligner 表现出 562–579 倍更高的不一致率,包括 Minimap2“all”模式的 52.0%(在 3,215,249 个比对read中)、BWA“all”模式的 51.9%(在 3,223,873 个比对read中)、Strobealign“default”模式的 51.4%(在 3,210,436 个比对read中)、Bowtie2“default”模式的 52.1%(在 3,209,422 个比对read中)、Minimap2“default”模式的 52.0%(在 3,215,249 个比对read中)和 BWA“default”模式的 52.0%(在 3,223,843 个比对read中)。此外,这些 aligner 的次优比对率比 X-Mapper 高 4–133 倍,分别为 1.6e − 04(Minimap2“all”模式)、1.6e − 04(BWA“all”模式)、3.9e − 05(Strobealign“default”模式)、5.6e − 06(Bowtie2“default”模式)、1.6e − 04(Minimap2“default”模式)和 6.5e − 06(BWA“default”模式)。
我们在 Bowtie2、Minimap2、BWA 和 X-Mapper 的默认设置下,以及多个 penalty 设置下测试了这些 aligner,结果始终发现 X-Mapper 提供了最一致的比对结果,并且次优比对的比例最少。这些观察结果表明,X-Mapper 在对复杂参考(如微生物宏基因组)进行比对时更加一致且更易解释。如果用户对次优比对感兴趣,可以通过 X-Mapper 中的单独参数(–max-penalty-span)进行控制,此处未测试。
图6 | X-Mapper在比对到复杂参考基因组时,表现出最低的比对一致性问题(reads中报告为不一致比对的比例)和最少的次优比对(reads中报告为次优比对的比例)。
比对一致性通过将一个全基因组测序(WGS)样本比对到88个细菌基因组(代表复杂参考)和比对到自身组装(代表简单参考)的结果进行比较来测试。若满足以下条件之一,则认为比对是一致的:(1) 对简单参考的最优比对在比对到复杂参考时也被报告,或(2) 复杂参考产生了比简单参考更好的比对(具有更低的penalty或更高的alignment score)。Bowtie2、Minimap2和BWA在其默认设置以及报告多个(甚至全部)比对的设置(标记为“all”)下运行。由于“all”设置通常会报告次优比对,因此还评估了每个比对工具中报告次优比对的reads比例。比对一致性使用四种不同的penalty设置进行了测试,包括(1) Bowtie2的默认设置,(2) BWA的默认设置,(3) X-Mapper的默认设置,以及(4) Minimap2的默认设置。分析中还测试了Bowtie2的seed mismatch配置。
讨论Discussion
X-Mapper 被设计为在过滤比对结果方面提供更大的灵活性,并生成更适合下游分析的输出,例如过滤后的 VCF 文件。在之前的工作中,我们发现诸如 SAMtools 和 BCFtools 等工具对比对结果的总结可能不完整,这显著影响了包括遗传变异调用在内的下游分析的准确性。具体来说,BCFtools(“bcftools call”命令)经常在重叠较少的配对末端read中排除 indels。这种 indels 的排除会导致主要变异的低估,从而提高次要变异的频率,进而导致点突变识别中的假阳性。因此,X-Mapper 包含一些特定设置(例如 –max-penalty-span),允许用户控制在指定的 penalty 范围内包含次优比对,或者完全排除它们。
为了改进内存使用,我们开发了 X-Mapper Next,使得在具有 16 GB RAM 的本地 PC 上对人类基因组进行比对成为可能,而 X-Mapper 的内存使用为 48 GB。对于细菌基因组,X-Mapper Next 将 441.4 Mbp 参考基因组的内存使用从 15–20 GB 减少到 3 GB,将 99.0 Mbp 参考基因组的内存使用从 3 GB 减少到 1 GB。
X-Mapper Next 在测试人类肠道微生物组测序数据时表现出相当的准确性、一致性和速度。次优比对率从0.05%降低到0.04%,read比对运行时间从41.17秒增加到59.50秒(30线程,penalty setting #1)。与此同时,不一致率从0.09%降低到0.08%,read比对运行时间也从100.99秒减少到84.79秒(30线程,penalty setting #1)。此外,X-Mapper Next 与 X-Mapper 相比,保持了高效的参考索引运行时间,在441.4 Mbp参考上为26.54秒,而 X-Mapper 为28.18秒(30线程);在99.0 Mbp参考上为10.06秒,而 X-Mapper 为7.21秒(30线程)。我们将继续改进 X-Mapper Next,并在未来对人类基因组进行更全面的测试。
作者简介
南洋理工大学章安妮组高级研究科学家Jeffry Gaston是本文的第一作者,麻省理工教授Eric Alm和南洋理工大学助理教授章安妮为本文的通讯作者。
Jeffry Gaston(第一作者)
Jeffry Gaston,前谷歌工程师,现任南洋理工大学章安妮组高级研究科学家。
Eric Alm (通讯作者)
Eric Alm,麻省理工教授。他的研究团队包括计算/理论和实验方法,用于理解和工程化人类微生物组。他们的研究重点是将基础科学发现迅速转化为临床实践,以改善患者的治疗效果。
章安妮 (通讯作者)
章安妮,南洋理工大学助理教授,博士生导师。她是一名拥有10年经验的生物信息学家和微生物学家,在算法开发和微生物组基因组学领域拥有丰富的研究成果和坚实的发表记录。她组建了一支跨学科团队,包括一名软件工程师和四名博士生:一名生物信息学家、一名化学工程师、一名人工智能专家以及一名微生物学家。她的实验室专注于设计针对生物数据的AI模型,以揭示微生物机制,使微生物对社会更有价值。
宏基因组推荐
1月10-12日,单菌基因组组装、注释、遗传表征、分子分型、系统进化和传播溯源
2月21-23日,家系、肿瘤临床基因组/外显子组数据分析
3月21-23日,高级转录组分析和R语言数据可视化
3月28-30日,第二届全国基因组信息学大会
4月11-13日,微生物组-扩增子16S分析
5月11-13日,微生物组-宏基因组分析
本公众号现全面开放投稿,希望文章作者讲出自己的科研故事,分享论文的精华与亮点。
投稿请联系小编(-genomics)
iMeta高引文章 fastp 复杂热图 ggtree 绘图imageGP 网络iNAP
iMeta网页工具 代谢组MetOrigin 美吉云乳酸化预测DeepKla
iMeta综述 肠菌菌群 植物菌群 口腔菌群 蛋白质结构预测
10000+:菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature
为鼓励读者交流快速解决科研困难,我们建立了“宏基因组”讨论群,己有国内外6000+ 科研人员加入。请添加主编微信meta-genomics带你入群,务必备注“姓名-单位-研究方向-职称/年级”。高级职称请注明身份,另有海内外微生物PI群供大佬合作交流。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍未解决群内讨论,问题不私聊,帮助同行。
来源:微生物组