南洋理工章安妮组Genome Biology：开发新的序列比对生信工具分析人类和微生物基因比现有工具准确11-24倍

成果简介

序列比对是许多生物信息学分析的基础。许多比对工具通过将序列分割为连续的固定长度片段（称为 k-mers ）开始。使用较长且独特的 seed 比对速度更快，而使用较短且避免突变的 seed 则更准确。在此，我们介绍了 X-Mapper ，旨在通过包含缺口的动态长度 seed （称为 gapped x-mers ）实现高速和高精度比对。分析人类参考序列时，我们观察到次优比对减少了 11–24 倍，而在细菌参考序列中不一致性降低了 3–579 倍，相较于其他比对工具，在非目标菌株和物种上的比对分别改进了 53% 和 30% 。其他基于 seed 的分析算法也可能受益于 gapped x-mers 。

短读长测序应用广泛，例如微生物组的宏基因组学研究，数据分析要求研究人员处理海量数据集并应对多样化的基因组。在处理短序列的所有步骤中，序列比对是计算最密集的任务之一，比对的质量直接影响下游分析的效率和准确性。因此，基于测序的生物学研究中的这一数据挑战需要一种具有高准确性、灵活性和速度的序列比对工具。

大多数序列比对工具基于将序列分割为固定长度为 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 ）。

图文导读

X-Mapper 算法

X ‑ Mapper algorithm

在本研究中，我们设计了一种基于 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 在不同复杂性样本中的比对准确性

Alignment accuracy of X‑Mapper in samples with various complexities

为了研究 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 在不同复杂性样本中的比对一致性

Alignment consistency of X‑Mapper in samples with various complexities

除了低比对 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 配置。

‘

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 ，并在未来对人类基因组进行更全面的测试。