CN114023389B - 宏基因组数据的分析方法 - Google Patents
宏基因组数据的分析方法 Download PDFInfo
- Publication number
- CN114023389B CN114023389B CN202210002820.9A CN202210002820A CN114023389B CN 114023389 B CN114023389 B CN 114023389B CN 202210002820 A CN202210002820 A CN 202210002820A CN 114023389 B CN114023389 B CN 114023389B
- Authority
- CN
- China
- Prior art keywords
- sequence
- read
- data
- mer
- length
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 31
- 238000000034 method Methods 0.000 claims abstract description 59
- 241000894007 species Species 0.000 claims abstract description 57
- 238000000605 extraction Methods 0.000 claims abstract description 40
- 238000007405 data analysis Methods 0.000 claims abstract description 37
- 244000005700 microbiome Species 0.000 claims abstract description 20
- 238000007781 pre-processing Methods 0.000 claims abstract description 12
- 238000012163 sequencing technique Methods 0.000 claims description 51
- 238000007671 third-generation sequencing Methods 0.000 claims description 20
- 238000001914 filtration Methods 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 10
- 241000894006 Bacteria Species 0.000 claims description 6
- 238000007672 fourth generation sequencing Methods 0.000 claims description 6
- 238000011144 upstream manufacturing Methods 0.000 claims description 6
- 241000700605 Viruses Species 0.000 claims description 5
- 238000012935 Averaging Methods 0.000 claims description 4
- 238000007619 statistical method Methods 0.000 claims description 4
- 241000233866 Fungi Species 0.000 claims description 3
- 244000045947 parasite Species 0.000 claims description 3
- 230000035945 sensitivity Effects 0.000 abstract description 8
- 238000001514 detection method Methods 0.000 description 10
- 108020004414 DNA Proteins 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 9
- 238000011160 research Methods 0.000 description 7
- 230000000813 microbial effect Effects 0.000 description 5
- 238000012545 processing Methods 0.000 description 5
- 241000588724 Escherichia coli Species 0.000 description 3
- 206010028980 Neoplasm Diseases 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 230000007613 environmental effect Effects 0.000 description 3
- 206010009944 Colon cancer Diseases 0.000 description 2
- 238000001712 DNA sequencing Methods 0.000 description 2
- 102000016928 DNA-directed DNA polymerase Human genes 0.000 description 2
- 108010014303 DNA-directed DNA polymerase Proteins 0.000 description 2
- 241000194032 Enterococcus faecalis Species 0.000 description 2
- 206010064571 Gene mutation Diseases 0.000 description 2
- 241000186840 Lactobacillus fermentum Species 0.000 description 2
- 241000191967 Staphylococcus aureus Species 0.000 description 2
- 201000011510 cancer Diseases 0.000 description 2
- 238000013135 deep learning Methods 0.000 description 2
- 238000012217 deletion Methods 0.000 description 2
- 230000037430 deletion Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 229940032049 enterococcus faecalis Drugs 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 239000012634 fragment Substances 0.000 description 2
- 230000002068 genetic effect Effects 0.000 description 2
- 208000015181 infectious disease Diseases 0.000 description 2
- 238000003780 insertion Methods 0.000 description 2
- 230000037431 insertion Effects 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 229940012969 lactobacillus fermentum Drugs 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 108090000623 proteins and genes Proteins 0.000 description 2
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- 208000001333 Colorectal Neoplasms Diseases 0.000 description 1
- 102000053602 DNA Human genes 0.000 description 1
- 102000004190 Enzymes Human genes 0.000 description 1
- 108090000790 Enzymes Proteins 0.000 description 1
- 241000282412 Homo Species 0.000 description 1
- 241000186781 Listeria Species 0.000 description 1
- 241000186779 Listeria monocytogenes Species 0.000 description 1
- 102000048850 Neoplasm Genes Human genes 0.000 description 1
- 108700019961 Neoplasm Genes Proteins 0.000 description 1
- 102100030569 Nuclear receptor corepressor 2 Human genes 0.000 description 1
- 101710153660 Nuclear receptor corepressor 2 Proteins 0.000 description 1
- 238000012408 PCR amplification Methods 0.000 description 1
- 206010061902 Pancreatic neoplasm Diseases 0.000 description 1
- 241000607142 Salmonella Species 0.000 description 1
- 241001138501 Salmonella enterica Species 0.000 description 1
- 238000000540 analysis of variance Methods 0.000 description 1
- 239000002246 antineoplastic agent Substances 0.000 description 1
- 238000003556 assay Methods 0.000 description 1
- 230000001580 bacterial effect Effects 0.000 description 1
- 244000052616 bacterial pathogen Species 0.000 description 1
- 239000008280 blood Substances 0.000 description 1
- 210000004369 blood Anatomy 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000000973 chemotherapeutic effect Effects 0.000 description 1
- 229940044683 chemotherapy drug Drugs 0.000 description 1
- 208000029742 colonic neoplasm Diseases 0.000 description 1
- 238000010205 computational analysis Methods 0.000 description 1
- 239000012228 culture supernatant Substances 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000001962 electrophoresis Methods 0.000 description 1
- 230000005714 functional activity Effects 0.000 description 1
- 230000002538 fungal effect Effects 0.000 description 1
- 210000001035 gastrointestinal tract Anatomy 0.000 description 1
- 230000007614 genetic variation Effects 0.000 description 1
- 238000012165 high-throughput sequencing Methods 0.000 description 1
- 230000002458 infectious effect Effects 0.000 description 1
- 230000000968 intestinal effect Effects 0.000 description 1
- 208000015486 malignant pancreatic neoplasm Diseases 0.000 description 1
- 230000011987 methylation Effects 0.000 description 1
- 238000007069 methylation reaction Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 201000002528 pancreatic cancer Diseases 0.000 description 1
- 208000008443 pancreatic carcinoma Diseases 0.000 description 1
- 244000052769 pathogen Species 0.000 description 1
- 238000002203 pretreatment Methods 0.000 description 1
- 102000004169 proteins and genes Human genes 0.000 description 1
- 239000001397 quillaja saponaria molina bark Substances 0.000 description 1
- 230000036632 reaction speed Effects 0.000 description 1
- 229930182490 saponin Natural products 0.000 description 1
- 150000007949 saponins Chemical class 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 239000003053 toxin Substances 0.000 description 1
- 231100000765 toxin Toxicity 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
- G16B40/20—Supervised data analysis
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Medical Informatics (AREA)
- Health & Medical Sciences (AREA)
- Evolutionary Biology (AREA)
- Theoretical Computer Science (AREA)
- Biophysics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Health & Medical Sciences (AREA)
- Data Mining & Analysis (AREA)
- Biotechnology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Software Systems (AREA)
- Public Health (AREA)
- Artificial Intelligence (AREA)
- Evolutionary Computation (AREA)
- Epidemiology (AREA)
- Databases & Information Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Bioethics (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明提供了宏基因组数据的分析方法。本发明提供了一种宏基因组数据的分析方法，所述方法包括：1）对原始数据预处理，获得具有期望质量的数据集；2）对步骤1）的数据集中的每一条read序列进行N次的K‑mer滑动提取，获得N个K‑mer序列；3）将所有的read序列中同一次K‑mer滑动提取获得的K‑mer序列归为一个K‑mer序列子集，得到N个K‑mer序列子集；4）将步骤3）得到的每个K‑mer序列子集分别进行宏基因组物种分析，得到N个数据分析结果；5）将步骤4）得到的N个数据分析结果合并，分析并获得宏基因组中各种微生物的信息。本发明的方法具有超高的灵敏度，并且能有效控制假阳性结果。
Description
技术领域
本发明属于宏基因组分析技术领域，尤其涉及一种宏基因组数据的分析方法，更具体地，本发明涉及一种基于三代测序的宏基因组数据的分析方法。
背景技术
宏基因组（Metagenome），也称微生物环境基因组，即环境中全部微小生物遗传物质的总和，目前主要指环境样品中的细菌和真菌的基因组总和。宏基因组学（Metagenomics），是一种以环境样品中的微生物群体基因组为研究对象，以功能基因筛选和/或测序分析为研究手段，以微生物多样性、种群结构、进化关系、功能活性、相互协作关系及与环境之间的关系为研究目的的新的微生物研究方法，一般包括从环境样品中提取基因组DNA进行高通量测序分析。
近年来，科研人员针对宏基因组学的研究越来越广泛，尤其是针对人体宏基因组学的研究更是如日中天，例如肠道菌群的宏基因组学研究以及肿瘤宏基因组学的研究，研究者不仅把宏基因组序列进行了测序及分类，而且还分析研究了其与人类疾病之间重要的相关性。例如2017年，《科学》杂志发表的一项研究表明了微生物如何侵入大多数胰腺癌患者，并分解这些患者使用的主要化学治疗药物（Geller LT, Barzily-Rokni M, Danino T,et al. Potential role of intratumor bacteria in mediating tumor resistance tothe chemotherapeutic drug gemcitabine. Science. 2017;357(6356):1156-1160.DOI: 10.1126/science.aah5043）；2020年《自然》杂志发表研究，首次将人体内的微生物与驱动癌症发展的遗传变异建立了直接联系，证实结直肠癌基因突变可由肠道菌群释放毒素所导致（Pleguezuelos-Manzano C, Puschhof J, Rosendahl Huber A, et al.Mutational signature in colorectal cancer caused by genotoxic pks+ E. coli.Nature. 2020;580(7802):269-273. DOI: 10.1038/s41586-020-2080-8）；同年，另一研究团队在《自然》杂志上发表的关于癌症微生物组的最新成果（Poore GD, Kopylova E, ZhuQ, et al. Microbiome analyses of blood and tissues suggest cancer diagnosticapproach. Nature. 2020;579(7800):567-574. DOI: 10.1038/s41586-020-2095-1），通过宏基因组分析显示在大多数主要癌症类型的组织和血液中存在独特的微生物DNA特征。
第三代测序技术，又称三代测序技术（Third generation sequencing）或单分子实时DNA测序技术，是一种在DNA测序时，不需要经过PCR扩增即可实现对每一条DNA分子的单独测序的技术。目前第三代测序技术原理主要分为以Pacbio的SMRT技术为代表的单分子荧光测序以及以牛津纳米孔公司和齐碳科技公司的纳米孔电泳技术为代表的纳米孔测序。三代测序的主要的技术特点之一是实现了DNA聚合酶内在自身的反应速度，一秒可以测10个碱基，测序速度是化学法测序的2万倍；其二是实现了DNA聚合酶内在自身的延续性，一个反应就可以测非常长的序列；二代测序可以测到上百个碱基，但是三代测序就可以测几千个碱基。
宏基因组测序能够从临床样本中同时检测几乎所有已知病原体。目前大多数宏基因组研究都采用了二代（Illumina）测序平台，测序运行时间超过16小时，总体样本到报告的周期为48-72小时；相比之下，三代（纳米孔）测序在有效去除宿主DNA背景后（例如利用皂素Charalampous T, Kay GL, Richardson H, et al. Nanopore metagenomics enablesrapid clinical diagnosis of bacterial lower respiratory infection. NatBiotechnol. 2019;37(7):783-792. DOI: 10.1038/s41587-019-0156-5），实时计算分析情况下可以在50分钟左右检测到微生物，并且整个检测周期可控制在6小时内 （Gu W,Deng X, Lee M, et al. Rapid pathogen detection by metagenomic next-generationsequencing of infected body fluids. Nat Med. 2021;27(1):115-124. doi:10.1038/s41591-020-1105-z），目前可广泛应用于病原菌及新兴病毒的基因组检测上，具有非常好的临床应用前景，例如疫情的定点实时监控等。
然而，目前基于三代测序检测宏基因组的方法还存在很大的挑战和问题。在实验层面，如何高效稳定地去除宿主（host）基因组或者富集宏基因组的DNA片段是面临的主要的问题。实验方法不当会造成用于测序的DNA中宏基因组DNA的相对占比较低，使得流入后续数据分析的有效宏基因组数据过少，同时为了弥补该缺陷，会使用更大得测序数据量，从而造成测序成本消耗，性价比极低。在数据分析层面，不仅要考虑如何较高灵敏度地检出微生物种类且有效地控制假阳性结果，同时也要考虑如何解决三代数据的独有特征导致的数据利用率不高的现状，例如纳米孔测序序列中的随机indel问题。现阶段基于三代测序数据检测宏基因组的方法，例如Kraken2、Bracken以及CNN等一些深度学习的方法，均没有同时解决上述存在的问题。
因此，对现有的宏基因组数据分析方法进行进一步的改进，在提高检测灵敏度的同时，解决数据利用率不高的问题，具有非常重要的意义。
发明内容
因此，本发明的目的是针对现有技术的不足，提供一种宏基因组数据的分析方法，本发明提供的方法能够在数据分析层面上良好地解决了上述问题，不仅从数据特征上有效的规避掉随机***缺失（indel）导致的数据利用率不高的问题，同时具有超高的灵敏度，即使在超高的宿主（host）背景噪音下也可以精确地检出微生物的种类组成，并且有效的控制了假阳性的结果，使得分析结果更加准确高效；同时，本发明的方法为实验技术层面减轻了压力，而且不需要过多的测序数据量，降低了测序成本，充分发挥了三代测序数据的优势。本发明的方法也可以很好地兼容目前二代测序数据分析的方法，扩大了二代测序分析宏基因组数据的手段。
本发明的目的是通过以下技术方案实现的。
一方面，本发明提供了一种宏基因组数据的分析方法，所述方法包括以下步骤：
1）对原始数据预处理，获得具有期望质量的数据集，所述具有期望质量的数据集包括多条read序列；
2）对步骤1）的数据集中的所述多条read序列中的每一条read序列进行N次的K-mer滑动提取，针对每一条read序列，获得N个K-mer序列；
其中，在所述的N次滑动提取中增加约束条件，使滑动提取的K-mer序列的起始位置在目标序列read上随机均匀分布，滑动提取的K-mer序列之间有碱基重合，并且最终N次提取后的所有K-mer序列覆盖其目标序列read至少80%的区域；
其中，所述N为整数；
3）将所有的read序列中同一次K-mer滑动提取获得的K-mer序列归为一个K-mer序列子集，得到N个K-mer序列子集；
4）将步骤3）得到的每个K-mer序列子集分别进行宏基因组物种分析，最终得到N个数据分析结果；
5）将步骤4）得到的N个数据分析结果合并，分析并获得宏基因组中各种微生物的信息。
根据本发明的宏基因组数据的分析方法，其中，在步骤1）中，所述原始数据为原始三代测序数据或长读长测序数据。
其中，长读长是指三代测序例如纳米孔测序获得的测序数据，其原始测序数据的长度可以达到例如10kb、20kb、150kb、160kb、200kb，300kb，1Mb，2Mb。本领域技术人员已知的是，随着测序技术的发展，后续可能得到更长读长的测序数据。但是不难理解，本发明的方法适用于任何长度的测序数据，并且越长的数据，越能体现本发明的方法的优越性。
根据本发明的宏基因组数据的分析方法，其中，在步骤1）中，原始数据为经纳米孔测序获得的长读长数据。
根据本发明的宏基因组数据的分析方法，其中，在步骤1）中，对原始数据的预处理包括去除其中的接头序列及条形码（barcode）序列，过滤低质量read序列以及过短的read序列。
在任选的实施方案中，所述低质量的阈值可以为Q5，也可以为Q7，这是本领域技术人员熟知的阈值。其中，Q表示测序read的平均质量值，即测序read中的每一个碱基的准确率求和取平均后获得的值。该阈值可以根据实际情况进行调整，具体的调整参数详见http://maq.sourceforge.net/fastq.shtml，该处通过引用将其并入本文。
在任选的实施方案中，过短的read序列的序列长度阈值包括但不限于100bp；例如所述阈值可以为50bp、200bp、300bp等。本领域技术人员可以根据实际情况进行调整，该阈值调整参数详见http://maq.sourceforge.net/fastq.shtml，该处通过引用将其并入本文。
在任选的实施方案中，步骤1）预处理使用Porechop软件和/或NanoFilt软件。这是本领域技术人员常用的数据处理软件，其使用方法和技术参数也是本领域技术人员熟知的，该处不一一赘述。
本领域技术人员已知的是，K-mer是指将read分成包含K个碱基的字符串，例如ATCG序列为4-mer，ATCGATCG序列为8-mer。
根据本发明的宏基因组数据的分析方法，其中，步骤2）是通过包括以下步骤的方法实现的：
① 对每一条read序列，根据长度K进行N次滑动提取且不分割read序列，长度K为任意整数，若滑动提取后的长度小于预设长度则过滤掉该提取序列；
②增加约束条件，其中所述约束条件包含：
滑动提取的K-mer序列在目标read序列上的位置随机；
滑动提取的K-mer序列分布在目标read序列的上游、中游及下游位点；
滑动提取的K-mer序列的起始位置在目标read序列上均匀分布；以及
最终N次提取后的所有K-mer序列有碱基重合地覆盖目标read序列至少80%的区域；
③针对每一条read序列，获得N个K-mer序列。
在任选的实施方案中，K的取值范围为75bp≤K≤500bp，100bp≤K≤200bp或k=200bp。本领域技术人员已知的是，K可以根据数据的质量特点或者分许所需设置。无需任何理论的限制，K可以是任意的整数值，例如K可以大于等于20bp、大于等于50bp、大于等于75bp、大于等于100bp、大于等于150bp、大于等于200bp、大于等于300bp、大于等于400bp等。
根据本发明宏基因组数据的分析方法，其中，在步骤2）中，N的取值范围为10≤N≤50，20≤N≤25。无需任何理论的限制，N可以是任意的整数值，例如N可以为5、6、7、8、9、10、11、12、13、14、15、16、17、18、19、20、21、22、23、24、25、26、27、28、29、30、31、32、33、34、35、36、37、38、39、40、41、42、43、44、45、46、47、48、49、50或100、或200、或300等。
根据本发明的宏基因组数据的分析方法，其中，在步骤3）中，所述分析结果包括可以使用Centrifuge或Metaphlan等。这些软件都是本领域技术人员可以知晓并使用的软件，也是进行宏基因组分析的常规软件。
根据本发明所述的宏基因组数据的分析方法，其中，所述步骤5）是通过包括以下步骤的方法实现的：
将步骤4）得到的N个数据分析结果合并，得到检出物种与其相对丰度的结果矩阵，计算出每一个检出物种的相对丰度的平均值、方差及标准差；给定置信水平，例如但不限于90%、95%或99%等，并计算该物种在该置信水平下的相对丰度值的置信区间，以置信区间的范围为界，将N次分析结果中的满足置信区间的相对丰度值取平均值，将该计算值作为该物种检出的最终相对丰度，且所有检出物种的相对丰度值的总和为100%。
在一个具体的实施方案中，本发明的方法包括以下步骤：
1）对原始三代测序数据进行数据预处理，利用Porechop软件以及NanoFilt软件去除建库中加入的接头及条形码序列，过滤低质量以及过短的read序列，其中低质量的阈值为Q7，过短序列长度阈值为100bp，得到数据集；
2）将步骤1）数据集的每一条read序列，根据长度K进行N次滑动提取且不分割read，设置100bp≤K≤200bp，20≤N≤25，若滑动提取后的长度小于预设长度则过滤掉该提取序列，
在滑动提取中增加约束条件，包括：
滑动提取的K-mer序列在其目标read序列上的位置随机；
滑动提取的K-mer序列随机分布在目标read上的上游、中游及下游；
滑动提取的K-mer序列起始位点均匀分布在目标read上且无随机偏好；以及
最终N次提取后所有的K-mer序列有碱基重合地覆盖其目标序列read至少80%的区域；
针对每一条read序列，获得N个K-mer序列；
3）将所有的read序列中同一次K-mer滑动提取获得的K-mer序列归为一个K-mer序列子集，得到N个K-mer序列子集；
即对于每一条read序列，获得N个固定长度K-mer序列，不同read序列的同一次提取获得的K-mer序列属于同一个K-mer序列子集，从而将原始长读长的测序数据分解成了根据固定长度滑动提取的短读长的测序数据集合；
4）将步骤3）得到的最终长度为K的K-mer序列子集分别进行宏基因组物种分析，其中物种包含细菌、真菌、寄生虫和/或病毒，最终得到N个数据分析结果；
5）将步骤4）得到的N个数据分析结果合并，得到检出物种与其相对丰度的结果矩阵，计算出每一个检出物种的相对丰度的平均值、方差及标准差；给定置信水平为95%，并计算该物种在该置信水平下的相对丰度值的置信区间，以置信区间的范围为界，将N次分析结果中的满足置信区间的相对丰度值取平均值，将该计算值作为该物种检出的最终相对丰度，且所有检出物种的相对丰度值的总和为100%。
所述方法还包括将步骤4）得到的N个数据分析结果合并，并进行以下分析：
基于物种丰度矩阵对样本多样性进行统计分析，例如但不限于alpha diversity，Shannon指数及Simpson指数等；
基于物种丰度矩阵对样本组间差异显著物种进行统计分析，例如但不限于LDAEffect Size （LefSe）分析或方差分析等。
本发明还提供了一种用于宏基因组数据分析的装置，其中，所述装置包括：
数据预处理模块，用于对所述序列预处理，获得具有期望质量的数据集；
数据处理模块，用于将数据集处理为N个K-mer序列子集；
微生物鉴定模块，用于基于微生物参考序列数据库，从所述期望的数据集中分析得到微生物物种分类数据；
数据结果整合模块，对物种分类数据分析结果进行合并；及
报告模块，用于根据数据输出结果。
所述数据预处理模块包括：Porechop软件和/或NanoFilt软件；
所述数据处理模块包括使用本发明所述的方法，以将数据集处理为N个K-mer序列子集；
所述微生物鉴定模块包括Metaphlan软件、Centrifuge软件和/或宏基因组数据库。
本发明提出了一种基于三代测序数据检测宏基因组的方法及装置。具体地，本发明的发明人发现，对原始长读长的read序列进行滑动提取，转变成短读长的测序序列，能够更加充分地发挥三代测序数据长读长的优势。
现有的宏基因组分析方法虽然不用改变原始测序read序列的结构及数据特征，但是存在检测灵敏度不足出现假阴性的问题，本发明方法虽然改变了原始的测序数据结构，但本发明的方法创造性地对于每一条read序列进行了多次有碱基重合的抽样重复，很大程度上保留了原始测序数据包含的信息，同时规避了随机indel的数据特征，从而提高了比对率及数据利用率，再加上后面置信区间的“过滤”条件，整体对假阳性及假阴性都做了很好的控制。
因此，使用本发明的方法进行宏基因组数据分析，不仅能从数据特征上有效的规避掉随机***缺失（indel）导致的数据利用率不高的问题，同时具有超高的灵敏度，即使在超高的宿主（host）背景噪音下也可以精确地检出微生物的种类组成，并且有效地控制了假阳性的结果；为实验技术层面减轻了压力，而且不需要过多的测序数据量，降低了测序成本，充分发挥了三代测序数据的优势。本发明的方法也可以很好地兼容目前二代测序数据分析的方法，扩大了二代测序分析宏基因组数据的手段；同时作为目前三代测序数据分析宏基因组分析方法，使得分析结果更加准确高效。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对本发明实施例中所需要使用的附图作简单地介绍，显而易见地，下面所描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1示出为本发明的一个实施方案中宏基因组数据分析方法的流程框架图；
图2示出为本发明的一个实施方案中宏基因组数据分析装置的结构框图。
具体实施方式
下面将详细描述本发明的各个方面的特征和示例性实施例。在下面的详细描述中，提出了许多具体细节，以便提供对本发明的全面理解。但是，对于本领域技术人员来说很明显的是，本发明可以在不需要这些具体细节中的一些细节的情况下实施。下面对实施例的描述仅仅是为了通过示出本发明的示例来提供对本发明的更好的理解。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合附图对实施例进行详细描述。
然而，目前基于三代测序检测宏基因组的方法还存在很大的挑战和问题。在实验层面，如何高效稳定地去除宿主（host）基因组或者富集宏基因组的DNA片段是面临的主要的问题。实验方法不当会造成用于测序的DNA中宏基因组DNA的相对占比较低，使得流入后续数据分析的有效宏基因组数据过少，同时为了弥补该缺陷，会使用更大得测序数据量，从而造成测序成本消耗，性价比极低。在数据分析层面，不仅要考虑如何较高灵敏度地检出微生物种类且有效地控制假阳性结果，同时也要考虑如何解决三代数据的独有特征导致的数据利用率不高的现状，例如纳米孔测序序列中的随机indel问题。现阶段基于三代测序数据检测宏基因组的方法，例如Kraken2、Bracken以及CNN等一些深度学习的方法，均没有同时解决上述存在的问题。
因此，对现有的宏基因组数据分析方法进行进一步的改进，在提高检测灵敏度的同时，解决数据利用率不高的问题，具有非常重要的意义。
结合本发明的图1和图2，本发明提供了一种基于三代测序数据进行宏基因组数据分析的分析方法，所述方法包括以下步骤：
S1：对原始三代测序数据进行数据预处理，利用Porechop软件以及NanoFilt软件去除建库中加入的接头及条形码序列，过滤低质量以及过短的read序列，其中低质量的阈值为Q7，过短序列长度阈值为100bp，得到数据集；
S2：将步骤S1数据集的每一条read序列，根据长度K进行N次滑动提取且不分割read，设置100bp≤K≤200bp，20≤N≤25，若滑动提取后的长度小于预设长度则过滤掉该提取序列，
在滑动提取中增加约束条件，包括：
滑动提取的K-mer序列在其目标read序列上的位置随机；
滑动提取的K-mer序列随机分布在目标read上的上游、中游及下游；
滑动提取的K-mer序列起始位点均匀分布在目标read上且无随机偏好；以及
最终N次提取后所有的K-mer序列有碱基重合地覆盖其目标序列read至少80%的区域；
针对每一条read序列，获得N个K-mer序列；
S3：将所有的read序列中同一次K-mer滑动提取获得的K-mer序列归为一个K-mer序列子集，得到N个K-mer序列子集；
即对于每一条read序列，获得N个固定长度K-mer序列，不同read序列的同一次提取获得的K-mer序列属于同一个K-mer序列子集，从而将原始长读长的测序数据分解成了根据固定长度滑动提取的短读长的测序数据集合；
S4：将步骤S3得到的最终读长为K的K-mer序列子集分别进行宏基因组物种分析，其中物种包含细菌、真菌、寄生虫和/或病毒，最终得到N个数据分析结果；
S5：将步骤S4得到的N个数据分析结果合并，得到检出物种与其相对丰度的结果矩阵，计算出每一个检出物种的相对丰度的平均值、方差及标准差；给定置信水平为95%，并计算该物种在该置信水平下的相对丰度值的置信区间，最终将落在该置信区间内的相对丰度值的均值作为该物种检出的最终相对丰度值，且所有检出物种的相对丰度值的总和为100%。
本发明的发明人发现，对原始长读长的read序列进行滑动提取，转变成短读长的测序序列，能够更加充分地发挥三代测序数据长读长的优势。
进一步地，如图2所示，本发明的一个实施方案中，提供了一种用于宏基因组数据分析的装置，其中，所述装置包括：数据预处理模块101，使用例如Porechop软件和/或NanoFilt软件对所述序列预处理，获得具有期望质量的数据集；数据处理模块102，用于使用本发明的方法将数据集处理为N个K-mer序列子集；微生物鉴定模块103，用于基于微生物参考序列数据库和/或宏基因组数据库，使用例如Metaphlan软件、Centrifuge软件从所述期望的数据集中分析得到微生物物种分类数据；数据结果整合模块104，对物种分类数据分析结果进行合并；及报告模块105，用于根据数据输出结果。
实施例1 使用本发明的方法分析数据
1. 将含有0%、30%、50%以及90%的人源宿主比例的4组宏基因组标准品，每种标准品中均包含粪肠球菌（Enterococcus faecalis）、大肠杆菌（Escherichia coli）、发酵乳杆菌（Lactobacillus fermentum）、李斯特菌（Listeria monocytogenes）、沙门氏菌（Salmonella enterica）、金黄色葡萄球菌（Staphylococcus aureus），且各个菌在每组标准品中的期望丰度已知，通过实验文库制备，利用型号为QNome-9604的三代纳米孔测序仪进行测序，得到原始的长读长测序数据，所述数据的最大读长达到16Kb,平均长度7.8kb；利用Porechop软件以及NanoFilt软件去除实验建库过程中加入的接头及barcode序列，过滤低质量Q5以下以及长度小于100bp的read序列；
2. 将步骤1得到的序列的数据集中的每一条read序列，根据长度K=100bp进行20次（N=20）的滑动提取且不分割read，并且若滑动提取后的长度小于100bp则过滤掉该提取的短读长序列；
增加约束条件，经过20次滑动提取后得到的短读长序列在其目标read上的位置随机，并且其起始位点要较为均匀地分布在其目标read上的上游、中游及下游且滑动提取出的20条短读长序列集合覆盖了其目标序列read的85%；
3. 将每次提取后得到的K-mer序列归入一个K-mer序列子集，最终得到20个固定长度（K=100bp）的滑动提取后的K-mer序列子集，即将原始长读长的测序数据分解成了根据固定长度K=100bp滑动提取的短读长的测序数据集合；
4. 将步骤3得到的最终20个固定长度（K=100bp）的短读长序列数据集合分别利用metaphlan3进行宏基因组物种分析，可以得到20个数据分析结果；
5. 将20个数据分析结果合并，得到检出物种与其相对丰度的结果矩阵，分别计算出每一个检出物种的相对丰度的平均值、方差及标准差，给定95%置信水平，并计算该物种在95%置信水平下的相对丰度值的置信区间，最终将落在该置信区间内的相对丰度值的均值作为该物种检出的相对丰度值，所有检出物种的相对丰度值的总和为100%。
结果统计如表1所示，可见，本发明方法不论在纯菌环境的数据中，还是在超高宿主背景噪音的数据条件下，都可以非常灵敏的将各菌种检出，且没有任何假阳性结果，故本发明的方法特别适用于三代测序数据的宏基因组分析。
表1. 本发明方法检出各菌种以及其相对丰度的结果统计
注：由于metaphlan3在分析时不考虑人源宿主并且只计算比对到菌种的相对丰度，所以人源宿主并没有结果，用“-”表示。
另外，本文中术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。
应理解，在本发明实施例中，“与A相应的B”表示B与A相关联，根据A可以确定B。但还应理解，根据A确定B并不意味着仅仅根据A确定B，还可以根据A和/或其它信息确定B。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。
Claims (11)
1.一种宏基因组数据的分析方法，其特征在于，所述方法包括以下步骤：
1）对原始数据预处理，获得具有期望质量的数据集，所述具有期望质量的数据集包括多条read序列；
2）对步骤1）的数据集中的所述多条read序列中的每一条read序列进行N次的K-mer滑动提取，针对每一条read序列，获得N个K-mer序列；
其中，在所述的N次滑动提取中增加约束条件，使滑动提取的K-mer序列的起始位置在目标序列read上随机均匀分布，滑动提取的K-mer序列之间有碱基重合，并且最终N次提取后的所有K-mer序列覆盖其目标序列read至少80%的区域；
其中，所述N为整数；
3）将所有的read序列中同一次K-mer滑动提取获得的K-mer序列归为一个K-mer序列子集，得到N个K-mer序列子集；
4）将步骤3）得到的每个K-mer序列子集分别进行宏基因组物种分析，最终得到N个数据分析结果；
5）将步骤4）得到的N个数据分析结果合并，分析并获得宏基因组中各种微生物的信息；
步骤5）具体包括以下步骤：
将步骤4）得到的N个数据分析结果合并，得到检出物种与其相对丰度的结果矩阵，计算出每一个检出物种的相对丰度的平均值、方差及标准差；给定置信水平，并计算该物种在该置信水平下的相对丰度值的置信区间，以置信区间的范围为界，将N次分析结果中的满足置信区间的相对丰度值取平均值，将该计算值作为该物种检出的最终相对丰度，且所有检出物种的相对丰度值的总和为100%。
2.根据权利要求1所述的宏基因组数据的分析方法，其特征在于，在步骤1）中，所述原始数据为原始三代测序数据或长读长测序数据。
3.根据权利要求1所述的宏基因组数据的分析方法，其特征在于，在步骤1）中，所述原始数据为经纳米孔测序获得的长读长数据。
4.根据权利要求1所述的宏基因组数据的分析方法，其特征在于，在步骤1）中，对原始数据的预处理包括去除其中的接头序列及条形码序列，过滤质量低于Q7的read序列以及长度低于100bp的read序列。
5.根据权利要求1所述的宏基因组数据的分析方法，其特征在于，所述步骤2）包括以下步骤：
① 对每一条read序列，根据长度K进行N次滑动提取且不分割read序列，长度K为任意整数，若滑动提取后的长度小于预设长度则过滤掉该提取序列；
②增加约束条件，其中所述约束条件包含：
滑动提取的K-mer序列在目标read序列上的位置随机；
滑动提取的K-mer序列分布在目标read序列的上游、中游及下游位点；
滑动提取的K-mer序列的起始位置在目标read序列上均匀分布；以及
最终N次提取后的所有K-mer序列有碱基重合地覆盖目标read序列至少80%的区域；
③针对每一条read序列，获得N个K-mer序列。
6.根据权利要求5所述的宏基因组数据的分析方法，其特征在于，所述K的取值范围为75bp≤K≤500bp。
7.根据权利要求5所述的宏基因组数据的分析方法，其特征在于，所述K的取值范围为100bp≤K≤200bp。
8.根据权利要求1所述的宏基因组数据的分析方法，其特征在于，在步骤2）中，所述N的取值范围为10≤N≤50。
9.根据权利要求1所述的宏基因组数据的分析方法，其特征在于，在步骤2）中，所述N的取值范围为20≤N≤25。
10.一种宏基因组数据的分析方法，其特征在于，所述方法包括以下步骤：
1）对原始三代测序数据进行数据预处理，利用Porechop软件以及NanoFilt软件去除建库中加入的接头及条形码序列，过滤低质量以及过短的read序列，其中低质量的阈值为Q7，过短序列长度阈值为100bp，得到数据集；
2）将步骤1）数据集的每一条read序列，根据长度K进行N次滑动提取且不分割read，设置100bp≤K≤200bp，20≤N≤25，若滑动提取后的长度小于预设长度则过滤掉该提取序列；
在所述滑动提取中增加约束条件，包括：
滑动提取的K-mer序列在其目标read序列上的位置随机；
滑动提取的K-mer序列随机分布在目标read上的上游、中游及下游；
滑动提取的K-mer序列起始位点均匀分布在目标read上且无随机偏好；以及
最终N次提取后所有的K-mer序列有碱基重合地覆盖其目标序列read至少80%的区域；
针对每一条read序列，获得N个K-mer序列；
3）将所有的read序列中同一次K-mer滑动提取获得的K-mer序列归为一个K-mer序列子集，得到N个K-mer序列子集；
即对于每一条read序列，获得N个固定长度K-mer序列，不同read序列的同一次提取获得的K-mer序列属于同一个K-mer序列子集，从而将原始长读长的测序数据分解成了根据固定长度滑动提取的短读长的测序数据集合；
4）将步骤3）得到的最终长度为K的K-mer序列子集分别进行宏基因组物种分析，其中物种包含细菌、真菌、寄生虫和/或病毒，最终得到N个数据分析结果；
5）将步骤4）得到的N个数据分析结果合并，得到检出物种与其相对丰度的结果矩阵，计算出每一个检出物种的相对丰度的平均值、方差及标准差；给定置信水平为95%，并计算该物种在该置信水平下的相对丰度值的置信区间，以置信区间的范围为界，将N次分析结果中的满足置信区间的相对丰度值取平均值，将该计算值作为该物种检出的最终相对丰度，且所有检出物种的相对丰度值的总和为100%。
11.根据权利要求10所述的宏基因组数据的分析方法，其特征在于，还包括将步骤4）得到的N个数据分析结果合并，并进行以下分析：
基于物种丰度矩阵对样本多样性进行统计分析；
基于物种丰度矩阵对样本组间差异显著物种进行统计分析。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210002820.9A CN114023389B (zh) | 2022-01-05 | 2022-01-05 | 宏基因组数据的分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210002820.9A CN114023389B (zh) | 2022-01-05 | 2022-01-05 | 宏基因组数据的分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114023389A CN114023389A (zh) | 2022-02-08 |
| CN114023389B true CN114023389B (zh) | 2022-03-25 |
Family
ID=80069300
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210002820.9A Active CN114023389B (zh) | 2022-01-05 | 2022-01-05 | 宏基因组数据的分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114023389B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107034279A (zh) * | 2017-05-05 | 2017-08-11 | 中山大学 | 结核病微生物标志物在制备诊断结核病的试剂中的应用 |
| CN113862382A (zh) * | 2020-06-30 | 2021-12-31 | 北京大学人民医院 | 肠道菌群的生物标志物在制备诊断成人免疫性血小板减少症的产品中的应用 |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106228035B (zh) * | 2016-07-07 | 2019-03-01 | 清华大学 | 基于局部敏感哈希和非参数化贝叶斯方法的高效聚类方法 |
| JP2020500004A (ja) * | 2016-09-15 | 2020-01-09 | サン ジェノミクス インコーポレイテッド | 試料における１つまたは複数のタイプの微生物の多様な集団からの核酸分子の普遍的抽出方法 |
| CN108334750B (zh) * | 2018-04-19 | 2019-02-12 | 江苏先声医学诊断有限公司 | 一种宏基因组数据分析方法及*** |
| WO2019232357A1 (en) * | 2018-05-31 | 2019-12-05 | Arizona Board Of Regents On Behalf Of The University Of Arizona | Methods for comparative metagenomic analysis |
| US11830581B2 (en) * | 2019-03-07 | 2023-11-28 | International Business Machines Corporation | Methods of optimizing genome assembly parameters |
| CN111933218B (zh) * | 2020-07-01 | 2022-03-29 | 广州基迪奥生物科技有限公司 | 一种优化的宏基因组binning分析微生物群落的方法 |
| CN113160882B (zh) * | 2021-05-24 | 2022-11-15 | 成都博欣医学检验实验室有限公司 | 一种基于三代测序的病原微生物宏基因组检测方法 |
-
2022
- 2022-01-05 CN CN202210002820.9A patent/CN114023389B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107034279A (zh) * | 2017-05-05 | 2017-08-11 | 中山大学 | 结核病微生物标志物在制备诊断结核病的试剂中的应用 |
| CN113862382A (zh) * | 2020-06-30 | 2021-12-31 | 北京大学人民医院 | 肠道菌群的生物标志物在制备诊断成人免疫性血小板减少症的产品中的应用 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114023389A (zh) | 2022-02-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Press et al. | Hi-C deconvolution of a human gut microbiome yields high-quality draft genomes and reveals plasmid-genome interactions | |
| US20190194727A1 (en) | Multitag sequencing ecogenomics analysis | |
| Shendure et al. | Next-generation DNA sequencing | |
| CN111816258B (zh) | 人体菌群16S rDNA高通量测序物种精确鉴定的优化方法 | |
| US20210403991A1 (en) | Sequencing Process | |
| CN109923217A (zh) | 宏基因组样品中病原体的鉴定和抗生素表征 | |
| CN105420375B (zh) | 一种环境微生物基因组草图的构建方法 | |
| CN107077537A (zh) | 用短读测序数据检测重复扩增 | |
| Quijada et al. | High-throughput sequencing and food microbiology | |
| CA3054487A1 (en) | Systems and methods for metagenomic analysis | |
| CN103981256A (zh) | 一种沙门氏菌crispr分型方法 | |
| CN115662516A (zh) | 一种基于二代测序技术的高通量预测噬菌体宿主的分析方法 | |
| CN107577923A (zh) | 一种高度相似微生物的鉴定和分类方法 | |
| JP6588536B2 (ja) | 異なる種属の微生物間の種類と存在比を比較するための人工外来性参照分子 | |
| CN114023389B (zh) | 宏基因组数据的分析方法 | |
| Kim | Application of metagenomic techniques: understanding the unrevealed human microbiota and explaining the in clinical infectious diseases | |
| Vandamme et al. | Genomic taxonomy and biodiversity of the Burkholderia cepacia complex | |
| CN116705160A (zh) | 一种基于纳米孔测序数据的病原宏基因组分析方法 | |
| DS et al. | Microbial abundance estimation and metagenomic studies of holcomb creosote superfund site soil sample | |
| 豊間根耕地 | Studies on identification and evaluation of CRISPR diversity on human skin microbiome for development of a new personal identification method | |
| Jiksing et al. | Draft genome sequence data of two Salmonella bacteria from serogroup C type | |
| Perumal et al. | Downstream Analysis and Visualization-Knowledge Discovery–Alpha and Beta Diversity | |
| CN116705156A (zh) | 基于决策树算法寻找病毒基因组分类决定性位点的方法 | |
| CN116779036A (zh) | 一种基于多重pcr的靶向病原体纳米孔测序快速分析方法 | |
| CN112735530A (zh) | 一种基于菌群结构进行样品溯源的方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| REG | Reference to a national code |
Ref country code: HK Ref legal event code: DE Ref document number: 40061506 Country of ref document: HK |