以乌苏图病毒为例验证全基因组纳米孔测序

摘要

我们此前已验证了一种基于扩增子的乌苏图病毒（USUV）全基因组测序方案，该方案在纳米孔测序平台上进行。本文中，我们对所使用的方法进行了更详细的描述，并确定了纳米孔R10流通池的错误率。

摘要

全基因组测序可用于病毒暴发的特征分析与溯源。已有多种基于纳米孔测序平台的全基因组测序方案被用于不同病毒的研究。这些方法采用基于重叠扩增子的策略，可靶向特定病毒或一组遗传关系相近的病毒。除了确认病毒的存在外，测序还可用于基因组流行病学研究，追踪病毒传播路径，并揭示其来源、宿主库及传播方式。对于此类应用，理解所用测序平台相关错误率的潜在影响至关重要。在临床和公共卫生领域的常规应用中，每次对方案进行重要修改时，都必须对这些影响进行记录。此前，已通过与Illumina测序直接比较的方法，验证了在纳米孔测序平台（R9.4流动池）上进行乌苏图病毒全基因组测序的方案。本文以R10流动池与Illumina测序的比较为例，描述了确定所需测序读段覆盖深度的方法。

引言

第三代测序技术的快速发展使我们能够在病毒暴发期间接近实时地进行测序。及时获得的遗传信息有助于确定病毒病原体的来源和进化过程。然而，下一代测序领域的金标准仍然是第二代测序仪。这些技术依赖于特定且耗时的步骤，例如乳液PCR中的克隆扩增或克隆桥接扩增。第三代测序仪成本更低、可手持操作，并且具有简化的文库制备方法。特别是测序设备的小型化和较低的购买价格，使其成为可部署、适用于现场测序的理想选择。例如，在塞拉利昂的埃博拉病毒暴发期间以及在巴西正在进行的虫媒病毒暴发调查中已体现出其应用价值^1,2,3。然而，据报道，纳米孔测序存在较高的错误率⁴，这可能会限制其应用范围。

纳米孔测序技术正在快速发展，市场上定期推出新产品。例如，1D²试剂盒能够对DNA分子的两条链进行测序，从而提高碱基识别的准确性⁵，以及R10流动池的开发，可在纳米孔内两个不同时间点测量电流变化⁶。此外，生物信息学工具的改进，如碱基识别算法的优化，也将提升碱基识别的准确性⁷。其中最常用的碱基识别软件之一（如Albacore）在9个月内至少更新了12次⁵。最近，制造商还发布了一种名为flip-flop的新型碱基识别器，并已集成至默认的纳米孔测序软件中⁸。这些改进的综合效应将带来更准确的测序结果，并降低纳米孔测序仪的错误率。

乌苏图病毒（USUV）是一种通过蚊子传播的节肢动物传播病毒，属于黄病毒科（Flaviviridae），其基因组为约11,000个核苷酸的单股正链RNA。USUV主要影响灰林鸮和乌鸫^9,10，尽管其他鸟类物种也易感USUV感染¹¹。最近，USUV也在啮齿类动物和鼩鼱中被发现，但它们在病毒传播中的潜在作用尚不清楚¹²。在人类中，已在献血者中发现无症状感染^13,14,15,16，同时也有报道指出USUV感染与脑炎或脑膜脑炎相关^17,18。在荷兰，USUV于2016年首次在野生鸟类中被检测到¹⁰，并于2018年在无症状献血者中发现¹⁴。自首次发现USUV以来，此后数年均报告了疫情暴发，目前正在进行包括全基因组测序在内的监测工作，以追踪该节肢动物传播病毒在既往无免疫人群中出现及传播的情况。

与其他病毒（如埃博拉病毒、寨卡病毒和黄热病毒）已报道的方法类似^3,19,20，我们开发了一套引物用于对USUV全长基因组进行测序²¹。这种基于聚合酶链式反应（PCR）的方法能够从宿主污染程度较高的样本（如脑组织样本）中回收全长USUV基因组，适用于Ct值高达约32的样本。基于扩增子的测序方法具有比宏基因组测序更高的灵敏度和特异性。然而，该方法的局限性在于：为覆盖所有毒株，所设计引物需基于序列相似性，且引物设计依赖于当前对病毒多样性的认知水平。

鉴于第三代测序技术的持续发展与改进，有必要定期评估测序仪的错误率。本文以 Usutu 病毒（USUV）为例，介绍一种将纳米孔测序结果直接与 Illumina 测序结果进行比对以评估其性能的方法。该方法应用于使用最新 R10 流通池生成的序列，并采用最新版本的翻转-翻转（flip-flop）碱基识别算法进行碱基识别。

方案

注意：将使用的软件工具列表：usearch v11.0.667；muscle v3.8.1551；porechop 0.2.4；cutadapt 2.5；minimap2 2.16-r922；samtools 1.9；trimmomatic 0.39；bbmap 38.33；spades v3.13.1；kma-1.2.8

1. 引物设计

1. 首先从公共或私人数据库中下载或获取一组相关的参考全基因组序列。例如，从NCBI数据库中获取所有全长的USUV基因组（taxid64286）²²。USUV编码的基因组长度约为11,000个核苷酸，因此仅获取序列长度在8,000至12,000个核苷酸之间的序列。使用以下搜索条目进行检索：
   - taxid64286[Organism:noexp] AND 8000[SLEN]:12000[SLEN]。
   1. 点击 Send to | Complete Record | File；选择格式为FASTA并生成文件。
2. 为减少参考序列的数量，需从数据集中去除重复序列或核苷酸序列相似性超过99%的序列。可使用usearch工具中的cluster_fast选项完成此操作²³。在命令行中输入：
   - usearch -cluster_fast All_USUV.fasta -id 0.99 -centroids All_USUV_dedup.fasta
3. 为设计引物，需对序列进行比对，此步骤使用MUSCLE软件完成²⁴。在命令行中输入：
   - muscle -in All_USUV_dedup.fasta -out All_USUV_dedup_aligned.fasta -log log_muscle.txt
   注意：必须手动检查比对结果以发现可能的错误。如有需要，可手动修正错误，并根据大多数全基因组序列的长度对两端进行裁剪。
4. 使用Primal工具初步筛选可用于全长扩增子测序的引物¹⁹。将比对好的序列上传至Primal网站（http://primal.zibraproject.org/），选择理想的扩增子长度以及不同扩增子之间的重叠长度。访问primal.zibraproject.org，填写方案名称，上传比对后的FASTA文件，设置扩增子长度、重叠大小，然后生成引物方案。
5. 对所有可获得的完整USUV序列（非缩减或去重后的序列集）进行比对。在命令行中输入：
   - muscle -in All_USUV.fasta -out All_USUV_aligned.fasta -log log_muscle.txt
   注意：将生成的引物比对到完整的序列比对结果中（不要使用去重后的比对结果），手动修正错误，并最多包含5个简并引物位点。

2. 多重 PCR

1. 使用设计的引物进行多重PCR，并进行纳米孔测序和Illumina测序。针对USUV的多重PCR按照先前所述方法进行^19,21。
2. 使用flip-flop版本3.0.6.6+9999d81进行碱基识别。

3. 利用纳米孔数据生成共识序列的数据分析

1. 可在单次纳米孔测序运行中对多个样本进行多重检测。完成测序后，对纳米孔数据进行解复用，可使用 Porechop²⁵ 实现。为防止污染并提高准确性，需使用 require_two_barcodes 参数。在命令行中输入：
   - porechop -i Run_USUV.fastq -o Run_USUV_demultiplex --require_two_barcodes
2. 解复用后，使用 cutadapt²⁶ 去除引物序列（正向和反向序列均在文件 Primers_Usutu.fasta 中标明）。此外，去除长度短于 75 个核苷酸的序列。必须去除引物序列，因其可能在共识序列中引入人为偏差。在命令行中输入：
   - cutadapt -b file:Primers_USUV.fasta -o BC01_trimmed.fastq BC01.fastq -m 75
3. 可使用 minimap2²⁷ 将解复用后的测序读段比对至一组不同的参考毒株，并使用 samtools²⁸ 生成共识序列。以下示例展示了一个样本（BC01）基于参考序列的比对及共识序列生成流程。在命令行中输入：
   - minimap2 -ax map-ont Random_Refs_USUV.fasta BC01_trimmed.fastq > BC01.bam
   - samtools sort BC01.bam > BC01_sorted.bam
   - bcftools mpileup -Ou -f Random_Refs_USUV.fasta BC01_sorted.bam | bcftools call -mv -Oz -o BC01.vcf.gz
   - bcftools index BC01.vcf.gz
   - cat Random_Refs_USUV.fasta | bcftools consensus BC01.vcf.gz > BC01_consensus.fasta
4. 在进行基于参考序列的比对时，必须使用亲缘关系密切的参考序列。因此，应对生成的共识序列进行 BlastN 搜索，以确定最接近的参考毒株。随后，使用该最近的参考毒株作为参考序列，重复进行基于参考的比对（步骤 3.3 和 3.4）。在命令行中输入：
   - minimap2 -ax map-ont Ref_USUV_BC01.fasta BC01_trimmed.fastq > BC01_ref.bam
   - samtools sort BC01_ref.bam > BC01_sorted_ref.bam
   - bcftools mpileup -Ou -f Ref_USUV_BC01.fasta BC01_sorted_ref.bam | bcftools call -mv -Oz -o BC01_ref.vcf.gz
   - bcftools index BC01_ref.vcf.gz
   - cat Ref_USUV_BC01.fasta | bcftools consensus BC01_ref.vcf.gz > BC01_ref_consensus.fasta

4. Illumina 数据分析

1. 这些序列在测序后会自动进行解复用。可使用 Trimmomatic 对读段进行质量控制。²⁹对于成对末端Illumina测序数据，使用常用的截断值——中位PHRED得分为33，并设置最小读长为75，以获得准确且高质量的测序读段。在命令行中输入：
   - trimmomatic PE -phred33 9_S9_L001_R1_001.fastq.gz 9_S9_L001_R2_001.fastq.gz 9_1P.fastq 9_1U.fastq 9_2P.fastq 9_2U.fastq LEADING:3 TRAILING:3 SLIDINGWINDOW:3:15 MINLEN:75
2. 使用 cutadapt 去除引物序列（引物序列信息见文件 Primers_Usutu.fasta 中的正向和反向序列），以避免引入人为偏差²⁶此外，使用以下命令删除长度短于75个核苷酸的序列。在命令行中输入：
   - cutadapt -b o 9_1P_trimmed.fastq -p 9_2P_trimmed.fastq 9_1P.fastq 9_2P.fastq -m 75
3. 在从头组装之前，可对序列读段进行标准化处理，以实现基因组范围内的均匀覆盖度。这一步骤至关重要，因为诸如 SPAdes 等从头组装软件在组装序列读段时会考虑读段覆盖度。使用 BBMap 软件包中的 BBNorm 将读段标准化至覆盖度为 50。³⁰在命令行中输入：
   - bbmap/bbnorm.sh target=50 in=9_1P_trimmed.fastq in2=9_2P_trimmed.fastq out=Sample9_FW_norm.fastq out2=Sample9_RE_norm.fastq
4. 使用 SPAdes 对标准化后的测序读长进行从头组装³¹使用所有不同的kmer值（21、33、55、77、99和127）进行组装时采用默认设置。在命令行中输入：
   - spades.py -k 21,33,55,77,99,127 -o Sample9 -1 Sample9.qc.f.fq -2 Sample9.qc.r.fq
5. 使用 minimap2 以及 Geneious、Bioedit 或 Ugene 等软件将质控后的测序读段比对至获得的共识序列，以校正比对结果。必须检查重叠群（contig）的起始端和末端。
   1. 使用 minimap2 将质控后的测序读段与获得的共识序列进行比对。
   2. 在 Geneious/Bioedit/UGene 中导入比对结果。
   3. 手动检查、校正和审编基因组的起始和终止区域。

5. 使用Illumina数据作为金标准，确定为补偿纳米孔测序错误特征所需的测序读长覆盖度

1. 选择比对到一个扩增子的序列读段，此处为扩增子26。随后使用 minimap2 将纳米孔测序读段比对至该扩增子。利用 Samtools 仅筛选比对到扩增子26的读段，并将 bam 文件转换为 fastq 格式。在命令行中输入：
   - minimap2 -ax map-ont -m 150 Amplicon26.fasta BC01_trimmed.fastq > BC01.bam
   - samtools view -b -F 4 BC01.bam > BC01_mapped.bam
   - samtools bam2fq BC01_mapped.bam | seqtk seq - -> BC01_mapped.fastq
2. 对例如200条序列读段的子集进行一千次随机抽样。例如，若将数值改为10，则将进行一千次每次随机抽取10条序列读段的子集。该脚本作为 补充文件1 提供。在命令行中输入：
   - python Random_selection.py
3. 将所有随机选取的序列读段比对至扩增子26。使用 KMA³² 进行序列读段比对，并立即生成一致性序列。采用纳米孔测序优化设置，即使用 -bcNano 标志。在命令行中输入：
   - kma index -i Amplicon26.fasta
   - for file in random_sample*; do
   - sampleID=${file%.fastq}
   - kma -i ${sampleID}.fastq -o ${sampleID} -t_db Amplicon26.fasta -mem_mode -mp 5 -mrs 0.0 -bcNano
   - done
4. 使用以下命令行指令检查生成的一致性序列：
   - cat *.fsa > All_genomes.fsa
   - minimap2 -ax map-ont Amplicon26.fasta All_genomes.fsa > All_genomes.bam
   - samtools sort All_genomes.bam > All_genomes_sorted.bam
   - samtools stats All_genomes_sorted.bam > stats.txt
   1. 错误率显示在 stats.txt 文件中标题为 error rate #mismatches / bases mapped 的部分。使用以下命令在屏幕上显示：
      - grep ^SN stats.txt | cut -f 2-
   2. 插入缺失（indel）的数量显示在标题为 #Indels per cycle 的部分。使用以下命令在屏幕上显示：
      - grep ^IC stats.txt | cut -f 2-

结果

最近，新型流动槽版本（R10）已发布，该版本改进了用于将电子电流信号转换为DNA序列的碱基识别软件（即所谓的“翻转-翻转”碱基识别器）。因此，我们对先前在R9.4流动槽和Illumina MiSeq测序仪上测序过的一只USUV阳性猫头鹰脑组织中的USUV进行了重新测序²¹。本文描述了通过与Illumina测序结果直接比较，确定可靠获得一致序列所需测序深度的方法。

通过使用新型流动池并结合翻转-翻转（flip-flop）碱基识别方法，我们证明，40x的测序深度可获得与Illumina测序相当的完全一致结果。30x测序深度的错误率为0.0002%，相当于每测序585,000个核苷酸出现一个错误；而20x测序深度则为每测序63,529个核苷酸出现一个错误。10x测序深度下，每测序3,312个核苷酸出现一个错误，这意味着在完整的USUV基因组中，平均有超过三个核苷酸被错误识别。当测序深度高于30x时，未观察到任何插入缺失（indel）。20x测序深度检测到一个indel位点，而10x测序深度则在29个位点上检测到indel。不同测序深度阈值下的错误率概览见表1。

覆盖度	第一次迭代错误数	第一次迭代错误率	插入缺失：	第二次迭代错误数	第二次迭代错误率	插入缺失：	第三次迭代错误数	第三次迭代错误率	插入缺失：
10×	100	0.0274%	4	116	0.0297%	18	110	0.0282%	7
20×	4	0.0010%	0	6	0.0015%	1	7	0.0018%	0
30x	2	0.0005%	0	0	0.0000%	0	0	0.0000%	0
40x	0	0.0000%	0	0	0.0000%	0	0	0.0000%	0
50×	0	0.0000%	0	0	0.0000%	0	0	0.0000%	0

表1：纳米孔测序错误率概述。 每次迭代代表一千个随机样本。

补充文件 1：随机选择。 请点击此处查看该文件（右键点击可下载）。

讨论

纳米孔测序技术不断发展，因此需要方法来监测其错误率。本文介绍了一种用于监测纳米孔测序仪错误率的工作流程。该方法在新型流动池发布后，或碱基识别软件有新版本发布时尤为有用。此外，对于希望建立并验证自身测序方案的用户而言，该方法同样具有实用价值。

不同的软件和比对工具可能会产生不同的结果³³。在本研究中，我们力求使用常见且可免费获取的软件包，这些软件具有清晰的文档说明。在某些情况下，可能会优先选择商业工具，这类工具通常具有更友好的用户界面，但需要付费使用。未来，如果测序技术或碱基识别软件发生重大更新，该方法可再次应用于同一样本。理想情况下，应在每次碱基识别软件或流通池更新后进行此操作，但鉴于当前技术发展的速度，也可仅在重大更新后进行。

测序错误率的降低使得可 multiplex（多重检测）的样本数量增加。因此，纳米孔测序正逐步接近取代传统的实时PCR用于诊断检测，流感病毒的诊断已实现这一替代。此外，错误率的降低也提高了该技术的可用性，例如用于稀有变异株的鉴定以及高通量、无偏倚的宏基因组测序。

该实验方案中的一个关键步骤是必须具备可靠且高度相近的参考序列。引物的设计基于当前对病毒多样性的认知，可能需要定期更新。在建立基于扩增子的测序方法时，另一个关键点是平衡引物浓度，以获得均匀的扩增子深度。这使得在一次测序运行中可对更多样本进行多重分析，从而显著降低检测成本。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Agencourt AMPure XP 磁珠	Beckman Coulter	A63881
dNTP 混合物	Qiagen	201900
FLO-MIN106 R10 流式芯片	Nanopore	R10 flowcell
KAPA Hyperplus 文库制备试剂盒	Roche	7962436001
文库加载磁珠试剂盒	Nanopore	EXP-LLB001
1D 连接测序试剂盒	Nanopore	SQK-LSK109
1D 原始条形码试剂盒 1-12	Nanopore	EXP-NBD103
1D 原始条形码试剂盒 13-24	Nanopore	EXP-NBD104
NEB 平末端/TA 连接酶主混合液	NEB	M0367S
NEB Next 快速连接模块	NEB	E6056
NEB Next Ultra II 末端修复 / dA 加尾模块	NEB	E7546S
Protoscript II 逆转录酶	NEB	M0368X
Q5 高保真聚合酶	NEB	M0491
Qubit dsDNA HS 检测试剂盒	Thermo Fisher	Q32851
随机引物	Promega	C1181
RNAsin 核糖核酸酶抑制剂	Promega	N2111