CN110517726B - 一种基于高通量测序数据的微生物成分及浓度检测方法 - Google Patents

一种基于高通量测序数据的微生物成分及浓度检测方法 Download PDF

Info

Publication number
CN110517726B
CN110517726B CN201910637328.7A CN201910637328A CN110517726B CN 110517726 B CN110517726 B CN 110517726B CN 201910637328 A CN201910637328 A CN 201910637328A CN 110517726 B CN110517726 B CN 110517726B
Authority
CN
China
Prior art keywords
species
data
variable region
score
sequencing
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
Application number
CN201910637328.7A
Other languages
English (en)
Other versions
CN110517726A (zh
Inventor
袁细国
王双
刘畅
张军英
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201910637328.7A priority Critical patent/CN110517726B/zh
Publication of CN110517726A publication Critical patent/CN110517726A/zh
Application granted granted Critical
Publication of CN110517726B publication Critical patent/CN110517726B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biophysics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明属于基因序列数据检测技术领域，公开了一种基于高通量测序数据的微生物成分及浓度检测方法；量化碱基匹配或错配指示度、比对产生的gap与特定可变区指示度信息，建立基于NGS数据的读段‑物种库的隶属关系体系，通过设置隶属分数阈值过滤掉潜在的干扰读段，为准确鉴定与估计样本中物种成分与浓度奠定基础；提取物种的覆盖率、比对的gap得分、可变区得分三个特征，使用机器学习方法进行样本中物种成分的鉴定，提高物种成分鉴定的准确度。本发明基于高通量测序数据，研究病变样本中的微生物成分与浓度，脱离了以微生物培养为核心的传统临床检测手段，实现了速度快、准确率高的临床致病菌检测。

Description

一种基于高通量测序数据的微生物成分及浓度检测方法
技术领域
本发明属于基因序列数据检测技术领域，尤其涉及一种基于高通量测序数据的微生物成分及浓度检测方法。
背景技术
目前，最接近的现有技术：识别样本中物种成分常用的技术有基于序列间的相似性的计算方法、基于共享系统发育树对序列的划分方法和基于序列伪比对方法。基于序列间相似性方法的核心在于计算样本中测序读段间的相似性或测序读段与物种的参考序列的相似性。早期，基于相似性的方法使用BLAST工具(Basic Local Alignment SearchTool)将测序读段局部比对到参考序列，并计算比对的相似性分数。将测序读段分配到相似性分数最高的物种下，直至测序读段全部分配完毕，物种下有测序读段表明样本中含有该物种。在此期间，很多学者在BLAST的基础上对算法加以改进，但BLAST的计算量较大，所以很多以BLAST为雏形的算法很难适应于参考数据库规模扩大或数据测序深度增大的情形。这些算法很大程度上已经被usearch、uclust[Edgar RC.Search and clustering ordersof magnitude faster than BLAST.Bioinformatics.2010；26:2460–2461.][EdgarRC.UPARSE:highly accurate OTU sequences frommicrobial amplicon reads.NatMethods.2013；10:996–998.]和其他一些基于相似性聚类的算法[Al-Ghalith GA,Montassier E,Ward HN,Knights D.NINJA-OPS:Fast Accurate Marker Gene AlignmentUsing Concatenated Ribosomes.PLoS Comput Biol.2016；12:e1004658.Albanese D,Fontana P,De Filippo C,Cavalieri D,Donati C.MICCA:a complete and accuratesoftware for taxonomic profiling of metagenomic data.Sci Rep.2015；5:9743.MaheF,Rognes T,Quince C,de Vargas C,Dunthorn M.Swarm:robust and fast clusteringmethod for amplicon-based studies.PeerJ.2014；2:e593.Kopylova E,Noe L,TouzetH.SortMeRNA:fast and accurate filtering of ribosomal RNAs inmetatranscriptomic data.Bioinformatics.2012；28:3211–3217.]所取代，它们的速度与准确率也呈现大幅度提升并且大致彼此相当[Al-Ghalith GA,Montassier E,Ward HN,Knights D.NINJA-OPS:Fast Accurate Marker Gene Alignment Using ConcatenatedRibosomes.PLoS Comput Biol.2016；12:e1004658.Kopylova E,Navas-Molina JA,Mercier C,Xu ZZ,Mahe′F,He Y,etal.Open-Source Sequence Clustering MethodsImprove the State Of the Art.mSystems.2016；1]。
基于共享系统发育树的方法通过分析样本，将测序读段准确地放置在系统发育树上，实现物种成分的识别。近年来，使用最大似然估计[Berger SA,Krompass D,StamatakisA.Performance,accuracy,and Web server for evolutionary placement of shortsequence reads under maximum likelihood.Syst Biol.2011；60:291–302]、贝叶斯后验概率[Matsen FA,Kodner RB,Armbrust EV.pplacer:linear time maximum-likelihoodand Bayesian phylogenetic placement of sequences onto a fixed referencetree.BMC Bioinformatics.2010；11:538.]或邻域连接[Price MN,Dehal PS,ArkinAP.FastTree:computing large minimum evolution trees with profiles instead ofadistance matrix.Mol Biol Evol.2009；26:1641–1650.]的系统发育算法相继被开发出来，但它们的物种识别灵敏度较低[Bazinet AL,CummingsMP.A comparative evaluationof sequence classification programs.BMC Bioinformatics.2012；13:92.]而且计算量庞大。
基于序列伪比对的物种成分识别算法有Kallisto[Reppell M,Novembre J,Mchardy A C.Using pseudoalignment and base quality to accurately quantifymicrobial community composition[J].PLoS Computational Biology,2018,14(4)]，通过提取测序读段间共享的k-碱基序列[Rosen G,Garbarine E,Caseiro D,Polikar R,Sokhansanj B.Metagenome fragment classification using N-mer frequencyprofiles.Adv Bioinformatics.2008；2008:205969.McHardy AC,Martin HG,Tsirigos A,Hugenholtz P,Rigoutsos I.Accurate phylogenetic classification of variable-length DNA fragments.Nat Methods.2007；4:63–72.Wang Q,Garrity GM,Tiedje JM,ColeJR.Naive Bayesian classifier for rapid assignment ofrRNA sequences intothe new bacterial taxonomy.Appl Environ Microbiol.2007；73:5261–5267.]构建deBruijn图，计算测序读段来自于某个特定物种的可能性，进而判断样本中的物种成分。
对于物种浓度的估计，现有方法均通过统计特定物种下比对到的读段量占样本的比值作为样本中物种的相对浓度估计结果。
综上所述，现有技术存在的缺陷是：
(1)现有技术的检测方法的数据比对过程中存在大量的多比对读段，即一条读段会比对到多条物种序列的现象，该现象往往因测序错误、未知物种的干扰或生物的多样性而产生，对物种的成分识别与浓度估计带来影响。传统的三种研究未采取有效手段对多比对读段作出处理，导致物种成分鉴定结果的假阳性偏高。
(2)现有的检测方法仅使用物种下比对的读段量作为物种成分鉴定与浓度估计的唯一标准，缺乏对比对结果的空间分布形式、比对的错误类型及生物固有的遗传信息等众多方面的考虑，以致浓度估计不够准确。基于序列间的相似性的计算方法对物种识别的灵敏度较低。
(3)现有技术的计算复杂度高且存储量过大，算法的效率过低：基于序列伪比对的物种成分识别算法在构建de Bruijn图过程中会占据巨大的存储空间，基于共享系统发育树的方法在构建系统发育树时时间复杂度过高。
本发明所解决的传统方法中的缺陷：
(1)解决多比对读段问题：从比对形式而言，考虑匹配/错配(Match/MisMatch)、插入/缺失(统记为gap)情况；从测序错误而言，考虑序列比对的质量分数；从物种的生物特性而言，落在物种可变区的读段属于该物种的机率越高。该环节从比对成功的数据中量化碱基匹配或错配指示度、比对产生的gap与读段来自可变区的指示度等信息，建立测序读段到物种库的隶属关系体系，通过过滤掉比对到物种下的低隶属分数的读段来矫正测序读段的比对结果，实现较高准确度的物种浓度估计。
(2)打破仅以比对的读段量作为物种成分鉴定与浓度估计的唯一标准：从数量特征而言，考虑物种下比对的读段量；从空间特征而言，量化物种的位点覆盖与比对结果下产生gap的情况；从物种的生物特性而言，可变区是物种类型鉴定的重要信号，应从矫正结果中量化物种的可变区得分。基于以上各种角度的特征分析，构建物种存在性判别的分类器，从而实现样本中物种成分的识别。
(3)本发明的时间复杂度与空间复杂度较低：本技术的时间复杂度仅取决于物种可变区的提取、高通量测序数据的比对、读段-物种库隶属分数的计算、分类器的分类决策过程；本技术使用的分类器为基于核的支持向量机，空间复杂度仅与支持向量的个数有关。
解决上述技术问题的难度：
(1)物种可变区的提取：使用EMBOSS工具和可变区的引物文件提取物种的可变区序列，可变区的引物序列如下：
编号 前向引物序列 反向引物序列 可变区
1 AGYGGCGNACGGGTGAGTAA TGCTGCCTCCCGTAGGAGT V2
2 CCTACGGGAGGCAGCAG ATTACCGCGGCTGCTGG V3
3 AYTGGGYDTAAAGNG TACNVGGGTATCTAATCC V4
4 AGGATTAGATACCCT CCGTCAATTCCTTTGAGTTT V5
5 TCGAtGCAACGCGAAGAA ACATtTCACaACACGAGCTGACGA V6
6 GYAACGAGCGCAACCC GTAGCRCGTGTGTMGCCC V7
7 ATGGCTGTCGTCAGCT ACGGGCGGTGTGTAC V8
8 AGYGGCGNACGGGTGAGTAA ATTACCGCGGCTGCTGG V23
9 CCTACGGGRSGCAGCAG GGGGTATCTAATCCC V34
10 CCTACGGGAGGCAGCAG CCGTCAATTCCTTTGAGTTT V35
11 AYTGGGYDTAAAGNG CCGTCAATTYYTTTRAGTTT V45
12 GGATTAGATACCCTGGTAGTC ACAGCCATGCAGCACCT V56
13 CAACGCGAAGAACCTTACC GTAGCRCGTGTGTMGCCC V67
14 AGGTGCTGCATGGCTGT GACGGGCGGTGWGTRCA V78
一个可变区提取命令如下：
fuzznuc-sequence fa_file.fa-pattern‘forward_primer’-outfileforward.fuzznuc
fuzznuc-sequence fa_file.fa-pattern‘reverse_primer’-outfilereverse.fuzznuc
(2)基于核的支持向量机分类器的参数设定：
svm.SVC(kernel＝’rbf’,C＝1.0,random_state＝0,gamma＝0.2)
解决上述技术问题的意义：
(1)理论意义：本发明脱离传统的基于同源性的鉴定技术，提出衡量各种比对形式与判断物种存在性的全新方法，为解决现存算法中存在的类型鉴定的准确性与浓度检测灵敏度低的问题提供新思路，对微生物组学的基础研究有重要的理论意义。
(2)生物意义：本发明可揭示微生物群落中物种组成与浓度大小，为微生物群落的多样性分析、物种在环境中的依赖关系、种群数量与结构变化的规律、人类微生物组等众多方面研究提供新视角。
(3)应用价值：本技术以精确识别病变样本中致病菌与浓度为出发点，可实现临床治疗中快速、准确的致病菌诊断，避免抗生素滥用带来的病原体扩展和流行，从而使得对感染疾病的有针对性、有效和低毒副作用的精准药物治疗成为可能。
发明内容
针对现有技术存在的问题，本发明提供了一种基于高通量测序数据的微生物成分及浓度检测方法。
本发明是这样实现的，一种基于高通量测序数据的微生物成分及浓度检测方法，所述基于高通量测序数据的微生物成分及浓度检测方法量化碱基匹配或错配指示度、比对产生的gap与特定可变区指示度信息；建立基于NGS数据的读段-物种库的隶属关系体系；通过设置隶属分数阈值过滤掉潜在的干扰读段，为准确鉴定与估计样本中物种成分与浓度奠定基础；提取物种的覆盖率、比对的gap得分、可变区得分三个特征，使用机器学习方法进行样本中物种成分的鉴定。
进一步，所述基于高通量测序数据的微生物成分及浓度检测方法具体包括：
第一步，病原菌数据库选取与NGS数据预处理，对于获取到的NGS测序样本，采用FastQC工具分析数据的质量情况；基于测序读段数据及参考序列，采用BWA进行比对，保留比对成功的读段序列；
第二步，构建读段-物种库的隶属关系体系，从比对成功的数据中量化碱基匹配或错配指示度、比对产生的gap与读段来自可变区的指示度等信息，建立测序读段到物种库的隶属关系体系，通过过滤掉比对到物种下的低隶属分数的读段来矫正测序读段的比对结果；
第三步，过滤后的数据，提取物种下比对结果的位点覆盖率coverage、物种比对的gap分数记作gapscore、物种的可变区得分记作HVRScore，用于构建物种成分判别的分类器。
进一步，所述第二步的物种的参考数据库集合为F＝{f1,f2,…,fH}，测序读段集合为R＝{r1,r2,…,rN}；rj的长度为Lj，rj的碱基序列为rj＝(rj[1],rj[2],…,rj[Lj])，rj比对到物种fi后，比对形式记作fi＝(fi[1],fi[2],…,fi[Lj])；由碱基的测序质量值与错误率的关系Q＝-10*lg(qj[i])，得rj的碱基测序错误率(qj[1],qj[2],…,qj[Lj])；
rj比对到fi产生的gap集合为X＝(x1,x2,…,xM)，其中，xi为产生gap的宽度，M为gap的数量；若读段中碱基间的测序错误是独立的，则从碱基的匹配形式与测序错误角度而言，rj来源于fk的指示度衡量如下：
Figure GDA0004231777670000061
使用EMBOSS工具获取物种序列的可变区范围，rj比对到fi的可变区指示度量化为wij：
Figure GDA0004231777670000062
综合序列的比对形式、测序错误与可变区的指示度，将rj比对到fi的隶属分数定义如下：
Figure GDA0004231777670000063
将所有的测序数据比对到物种库后就可构建出读段-物种库的隶属关系体系；score(rj,fi)表示读段rj隶属于fi的程度。
进一步，所述第三步的成分鉴定与浓度估计具体包括：
a)从物种的位点覆盖率而言，位点覆盖率越高则该物种存在的可能性越大，定义如下：
Figure GDA0004231777670000071
b)从空间特征而言，若物种比对后得到的gap集合为G＝{g1,g2,…,gn}，gi表示第i个gap的长度，则物种比对结果下产生gap的得分如下：
Figure GDA0004231777670000072
c)从物种的生物特性而言，16S共包含9个可变区HVR，可变区下比对到的读段越多，则该物种存在的可能性越大；HVRi表示物种的第i个可变区，令xi为物种比对结果下HVRi产生的gap长度，令yi为物种HVRi的长度，从矫正结果中量化物种的可变区得分如下：
Figure GDA0004231777670000073
通过构建800×的仿真数据，计算物种库中每个物种的特征向量coverage,gapscore,HVRScore，用作SVM分类器的训练集，实现对病变样本中物种成分的鉴定。
本发明的另一目的在于提供一种应用所述基于高通量测序数据的微生物成分及浓度检测方法的基因序列数据检测系统。
本发明的另一目的在于提供一种应用所述基于高通量测序数据的微生物成分及浓度检测方法的信息数据处理终端。
综上所述，本发明的优点及积极效果为：本发明量化碱基匹配或错配指示度、比对产生的gap与特定可变区指示度信息，建立基于NGS数据的读段-物种库的隶属关系体系，通过设置隶属分数阈值过滤掉潜在的干扰读段，为准确鉴定与估计样本中物种成分与浓度奠定基础；提取物种的覆盖率、比对的gap得分、可变区得分三个特征，使用机器学习方法进行样本中物种成分的鉴定，提高物种成分鉴定的准确度。
本发明能够解决现有技术对病原菌成分检测灵敏度低的问题；能够解决现有技术对病原菌浓度估计偏差过大的问题；能够解决临床检测流程中存在的速度慢的问题。
本发明基于高通量测序数据，研究病变样本中的微生物成分与浓度，脱离了以微生物培养为核心的传统临床检测手段，实现了速度快、准确率高的临床致病菌检测；本发明构建读段-物种库的隶属关系体系，为多比对读段的重分配、传统比对结果的矫正提供新思路，这是目前众多方法中未曾考虑的；本发明构建NGS样本的预处理→构建读段-物种库隶属关系体系→物种成分鉴定与浓度估计的流程，不以物种下比对的读段量作为物种成分鉴定与浓度估计的唯一标准，脱离了传统的基于同源性比对的思想，对比对结果的空间分布形式、比对的错误类型及生物固有的遗传信息等众多方面的考虑，提高了物种成分鉴定与浓度估计的准确率。
附图说明
图1是本发明实施例提供的基于高通量测序数据的微生物成分及浓度检测方法流程图。
图2是本发明实施例提供的基于高通量测序数据的微生物成分及浓度检测方法实现流程图。
图3是本发明实施例提供的读段的多比对形式示意图。
图4是本发明实施例提供的对仿真数据中物种类型的类型与浓度检测结果示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
针对现有技术的检测方法的数据比对过程中存在大量的多比对读段，即一条读段会比对到多条物种序列的现象，该现象往往因测序错误、未知物种的干扰或生物的多样性而产生，对物种的成分识别与浓度估计带来影响。传统方法未采取有效手段对多比对读段作出处理，导致物种成分鉴定结果的假阳性偏高；现有的检测方法仅使用物种下比对的读段量作为物种成分鉴定与浓度估计的唯一标准，缺乏对比对结果的空间分布形式、比对的错误类型及生物固有的遗传信息等众多方面的考虑，以致浓度估计不够准确。本发明从物种参考序列的层面而言，提取物种的16S rDNA下读段的分布形式、位点覆盖率与可变区得分特征，构建判别物种存在性的分类器，实现样本中物种成分的精准鉴定。从测序读段的层面而言，在比对成功的数据中量化碱基匹配或错配指示度、比对产生的gap与特定可变区指示度等信息，建立测序读段到物种库的隶属关系体系，通过矫正测序读段的比对结果，实现较高准确度的物种浓度估计。
下面结合附图对本发明的应用原理作详细的描述。
如图1所示，本发明实施例提供的基于高通量测序数据的微生物成分及浓度检测方法包括以下步骤：
S101：病原菌数据库选取与NGS数据预处理，对于获取到的NGS测序样本，采用FastQC工具分析数据的质量情况；基于测序读段数据及参考序列，采用BWA进行比对，保留比对成功的读段序列；
S102：构建读段-物种库的隶属关系体系，从比对成功的数据中量化碱基匹配或错配指示度、比对产生的gap与读段来自可变区的指示度等信息，建立测序读段到物种库的隶属关系体系，通过过滤掉比对到物种下的低隶属分数的读段来矫正测序读段的比对结果；
S103：过滤后的数据，提取物种下比对结果的位点覆盖率coverage、物种比对的gap分数(记作gapscore)、物种的可变区得分(记作HVRScore)，用于构建物种成分判别的分类器。
下面结合附图对本发明的应用原理作进一步的描述。
本发明实施例提供的基于高通量测序数据的微生物成分及浓度检测方法以病原菌的16S rDNA序列为鉴定对象，从高通量测序数据出发，对病变样本中微生物成分与浓度进行准确地鉴定与估计。该方法涉及以下概念与技术：
16S rDNA：病变样本的病原体主要源于细菌，细菌rRNA按沉降系数分为3种，分别为5S、16S和23S rRNA。16S rDNA(简称作16S)是细菌染色体上编码16S rRNA相对应的DNA序列，存在于所有细菌染色体基因中。
NGS技术与高通量测序数据：NGS技术可一次性对几十万到几百万条DNA分子进行序列测定，其产生的DNA序列读段称为高通量测序数据。
序列比对技术：将DNA测序读段比对到物种的参考基因组上，可使用BWA工具实现序列比对。
16S rDNA的可变区提取：细菌的16S rDNA包含10个高度保守区和9个可变区。可变区序列因细菌不同而异，保守区序列基本保守，可通过EMBOSS工具实现细菌可变区的提取。
如图2所示，本发明实施例提供的基于高通量测序数据的微生物成分及浓度检测方法具体包括以下步骤：
(1)病原菌数据库选取与NGS数据预处理
本发明采用的病原菌参考数据库涵盖257个临床致病菌物种，除了临床常见的葡萄球菌、链球菌、绿脓杆菌、变形杆菌以及沙门氏菌等外，还包含有各种较为少见的病原菌，分布于9个门、18个纲、35个目、52个科及86个属中。对于获取到的NGS测序样本，采用FastQC工具分析数据的质量情况，若数据中含有测序接头与低质量序列，使用Trimmomatic工具去除测序接头和低质量序列。基于测序读段数据及参考序列，采用BWA进行比对，保留比对成功的读段序列，以便后续的物种成分与浓度分析。
(2)构建读段-物种库的隶属关系体系
对(1)得到的数据比对结果中往往含有大量的多比对读段，通常因测序错误、未知物种干扰或生物多样性的现象而造成，而这些现象直接反映在序列的比对形式与质量分数上。从比对形式而言，考虑匹配/错配(Match/MisMatch)、插入/缺失(统记为gap)情况；从测序错误而言，考虑序列比对的质量分数；从物种的生物特性而言，落在物种可变区的读段属于该物种的机率越高。该发明从比对成功的数据中量化碱基匹配或错配指示度、比对产生的gap与读段来自可变区的指示度等信息，建立测序读段到物种库的隶属关系体系，通过过滤掉比对到物种下的低隶属分数的读段来矫正测序读段的比对结果，实现较高准确度的物种浓度估计。
物种的参考数据库集合为F＝{f1,f2,…,fH}，测序读段集合为R＝{r1,r2,…,rN}。令rj的长度为Lj，rj的碱基序列为rj＝(rj[1],rj[2],…,rj[Lj])，rj比对到物种fi后，比对形式记作fi＝(fi[1],fi[2],…,fi[Lj])。由于测序仪产生的测序错误是无法避免的，而碱基质量值能够定量描述碱基测序的准确性，即可表征测序的可信性。由碱基的测序质量值与错误率的关系Q＝-10*lg(qj[i])，可得rj的碱基测序错误率(qj[1],qj[2],…,qj[Lj])。
表1
Figure GDA0004231777670000111
测序读段比对到物种的参考序列的结果如表1所示，碱基的匹配与错配、碱基的插入与缺失(统称为gap)，可直观地指示读段rj来自于物种fi的的可能性。比对结果中碱基的错配数量越少，产生gap的个数与gap的宽度越小，则该读段比对到该物种的可信程度越高。令rj比对到fi产生的gap集合为X＝(x1,x2,…,xM)，其中，xi为产生gap的宽度，M为gap的数量。若读段中碱基间的测序错误是独立的，则从碱基的匹配形式与测序错误角度而言，rj来源于fk的指示度衡量如下：
Figure GDA0004231777670000112
由于微生物物种间的相似性极高，所以区分物种类型的最重要的信号就在于物种的可变区。对于一个多比对读段，若该读段比对到物种A的可变区内，同时比对到其他物种B、C的保守区，则基于物种的生物特性，认为该读段更有可能来源于物种A。读段的多比对形式如图3所示，本发明使用EMBOSS工具获取物种序列的可变区范围，rj比对到fi的可变区指示度量化为wij：
Figure GDA0004231777670000121
综合序列的比对形式、测序错误与可变区的指示度，将rj比对到fi的隶属分数定义如下：
Figure GDA0004231777670000122
将所有的测序数据比对到物种库后就可构建出读段-物种库的隶属关系体系，该体系可用于比对结果的矫正。score(rj,fi)表示读段rj隶属于fi的程度，该值越大则读段来自该物种的可能性越高，也就是比对的正确率越高。病变样本中可能含有未知的微生物或人类基因的干扰，来自这些干扰物种的读段虽然会错误地比对到物种库的某些物种上，但是它们往往存在很多的碱基错配、大量的gap，并且以很低的概率落在可变区内，以致读段到该物种的隶属分数很低。本发明通过大量仿真实验探索，将隶属分数阈值设置为0.44，过滤掉错误比对的读段，排除多比对读段的错误比对部分，实现比对结果的矫正，这对物种的成分鉴定极为重要。
(3)成分鉴定与浓度估计
针对(2)中过滤后的数据，提取物种下比对结果的位点覆盖率coverage、物种比对的gap分数(记作gapscore)、物种的可变区得分(记作HVRScore)，用于构建物种成分判别的分类器。
a)从物种的位点覆盖率而言，位点覆盖率越高则该物种存在的可能性越大，其定义如下：
Figure GDA0004231777670000123
b)从空间特征而言，若物种比对后得到的gap集合为G＝{g1,g2,…,gn}，gi表示第i个gap的长度，则物种比对结果下产生gap的得分如下：
Figure GDA0004231777670000124
c)从物种的生物特性而言，16S共包含9个可变区(HighVariable Region,HVR)，可变区下比对到的读段越多，则该物种存在的可能性越大。令HVRi表示物种的第i个可变区，令xi为物种比对结果下HVRi产生的gap长度，令yi为物种HVRi的长度，从矫正结果中量化物种的可变区得分如下：
Figure GDA0004231777670000131
通过构建800×的仿真数据，计算物种库中每个物种的特征向量(coverage,gapscore,HVRScore)，用作SVM分类器的训练集，可实现对病变样本中物种成分的鉴定。
由于在隶属分数过滤过程中，可能会错误地过滤掉某些多比对读段，这对物种的浓度估计会造成干扰，所以需要回收过滤掉的多比对读段，并将这些读段重新比对到已鉴定出的物种上，统计每个物种下比对到的读段量，计算它们在测序数据中的相对浓度，作为物种浓度的估计结果。
下面结合实验对本发明的技术效果作详细的描述。
本发明使用仿真工具ART生成800×，Insertion/deletion均为0.1的11组物种种类、浓度、干扰程度均不同仿真数据，数据内容如表1：
表1
Figure GDA0004231777670000132
[软件来源ARThttps://www.niehs.nih.gov/research/resources/software/biostatistics/art/]
本发明命名为PGMicroD工具，与现存的五种工具：Bwa、Harp[Kessner D,Turner TL,Novembre J.Maximum Likelihood Estimation of Frequencies of Known Haplotypesfrom Pooled Sequence Data[J].MolecularBiology and Evolution,2013,30(5):1145-1158.]、Kallisto[Bray N L,Pimentel H,Melsted,Páll,et al.Near-optimalprobabilistic RNA-seq quantification[J].Nature Biotechnology,2016.]、Karp[Reppell M,Novembre J,Mchardy A C.Using pseudoalignment and base quality toaccurately quantify microbial community composition[J].PLoS ComputationalBiology,2018,14(4).]、Mothur[Zhang Y M,Tian C F,Sui X H,et al.Robust MarkersReflecting Phylogeny and Taxonomy ofRhizobia[J].PLOS ONE,2012,7.]
对仿真数据中物种类型的类型与浓度检测结果如图4：
Figure GDA0004231777670000141
其中，ti是物种的估计浓度，τi是物种的真实浓度；
MaxRE＝max|(ti-τi)|，其中，ti是物种的估计浓度，τi是物种的真实浓度。PGMicroD与其他五种工具对来自血液、尿液、脑脊液的25组临床样本的检测结果；在仿真数据和真实数据上的作用效果，PGMicroD均表现出较佳状态。
本实验处列举六种工具在第七组仿真数据上的物种类型与浓度检测结果如下所示，PGMicroD的结果明显优于其他五种工具。
表2真实物种类型编号与浓度
物种编号 浓度
64 0.035333
65 0.025563
66 0.007271
67 0.007195
68 0.00818
69 0.003295
70 0.001161
71 0
72 0
73 0
74 0
75 0
76 0
77 0
78 0
79 0
80 0
61 0.113271
62 0.091458
63 0.094531
表3 Harp检测的物种类型编号与浓度
物种编号 浓度
61 0.168566
62 0.135826
63 0.136898
64 0.074625
65 0.064214
66 0.010661
67 0.010383
68 0.012053
69 0.005206
70 0.003389
71 0.086029
72 0.051235
73 0.010019
74 0.054908
75 0.000238
76 0.00087
77 0.078043
78 0.096719
79 0.000118
80 0
表4 Karp检测的物种类型编号与浓度
物种编号 浓度
66 0.002556
65 0.015976
72 0.012307
70 0.000795
77 0.018754
68 0.002934
62 0.032872
79 7.78E-05
69 0.00129
76 0.000244
75 3.16E-05
61 0.041369
78 0.022849
63 0.033881
64 0.017379
71 0.02214
67 0.002487
74 0.013867
73 0.002438
表5 Kallisto检测的物种类型编号与浓度
物种编号 浓度
61 0.135069
62 0.091999
63 0.09434
64 0.046987
65 0.045036
66 0.008123
67 0.008035
68 0.008252
69 0.004709
70 0.010488
71 0.073654
72 0.042915
73 0.028636
74 0.06399
75 0.010298
76 0.010079
77 0.049311
78 0.084087
79 0.006296
表6 Bwa检测的物种类型编号与浓度
物种编号 浓度
64 0.048531
65 0.051353
66 0.011405
67 0.008445
68 0.010976
69 0.008174
70 0.015293
71 0.05904
72 0.046884
73 0.047995
74 0.060214
75 0.013987
76 0.023789
77 0.052142
78 0.093553
79 0.019863
80 0.010427
61 0.138985
表7 Mothur检测的物种类型编号与浓度
Figure GDA0004231777670000161
Figure GDA0004231777670000171
表8 PGMicroD检测的物种类型编号与浓度
物种编号 浓度
71 0.061754
61 0.192345
67 0.013173
77 0.053082
62 0.107849
63 0.109251
65 0.05252
66 0.011885
70 0.018948
68 0.011241
64 0.049617
69 0.012523
72 0.052129
73 0
74 0
75 0
76 0
78 0
79 0
80 0
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于高通量测序数据的微生物成分及浓度检测方法，其特征在于，所述基于高通量测序数据的微生物成分及浓度检测方法量化碱基匹配或错配指示度、比对产生的gap与读段来自可变区指示度信息；建立基于NGS数据的读段-物种库的隶属关系体系；通过设置隶属分数阈值过滤掉潜在的干扰读段，为准确鉴定与估计样本中物种成分与浓度奠定基础；提取物种的覆盖率、比对的gap得分、可变区得分三个特征，使用机器学习方法进行样本中物种成分的鉴定，实现物种浓度估计；
所述基于高通量测序数据的微生物成分及浓度检测方法具体包括：
第一步，病原菌数据库选取与NGS数据预处理，对于获取到的NGS测序样本，采用FastQC工具分析数据的质量情况；基于测序读段数据及参考序列，采用BWA进行比对，保留比对成功的读段序列；
第二步，构建读段-物种库的隶属关系体系，从比对成功的数据中量化碱基匹配或错配指示度、比对产生的gap与读段来自可变区的指示度信息，建立测序读段到物种库的隶属关系体系，通过过滤掉比对到物种下的低隶属分数的读段来矫正测序读段的比对结果，实现物种浓度估计；
第三步，过滤后的数据，提取物种下比对结果的位点覆盖率coverage、物种比对的gap分数记作gapscore、物种的可变区得分记作HVRScore，用于构建物种成分判别的分类器。
2.如权利要求1所述的基于高通量测序数据的微生物成分及浓度检测方法，其特征在于，第二步的物种的参考物种库集合为F＝{f1,f2,…,fH}，测序读段集合为R＝{r1,r2,…,rN}；rj的长度为Lj，rj的碱基序列为rj＝(rj[1],rj[2],…,rj[Lj])，rj，j＝1，2，...，N，比对到物种fi，i＝1，2，...，H后，比对形式记作fi＝(fi[1],fi[2],…,fi[Lj])；由碱基的测序质量值与错误率的关系Q＝-10*lg(qj[D])，D＝1，2，...，Lj，得rj的碱基测序错误率(qj[1],qj[2],…,qj[Lj])；
rj比对到fi产生的gap集合为X＝(x1,x2,…,xM)，其中，xb为产生gap的宽度，b＝1,2,...,M，M为gap的数量；若读段中碱基间的测序错误是独立的，则从碱基的匹配形式与测序错误角度而言，rj来源于fk的指示度衡量如下：
Figure FDA0004185719690000021
使用EMBOSS工具获取物种序列的可变区范围，rj比对到fi的可变区指示度量化为wij：
Figure FDA0004185719690000022
综合序列的比对形式、测序错误与可变区的指示度，将rj比对到fi的隶属分数定义如下：
Figure FDA0004185719690000023
将所有的测序数据比对到物种库后就可构建出读段-物种库的隶属关系体系；score(rj,fi)表示读段rj隶属于fi的程度。
3.如权利要求1所述的基于高通量测序数据的微生物成分及浓度检测方法，其特征在于，第三步的成分鉴定具体包括：
a)从物种的位点覆盖率而言，位点覆盖率越高则该物种存在的可能性越大，定义如下：
Figure FDA0004185719690000024
b)从空间特征而言，若物种比对后得到的gap集合为G＝{g1,g2,…,gn}，ga表示第a个gap的长度，a＝1，2，...，n，则物种比对结果下产生gap的得分如下：
Figure FDA0004185719690000031
c)从物种的生物特性而言，16S共包含9个可变区HVR，可变区下比对到的读段越多，则该物种存在的可能性越大；HVRi表示物种的第i个可变区，令x'i为物种比对结果下HVRi产生的gap长度，令yi为物种HVRi的长度，从矫正结果中量化物种的可变区得分如下：
Figure FDA0004185719690000032
通过构建800×的仿真数据，计算物种库中每个物种的特征向量coverage,gapscore,HVRScore，用作SVM分类器的训练集，实现对病变样本中物种成分的鉴定。
4.一种应用权利要求1～3任意一项所述基于高通量测序数据的微生物成分及浓度检测方法的基因序列数据检测系统。
5.一种应用权利要求1～3任意一项所述基于高通量测序数据的微生物成分及浓度检测方法的信息数据处理终端。
CN201910637328.7A 2019-07-15 2019-07-15 一种基于高通量测序数据的微生物成分及浓度检测方法 Active CN110517726B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910637328.7A CN110517726B (zh) 2019-07-15 2019-07-15 一种基于高通量测序数据的微生物成分及浓度检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910637328.7A CN110517726B (zh) 2019-07-15 2019-07-15 一种基于高通量测序数据的微生物成分及浓度检测方法

Publications (2)

Publication Number Publication Date
CN110517726A CN110517726A (zh) 2019-11-29
CN110517726B true CN110517726B (zh) 2023-07-04

Family

ID=68623632

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910637328.7A Active CN110517726B (zh) 2019-07-15 2019-07-15 一种基于高通量测序数据的微生物成分及浓度检测方法

Country Status (1)

Country Link
CN (1) CN110517726B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111881324B (zh) * 2020-07-30 2023-12-15 苏州工业园区服务外包职业学院 高通量测序数据通用存储格式结构、其构建方法及应用
CN112509701A (zh) * 2021-02-05 2021-03-16 中国医学科学院阜外医院 急性冠脉综合征的风险预测方法及装置
CN114373508B (zh) * 2022-01-24 2024-02-02 浙江天科高新技术发展有限公司 一种基于16S rDNA序列的菌种鉴定方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106868116A (zh) * 2017-01-24 2017-06-20 华南农业大学 一种桑树病原菌高通量鉴定及种属分类方法及其应用
CN107292123A (zh) * 2016-03-31 2017-10-24 苏州普瑞森基因科技有限公司 一种基于高通量测序的微生物群落组成的方法和装置
WO2018218150A1 (en) * 2017-05-26 2018-11-29 President And Fellows Of Harvard College Systems and methods for high-throughput image-based screening
CN109081867A (zh) * 2017-06-13 2018-12-25 北京大学 癌症特异性tcr及其分析技术和应用
CN109273053A (zh) * 2018-09-27 2019-01-25 华中科技大学鄂州工业技术研究院 一种高通量测序的微生物数据处理方法
CN109757103A (zh) * 2016-07-14 2019-05-14 百时美施贵宝公司 针对tim3的抗体及其用途

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111133115A (zh) * 2017-09-20 2020-05-08 瑞泽恩制药公司 用于其肿瘤携带高过客基因突变负荷的患者的免疫治疗方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107292123A (zh) * 2016-03-31 2017-10-24 苏州普瑞森基因科技有限公司 一种基于高通量测序的微生物群落组成的方法和装置
CN109757103A (zh) * 2016-07-14 2019-05-14 百时美施贵宝公司 针对tim3的抗体及其用途
CN106868116A (zh) * 2017-01-24 2017-06-20 华南农业大学 一种桑树病原菌高通量鉴定及种属分类方法及其应用
WO2018218150A1 (en) * 2017-05-26 2018-11-29 President And Fellows Of Harvard College Systems and methods for high-throughput image-based screening
CN109081867A (zh) * 2017-06-13 2018-12-25 北京大学 癌症特异性tcr及其分析技术和应用
CN109273053A (zh) * 2018-09-27 2019-01-25 华中科技大学鄂州工业技术研究院 一种高通量测序的微生物数据处理方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Error Correction and Clustering Algorithms for Next Generation Sequencing;Xiao Yang;《2011 IEEE International Symposium on Parallel and Distributed Processing Workshops and Phd Forum》;20110901;全文 *
噬菌体展示的生物信息学研究与应用;贺碧芳;《中国博士学位论文全文数据库 (基础科学辑)》;20180930(第9期);全文 *

Also Published As

Publication number Publication date
CN110517726A (zh) 2019-11-29

Similar Documents

Publication Publication Date Title
CN110517726B (zh) 一种基于高通量测序数据的微生物成分及浓度检测方法
US20190130999A1 (en) Latent Representations of Phylogeny to Predict Organism Phenotype
US12100479B2 (en) Systems and methods for metagenomic analysis
US20180137243A1 (en) Therapeutic Methods Using Metagenomic Data From Microbial Communities
JP2019537780A (ja) メタゲノム試料中の病原体の同定と抗生物質の特徴づけ
CN108064272B (zh) 用于类风湿性关节炎的生物标记物及其用途
CN115719616B (zh) 一种病原物种特异性序列的筛选方法及系统
CN113160882A (zh) 一种基于三代测序的病原微生物宏基因组检测方法
CN111599413B (zh) 一种测序数据的分类单元组分计算方法
BE1024766A1 (nl) Werkwijze voor het typeren van nucleïnezuur- of aminozuursequenties op basis van sequentieanalyse
CN112331268B (zh) 目标物种特有序列的获取方法及目标物种检测方法
WO2014136106A1 (en) Method and system for analyzing the taxonomic composition of a metagenome in a sample
JP2016518822A (ja) アセンブルされていない配列情報、確率論的方法、及び形質固有（ｔｒａｉｔ−ｓｐｅｃｉｆｉｃ）のデータベースカタログを用いた生物材料の特性解析
Bhat et al. A critical analysis of state-of-the-art metagenomics OTU clustering algorithms
Marchiori et al. Skraken: Fast and sensitive classification of short metagenomic reads based on filtering uninformative k-mers
Chandrasiri et al. CH-Bin: A convex hull based approach for binning metagenomic contigs
US20220259657A1 (en) Method for discovering marker for predicting risk of depression or suicide using multi-omics analysis, marker for predicting risk of depression or suicide, and method for predicting risk of depression or suicide using multi-omics analysis
CN114613440A (zh) 一种基于长读长测序数据的病原微生物分析方法
CN115331737A (zh) 一种分析肠道菌群中致病菌和量化菌群地域特征的方法
EP4305191A1 (en) Systems and methods for identifying microbial biosynthetic genetic clusters
Aldawiri et al. A Novel Approach for Mapping Ambiguous Sequences of Transcriptomes
Nguyen et al. Efficient and accurate OTU clustering with GPU-based sequence alignment and dynamic dendrogram cutting
Williams Application of Exact Alignments with an In-memory Core Gene Database for an Improved Metagenomic Taxonomic Classification
CN112634983B (zh) 病原物种特异pcr引物优化设计方法
RU2709815C1 (ru) Способ поиска молекулярных маркеров патологического процесса для дифференциальной диагностики, мониторинга и таргетной терапии

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