一种基于纳米孔测序数据的病原宏基因组分析方法与流程

1.本发明涉及基因组测序技术领域，具体涉及一种基于纳米孔测序数据的病原宏基因组分析方法。


背景技术：

2.临床样本在采集过程中不可避免会有定植菌，实验过程中也会引入试剂背景菌，这些细菌都可能会引入耐药基因，造成假阳性，给临床精准用药带来困扰。基于二代测序短读长数据进行耐药分析，很难区分哪些耐药基因是定植菌/试剂背景菌携带的，哪些是致病菌携带的，难以建立病原微生物和耐药基因的关系。虽有各种算法模型来分析病原微生物和耐药基因的归属，但总不如纳米孔测序长读长数据直接比对来的准确。
3.病原宏基因组检测产生大量的测序数据，既包含病原体的核酸信息，也包含大量的人源序列，此外还不可避免的混有背景(如环境、工程试剂菌)污染序列。现有的二代病原宏基因组测序数据分析方法易出现多重比对和错误比对，有假阳性风险，难以保证鉴定结果的准确性。
4.数据库的建设是保证鉴定结果准确性的一个重要方面，现阶段在世界范围内已建立多个基因组序列信息数据库，如genbank、sra、ena等，这些数据库包含了大量病原微生物基因组dna序列及注释信息，但这些数据库中的部分数据来源于用户提交的序列信息，因此需要对数据库进行梳理、确证、分类、提取和整合。此外，还可自建背景污染序列等数据库。数据库建设程度也会限制鉴定结果的准确性。单个分析流程则很难实现在保证鉴定结果准确性基础上，大幅加快分析速度的要求。多个分析流程并行运行，并使用不同的数据库进行基因序列比对，才能保证鉴定结果的准确性和及时性。
5.基于bwt转换算法的bwa算法主要针对读长＜300bp的短读长数据进行分析，不能比对分析纳米孔测序的长读长(》2kb)数据，而blast算法对大量数据进行比对时速度较慢。因此，探索纳米孔测序长读长数据的分析方法是当下利用好长读长数据至关重要的一环。
6.鉴于上述缺陷，本发明创作者经过长时间的研究和实践终于获得了本发明。


技术实现要素：

7.本发明的目的在于解决现有的序列比对软件用于纳米孔测序，会降低分析的时效性、准确性和灵敏性的问题，提供了一种基于纳米孔测序数据的病原宏基因组分析方法。
8.为了实现上述目的，本发明公开了一种基于纳米孔测序数据的病原宏基因组分析方法，包括以下步骤：
9.s1，读取下机数据，采用国际通用的测序数据标准fastq格式，支持gzip方法进行数据压缩以减少存储占用；
10.s2，根据数据类型，使用nanofilt软件和fastqc软件和对应参数进行低质量数据过滤；
11.s3，使用samtools软件提取通过数据质控的高质量测序数据，用于后续分析；
12.s4，第一轮快速比对，通过并行kraken2、centrifuge、pandora和minimap2四种快速比对软件和算法进行物种初步判别；
13.s5，基于步骤s3和步骤s4的结果，使用bcftools软件提取疑似病原的读段并按初步鉴定的物种进行拆分；
14.s6，将步骤s5拆分的疑似病原读段，与各物种参考序列进行blast比对，提取各读段最优比对，并按照比对长度大于读段长度80％、序列相似度高于90％，比对期望值小于1e-5，进行筛选；
15.s7，将步骤s6验证的病原读段按病原类型和基因组类型行各分类水平的物种统计分析，明确该份样本中存在的病原类型、基因组类型属水平和种水平的读段数量及构成比；
16.s8，将步骤s3中生成数据进行耐药基因与突变数据库快速比对，通过使用minimap2快速比对算法进行物种初步判别；
17.s9，基于步骤s3和s8的结果，使用samtools软件提取疑似病原体耐药基因的读段并按初步鉴定结果生成fasta序列，随后使用bcftools获得一致性序列；
18.s10，将步骤s9中获得的疑似病原体耐药基因读段，与各物种参考序列进行blast比对，提取各读段最优比对，并按照比对长度大于读段长度80％、序列相似度高于90％，比对期望值小于1e-5，进行筛选；
19.s11，将步骤s10中验证的病原体耐药基因读段结合步骤七生成的读段taxid号码进行关联，并按物种进行丰度统计分析，明确该份样本中存在的耐药菌、耐药基因基因及其点突变的类型，读段数量及构成比；
20.s12，根据步骤s7和步骤s11中获得的样本中病原体统计信息，及耐药菌、耐药基因信息，补充添加病原体注释信息、样本信息和患者信息等，生成检测报告。
21.所述步骤s1中使用cat或zcat命令合并来自于同一份样本的多个fastq数据文件或压缩的fastq.gz文件。
22.所述步骤s2中针对三代测序的fastq测序数据文件进行质控检测和过滤，使用nanofilt软件，筛选平均测序数据质量大于q7的测序读段。
23.所述步骤s2中针对二代数据使用fastqc软件，筛选平均测序数据质量大于q30的测序读段。
24.所述步骤s4中kraken2为基于物种特异性k-mer序列的细菌、真菌、病毒和寄生虫识别的算法，centrifuge为快速宏基因组序列分类方法，用于细菌、真菌、病毒和寄生虫识别，pandora为快速并对序列识别方法，minimap2为针对病毒等基因组变异较大病原体的补充比对，基于minimizer哈希表搜索、chaining算法和动态规划算法。
25.所述步骤s7中病原类型包括细菌、真菌、病毒和寄生虫，基因组类型包括dna和rna。
26.与现有技术比较本发明的有益效果在于：本发明针对纳米孔测序产生的长读长数据，整合了多个数据库，构建了多个序列比对分析的病原识别并行分析的算法，保证了病原体鉴定和耐药基因分析的准确性和时效性。首先在第一轮比对中，通过并行四种快速比对算法实现了针对典型数据(1gb纳米孔测序数据)的物种初步判别；在第二轮和第三轮验证分析中，使用精确比对算法blast分别将上一轮判别的各个物种的读段与该物种参考序列及本样本识别到的其他病原参考序列进行交叉比对验证，有效减少了精确比对的计算量。
通过minimap2与耐药基因/点突变数据库进行比对，并通过blast精确比对算法进行耐药基因/点突变的识别。通过构建该方法和系统，实现了典型纳米孔测序病原宏基因组数据在30min内的准确分析和报告生成，并有效兼容二代测序数据分析。
附图说明
27.图1为本发明流程图；
28.图2为本发明测试的计算环境和耗时统计；
29.图3为本发明与常见二代分析软件的耗时比对；
30.图4为实施例1的检测质控图；
31.图5为实施例2的检测质控图。
具体实施方式
32.以下结合附图，对本发明上述的和另外的技术特征和资源消耗作更详细的说明。
33.数据检测过程：
34.1、核酸提取
35.将样本转入50ml无菌离心管中，在台式低速离心机中进行离心，条件为3500rpm离心10min，离心结束后，移除上清移除，用无菌1ml塑料滴管吹打混匀离心得到的沉淀，转移到新的1.5ml离心管中，用试剂盒提取核酸。
36.2、用qubit 4.0測定核酸浓度，具体操作为
37.1)将缓冲液与核酸染液按照199:1的比例配制工作液；
38.2)将工作液与标准品按照190:10的比例配制标准液1和2；
39.3)将工作液与样品按照199:1的比例配制样品液；
40.4)用标准液1和2绘制标准曲线；；
41.5)在步骤4)的标准曲线下测定样品液中样品的浓度。
42.3、宏基因组测序文库的构建，具体为：
43.1)文库制备，在pcr管中，加入3μl模板dna、1μl fragmentation混合液、20μl无核酸酶水、1μl rlb(01-12a)和25μl 2x master mix，将pcr管放入pcr仪，程序如下：98℃3min，(98℃15s
‑‑‑
54℃20s
‑‑‑
75℃3min)25cycle
‑‑‑
75℃ 7min；
44.2)磁珠纯化：将pcr产物转入新的1.5ml离心管，加入40ul磁珠并混匀，室温吸附5min，瞬时离心后将离心管放置在磁力架上，室温静置10min，待溶液清亮后用移液器移除上清，75％乙醇进行漂洗并移除上清,室温晾晒10min，加入dnase/rnase-free h2o溶解，取1μl用qubit定量；
45.3)上机测序：将barcode与样本进行混合，离心后，吸取200ul混合样本，缓慢加入芯片中，盖上盖子，选择程序参数，设置文件名、样本名和样本类型，选择实时basecalling进行实时测序。
46.4、纳米孔病原宏基因组分析系统进行分析
47.如图1和图2所示，本发明提出的基于纳米孔测序数据的病原宏基因组分析方法以及固定计算资源下的时间消耗如下：
48.(1)读取下机数据，采用国际通用的测序数据标准fastq格式，支持gzip方法进行
数据压缩以减少存储占用；使用cat或zcat命令合并来自于同一份样本的多个fastq数据文件或压缩的fastq.gz文件；采用数据分批上传和分析的方式，该步骤耗时约2分钟；
49.(2)根据数据类型，使用不同数据质控软件和参数进行低质量数据过滤：
50.①
针对三代测序的fastq测序数据文件进行质控检测和过滤，使用nanofilt软件，筛选平均测序数据质量大于q7的测序读段；
51.②
针对二代数据使用fastqc软件，筛选平均测序数据质量大于q30的测序读段；
52.(3)提取通过数据质控的高质量测序数据，用于后续分析；
53.(4)第一轮快速比对，通过使用以下四种快速比对算法进行物种初步判别：
54.①
kraken2，基于物种特异性k-mer序列的细菌、真菌、病毒和寄生虫识别；
55.②
centrifuge，快速宏基因组序列分类方法，用于细菌、真菌、病毒和寄生虫识别；
56.③
pandora，快速并对序列识别方法；
57.④
minimap2，针对病毒等基因组变异较大病原体的补充比对，基于minimizer(一段序列内最小哈希值的种子)哈希表搜索、chaining算法和动态规划算法；
58.⑤
采用并行方式运行第一轮四种快速比对，该步骤耗时约10分钟；
59.(5)基于步骤(3)、(4)结果，使用bcftools提取疑似病原的读段并按初步鉴定的物种进行拆分，该提取目标序列步骤约耗时5分钟；
60.(6)将步骤(5)拆分的疑似病原读段，与各物种参考序列进行blast比对，提取各读段最优比对，并按照比对长度大于读段长度80％、序列相似度高于90％，比对期望值(e value)小于1e-5，进行筛选；
61.(7)将步骤(6)验证的病原读段按物种进行统计分析，明确该份样本中存在的病原类型(细菌、真菌、病毒、寄生虫)、基因组类型(dna、rna)、属水平和种水平的读段数量及构成比；步骤(6)和步骤(7)再次进行blast与目标序列的比对耗时约为10分钟；
62.(8)将步骤(3)生成数据进行耐药基因与突变数据库快速比对，通过使用minimap2快速比对算法进行物种初步判别，该步骤与步骤(4)的物种比对并行运行，耗时约10分钟；
63.(9)基于步骤(3)、(8)结果，提取疑似病原体耐药基因的读段并按初步鉴定结果生成fasta序列，以及运用bcftools获得一致性序列，该步骤与步骤(5)并行运行，耗时5分钟；
64.(10)将步骤(9)获得的疑似病原体耐药基因读段，与各物种参考序列进行blast比对，提取各读段最优比对，并按照比对长度大于读段长度80％、序列相似度高于90％，比对期望值(e value)小于1e-5，进行筛选；
65.(11)将步骤(10)验证的病原体耐药基因读段结合步骤七生成的读段taxid号码进行关联，并按物种进行统计分析，明确该份样本中存在的耐药菌、耐药基因基因及其点突变的类型，读段数量及构成比；步骤(9)和步骤(10)，与步骤(6)和(7)并行运行，耗时约为10分钟
66.(12)根据步骤(7)和步骤(11)获得的样本中病原体统计信息，及耐药菌、耐药基因信息，补充添加病原体注释信息、样本信息和患者信息等，生成检测报告，该统计和自动产出报告过程耗时约3分钟。
67.图3展示了对同一批数据，使用本发明和常用软件metaphlan3进行分析的耗时比较。在dell t7920工作站中，我们配置了512g内存，intel 5220r cpu 2块，2t的ssd存储10块，nvidia rtx a2000(12g)显卡一块，在该计算资源下，本发明对1g的纳米孔数据进行分
析耗时在25-30分钟之间，显著优于metaphlan3软件的45-50分钟。
68.实施例1
69.2021年7月23日采用这种基于纳米孔测序数据的病原宏基因组分析方法对从医院采集的肺部感染患者的肺泡灌洗液进行检测，检测结果如下：
70.1、检出病毒列表
71.(1)dna病毒列表
[0072][0073]
(2)rna病毒列表
[0074][0075]
(3)检出结核/非结核分枝杆菌列表
[0076][0077]
检测质控图如图4所示，由图4可知，此次检测测序数据质量合格，结果可信。
[0078]
实施例2
[0079]
2021年9月1日采用这种基于纳米孔测序数据的病原宏基因组分析方法对从医院采集的肺部感染患者的肺泡灌洗液进行检测，检测结果如下：
[0080]
1、检出细菌列表
[0081][0082]
[0083]
2、检出真菌列表
[0084][0085]
3、检出病毒列表
[0086]
(1)dna病毒列表
[0087][0088]
(2)rna病毒列表
[0089][0090]
4、检出结核/非结核分枝杆菌列表
[0091][0092]
5、检出非典型病原体列表
[0093][0094]
6、检出寄生虫列表
[0095][0096][0097]
7、人体口腔及呼吸道菌群列表
[0098][0099]
8、耐药基因检出列表
[0100]
[0101][0102]
检测质控图如图5所示，由图5可知，此次检测测序数据质量合格，结果可信。
[0103]
以上所述仅为本发明的较佳实施例，对本发明而言仅仅是说明性的，而非限制性的。本专业技术人员理解，在本发明权利要求所限定的精神和范围内可对其进行许多改变，修改，甚至等效，但都将落入本发明的保护范围内。

技术特征：
1.一种基于纳米孔测序数据的病原宏基因组分析方法，其特征在于，包括以下步骤：s1，读取下机数据，采用国际通用的测序数据标准fastq格式，支持gzip方法进行数据压缩以减少存储占用；s2，根据数据类型，使用nanofilt软件和fastqc软件和对应参数进行低质量数据过滤；s3，使用samtools软件提取通过数据质控的高质量测序数据，用于后续分析；s4，第一轮快速比对，通过并行kraken2、centrifuge、pandora和minimap2四种快速比对软件和算法进行物种初步判别；s5，基于步骤s3和步骤s4的结果，使用samtools软件提取疑似病原的读段并按初步鉴定的物种进行拆分；s6，将步骤s5拆分的疑似病原读段，与各物种参考序列进行blast比对，提取各读段最优比对，并按照比对长度大于读段长度80％、序列相似度高于90％，比对期望值小于1e-5，进行筛选；s7，将步骤s6验证的病原读段按病原类型和基因组类型行各分类水平的物种统计分析，明确该份样本中存在的病原类型、基因组类型属水平和种水平的读段数量及构成比；s8，将步骤s3中生成数据进行耐药基因与突变数据库快速比对，通过使用minimap2快速比对算法进行物种初步判别；s9，基于步骤s3和s8的结果，使用samtools软件提取疑似病原体耐药基因的读段并按初步鉴定结果生成fasta序列，随后使用bcftools获得一致性序列；s10，将步骤s9中获得的疑似病原体耐药基因读段，与各物种参考序列进行blast比对，提取各读段最优比对，并按照比对长度大于读段长度80％、序列相似度高于90％，比对期望值小于1e-5，进行筛选；s11，将步骤s10中验证的病原体耐药基因读段结合步骤七生成的读段taxid号码进行关联，并按物种进行丰度统计分析，明确该份样本中存在的耐药菌、耐药基因基因及其点突变的类型，读段数量及构成比；s12，根据步骤s7和步骤s11中获得的样本中病原体统计信息，及耐药菌、耐药基因信息，补充添加病原体注释信息、样本信息和患者信息等，生成检测报告。2.如权利要求1所述的一种基于纳米孔测序数据的病原宏基因组分析方法，其特征在于，所述步骤s1中使用cat或zcat命令合并来自于同一份样本的多个fastq数据文件或压缩的fastq.gz文件。3.如权利要求1所述的一种基于纳米孔测序数据的病原宏基因组分析方法，其特征在于，所述步骤s2中针对三代测序的fastq测序数据文件进行质控检测和过滤，使用nanofilt软件，筛选平均测序数据质量大于q7的测序读段。4.如权利要求1所述的一种基于纳米孔测序数据的病原宏基因组分析方法，其特征在于，所述步骤s2中针对二代数据使用fastqc软件，筛选平均测序数据质量大于q30的测序读段。5.如权利要求1所述的一种基于纳米孔测序数据的病原宏基因组分析方法，其特征在于，所述步骤s4中kraken2为基于物种特异性k-mer序列的细菌、真菌、病毒和寄生虫识别的算法，centrifuge为快速宏基因组序列分类方法，用于细菌、真菌、病毒和寄生虫识别，pandora为快速并对序列识别方法，minimap2为针对病毒等基因组变异较大病原体的补充
比对，基于minimizer哈希表搜索、chaining算法和动态规划算法。6.如权利要求1所述的一种基于纳米孔测序数据的病原宏基因组分析方法，其特征在于，所述步骤s7中病原类型包括细菌、真菌、病毒和寄生虫，基因组类型包括dna和rna。

技术总结
本发明涉及基因组测序技术领域，具体涉及一种基于纳米孔测序数据的病原宏基因组分析方法，首先在第一轮比对中，通过并行Kraken2、Centrifuge、Pandora和Minimap2四种快速比对软件和算法实现了针对典型数据(1Gb纳米孔测序数据)的物种初步判别；在第二轮和第三轮验证分析中，使用精确比对算法BLAST分别将上一轮判别的各个物种的读段与该物种参考序列及本样本识别到的其他病原参考序列进行交叉比对验证，有效减少了精确比对的计算量。通过Minimap2与耐药基因/点突变数据库进行比对，并通过blast精确比对算法进行耐药基因/点突变的识别。通过构建该方法和系统，实现了典型纳米孔测序病原宏基因组数据在30min内的准确分析和报告生成，并有效兼容二代测序数据分析。析。析。


技术研发人员：余乐 李寅虎 武志慧 高军涛 杜鹏程 刘树青
受保护的技术使用者：北京源生康泰基因科技有限公司
技术研发日：2023.06.16
技术公布日：2023/9/7