Systematic benchmarking of small variant calling pipelines for long-read RNA sequencing data

技术摘要：长读段 RNA 测序数据中小型变异检测流水线的系统性基准测试

问题陈述
长读段 RNA 测序（lrRNA-seq）能够实现转录本解析的变异检测，将遗传变化与特定异构体联系起来，并促进对等位基因特异性事件和剪接的研究。然而，针对 lrRNA-seq 的小型变异（SNVs 和 indels）检测流水线的系统性且中立的评估仍然有限。现有工具在不同测序技术（Oxford Nanopore Technologies [ONT] 和 Pacific Biosciences [PacBio]）、比对策略、变异检测器选择、基因组上下文以及下游单倍型相位化（haplotype phasing）方面的性能尚未得到充分理解。此外，对于专门为 lrRNA-seq 优化的变异检测器而言，比对文件转换（将剪接的 RNA 比对转换为连续的类 DNA 表示）是否仍然必要，目前尚不明确。

方法论
作者使用 Genome in a Bottle (GIAB) 数据集进行了系统性基准测试，该数据集包含三种细胞系（HG002、HG004、HG005），分别采用四种 ONT 文库制备方案（cDNA R9.4.1、cDNA R10.4.1、dRNA002、dRNA004）和两种 PacBio 方案（Iso-Seq 和 Mas-Seq/Kinnex）进行测序。研究还扩展到了最近的 LongBench 数据集（H211 和 H526 肺癌细胞系）。

基准测试工作流包括：

1. 比对： 使用 minimap2 和 pbmm2（一种针对 PacBio 优化的包装器）进行序列比对。研究评估了比对转换（使用 SplitCigarReads 来模拟连续的 DNA 比对）与使用原始剪接比对对结果的影响。
2. 变异检测： 评估了五种 lrRNA-seq 特有的检测器：Clair3-RNA、DeepVariant（针对 PacBio 使用 MASSEQ 模型，针对 ONT 使用 ONT DNA 模型）、isoLASER、longcallR 以及 longcallR-nn。同时纳入了 GATK HaplotypeCaller 作为短读段基准。
3. 单倍型相位化： 使用 WhatsHap、LongPhase、HapCUT2、HiPhase 以及 isoLASER 和 longcallR 的原生相位化输出来评估相位化性能。
4. 评估指标： 通过在高置信度区域内针对真实值（ground-truth）调用的精确度（precision）、召回率（recall）和 F1 分数进行量化，并根据测序深度、变异类型（SNV/indel）和基因组上下文（同聚物、串联重复序列、GC 含量）进行分层。相位化通过转换错误率（switch error rate）、相位覆盖度（phasing coverage）和相位块 N50 进行评估。

核心贡献与结果

• 测序特性是主要决定因素： 研究表明，测序特性（化学性质、错误率、吞吐量和可比对性）是变异检测性能的最强决定因素，其对 F1 分数和真阳性计数解释的方差大于变异检测器的选择。PacBio 数据集（尤其是 Mas-Seq）由于更高的测序准确度和更多的可调用碱基，表现持续优于 ONT 数据集。在 ONT 化学性质中，dRNA004 和 cDNA R10 表现最好，而 dRNA002 表现最差。
• 变异检测器性能：
  • Clair3-RNA 在针对两种 ONT 和 PacBio 化学性质的 SNV 和 indel 检测中，展现出了最广泛且最稳健的适用性。
  • DeepVariant 在 PacBio 数据上表现强劲，尤其是在 indel 检测方面，但在 ONT 数据上表现较差，原因是缺乏专门的 ONT RNA-seq 模型。
  • longcallR 在无需依赖预训练模型的情况下，在 ONT 数据的 SNV 检测方面具有竞争力，并提供直接的相位化输出。
  • isoLASER 显示出高精确度但较低的召回率，特别是在 ONT 数据上，且在变异检测方面较为保守。
  • longcallR-nn 在特定数据集上表现出竞争力的性能，但在 dRNA002 上表现挣扎。
• 比对策略：
  • 比对器选择： minimap2 通常比 pbmm2 产生更高的 F1 分数，因为 pbmm2 往往会引入更多的假阳性，特别是在剪接区域，这可能是由于参数预设的差异以及缺乏剪接位点注释输入所致。
  • BAM 转换： 对于 lrRNA-seq 特有的检测器，比对转换（在外显子-内含子接点处拆分读取段）并不会提高性能，在某些情况下反而会降低准确度，因为它破坏了 RNA 特有的比对特征（如剪接感知 CIGAR 图案）和单倍型信息。仅当在 PacBio 数据集上使用 DeepVariant 进行 indel 检测时，转换才是受益的。
• 基因组上下文与 RNA 编辑：
  • 性能在挑战性的基因组上下文中会发生退化，特别是在同聚物（尤其是 A/T）和极端 GC 区域。Indel 检测对同聚物长度的敏感度高于 SNV 检测。
  • RNA 编辑（A-to-I）导致了假阳性。集成 RNA 编辑数据库（例如通过 REDIportal 对 Clair3-RNA 进行位点标记）或使用迭代相位化（longcallR）的检测器能更好地区分真实的变异与编辑事件。
• 单倍型相位化： 没有一种单一的相位化方法能在所有指标上都达到最优。WhatsHap 倾向于相位覆盖度和连续性（更长的相位块），但以较高的转换错误率为代价；而 HapCUT2 则优先考虑准确性，但覆盖度略低。longcallR 为联合变异检测和相位化提供了强有力的替代方案。
• 泛化能力： 在 GIAB 数据集中观察到的性能趋势在 LongBench 癌症细胞系数据集中也是一致的，这进一步证实了测序特性和检测器相对性能层级的支配地位。

意义
本文为 lrRNA-seq 变异检测建立了基础性的基准测试，强调虽然工具选择很重要，但底层的测序技术和文库制备是成功的关键因素。它为构建流水线提供了具体的指导建议：

1. 推荐将 Clair3-RNA 作为通用 lrRNA-seq 变异检测中最具通用性的工具。
2. DeepVariant 是处理 PacBio 数据（尤其是 indel）的强力选择。
3. 对于 lrRNA-seq 特有的检测器，应避免使用比对转换，以保留长程信息，除非是在 PacBio 数据集上使用 DeepVariant。
4. WhatsHap 和 HapCUT2 在下游相位化中提供互补的权衡，而 longcallR 适用于集成的检测与相位化。

研究结论指出，未来的灵敏度提升可能需要针对性的富集策略，以克服表达依赖性的限制，并且需要进一步的工作来对转录本水平的变异分配和单细胞应用进行基准测试。