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

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

Info

Publication number
CN110517726A
CN110517726A CN201910637328.7A CN201910637328A CN110517726A CN 110517726 A CN110517726 A CN 110517726A CN 201910637328 A CN201910637328 A CN 201910637328A CN 110517726 A CN110517726 A CN 110517726A
Authority
CN
China
Prior art keywords
species
read
sequencing
score
data
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.)
Granted
Application number
CN201910637328.7A
Other languages
English (en)
Other versions
CN110517726B (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.)
Xian University of Electronic Science and Technology
Original Assignee
Xian University of Electronic Science and Technology
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 Xian University of Electronic Science and Technology filed Critical Xian University of Electronic Science and Technology
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

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 clusteringorders of magnitude faster than BLAST.Bioinformatics. 2010；26:2460–2461.][Edgar RC.UPARSE:highly accurate OTU sequences from microbial ampliconreads.Nat Methods.2013；10:996–998.]和其他一些基于相似 性聚类的算法[Al-GhalithGA,Montassier E,Ward HN,Knights D.NINJA-OPS: Fast Accurate Marker GeneAlignment Using Concatenated Ribosomes.PLoS Comput Biol.2016；12:e1004658.Albanese D,Fontana P,De Filippo C,Cavalieri D, Donati C.MICCA:acomplete and accurate software for taxonomic profiling of metagenomicdata.Sci Rep.2015；5:9743.Mahe F,Rognes T,Quince C,de Vargas C,DunthornM.Swarm:robust and fast clustering method for amplicon-based studies.PeerJ.2014；2:e593.Kopylova E,Noe L,Touzet H.SortMeRNA:fast and accuratefiltering of ribosomal RNAs in metatranscriptomic data.Bioinformatics.2012；28:3211–3217.]所取代，它们的速度与准确率也呈现大幅度提升并且大致彼此相 当[Al-Ghalith GA,Montassier E,Ward HN,Knights D.NINJA-OPS:Fast Accurate Marker GeneAlignment Using Concatenated Ribosomes.PLoS Comput Biol.2016； 12:e1004658.Kopylova E,Navas-Molina JA,Mercier C,Xu ZZ,Mahe′F,He Y,et al.Open-Source Sequence Clustering Methods Improve the State Of the Art.mSystems.2016；1]。
基于共享系统发育树的方法通过分析样本，将测序读段准确地放置在系统 发育树上，实现物种成分的识别。近年来，使用最大似然估计[Berger SA, Krompass D,Stamatakis A.Performance,accuracy,and Web server for evolutionary placementof short sequence reads under maximum likelihood.Syst Biol.2011； 60:291–302]、贝叶斯后验概率[Matsen FA,Kodner RB,Armbrust EV.pplacer: linear time maximum-likelihood and Bayesian phylogenetic placement of sequences onto a fixedreference tree.BMC Bioinformatics.2010；11:538.]或邻域连接[Price MN,Dehal PS,Arkin AP.FastTree:computing large minimum evolution trees with profilesinstead ofa distance matrix.Mol Biol Evol.2009；26:1641–1650.]的系统发 育算法相继被开发出来，但它们的物种识别灵敏度较低[Bazinet AL,Cummings MP.Acomparative evaluation of sequence classification programs.BMCBioinformatics.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,TsirigosA, Hugenholtz P,Rigoutsos I.Accurate phylogenetic classification of variable-length DNA fragments.Nat Methods.2007；4:63–72.Wang Q,Garrity GM,Tiedje JM,Cole JR.Naive Bayesian classifier for rapid assignment of rRNA 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的指示度衡量如下：
使用EMBOSS工具获取物种序列的可变区范围，rj比对到fi的可变区指示 度量化为wij：
综合序列的比对形式、测序错误与可变区的指示度，将rj比对到fi的隶属分 数定义如下：
将所有的测序数据比对到物种库后就可构建出读段-物种库的隶属关系体 系；score(rj,fi)表示读段rj隶属于fi的程度。
进一步，所述第三步的成分鉴定与浓度估计具体包括：
a)从物种的位点覆盖率而言，位点覆盖率越高则该物种存在的可能性越大， 定义如下：
b)从空间特征而言，若物种比对后得到的gap集合为G＝{g1,g2,…,gn}， gi表示第i个gap的长度，则物种比对结果下产生gap的得分如下：
c)从物种的生物特性而言，16S共包含9个可变区HVR，可变区下比对到 的读段越多，则该物种存在的可能性越大；HVRi表示物种的第i个可变区，令 xi为物种比对结果下HVRi产生的gap长度，令yi为物种HVRi的长度，从矫正结 果中量化物种的可变区得分如下：
通过构建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
测序读段比对到物种的参考序列的结果如表1所示，碱基的匹配与错配、 碱基的插入与缺失(统称为gap)，可直观地指示读段rj来自于物种fi的的可能 性。比对结果中碱基的错配数量越少，产生gap的个数与gap的宽度越小，则该 读段比对到该物种的可信程度越高。令rj比对到fi产生的gap集合为X＝(x1, x2,…,xM)，其中，xi为产生gap的宽度，M为gap的数量。若读段中碱基间的 测序错误是独立的，则从碱基的匹配形式与测序错误角度而言，rj来源于fk的 指示度衡量如下：
由于微生物物种间的相似性极高，所以区分物种类型的最重要的信号就在 于物种的可变区。对于一个多比对读段，若该读段比对到物种A的可变区内， 同时比对到其他物种B、C的保守区，则基于物种的生物特性，认为该读段更有 可能来源于物种A。读段的多比对形式如图3所示，本发明使用EMBOSS工具 获取物种序列的可变区范围，rj比对到fi的可变区指示度量化为wij：
综合序列的比对形式、测序错误与可变区的指示度，将rj比对到fi的隶属分 数定义如下：
将所有的测序数据比对到物种库后就可构建出读段-物种库的隶属关系体 系，该体系可用于比对结果的矫正。score(rj,fi)表示读段rj隶属于fi的程度，该 值越大则读段来自该物种的可能性越高，也就是比对的正确率越高。病变样本 中可能含有未知的微生物或人类基因的干扰，来自这些干扰物种的读段虽然会 错误地比对到物种库的某些物种上，但是它们往往存在很多的碱基错配、大量 的gap，并且以很低的概率落在可变区内，以致读段到该物种的隶属分数很低。 本发明通过大量仿真实验探索，将隶属分数阈值设置为0.44，过滤掉错误比对 的读段，排除多比对读段的错误比对部分，实现比对结果的矫正，这对物种的 成分鉴定极为重要。
(3)成分鉴定与浓度估计
针对(2)中过滤后的数据，提取物种下比对结果的位点覆盖率coverage、 物种比对的gap分数(记作gapscore)、物种的可变区得分(记作HVRScore)， 用于构建物种成分判别的分类器。
a)从物种的位点覆盖率而言，位点覆盖率越高则该物种存在的可能性越大， 其定义如下：
b)从空间特征而言，若物种比对后得到的gap集合为G＝{g1,g2,…,gn}， gi表示第i个gap的长度，则物种比对结果下产生gap的得分如下：
c)从物种的生物特性而言，16S共包含9个可变区(HighVariable Region, HVR)，可变区下比对到的读段越多，则该物种存在的可能性越大。令HVRi表 示物种的第i个可变区，令xi为物种比对结果下HVRi产生的gap长度，令yi为 物种HVRi的长度，从矫正结果中量化物种的可变区得分如下：
通过构建800×的仿真数据，计算物种库中每个物种的特征向量(coverage,gapscore,HVRScore)，用作SVM分类器的训练集，可实现对病变样本中物种成 分的鉴定。
由于在隶属分数过滤过程中，可能会错误地过滤掉某些多比对读段，这对 物种的浓度估计会造成干扰，所以需要回收过滤掉的多比对读段，并将这些读 段重新比对到已鉴定出的物种上，统计每个物种下比对到的读段量，计算它们 在测序数据中的相对浓度，作为物种浓度的估计结果。
下面结合实验对本发明的技术效果作详细的描述。
本发明使用仿真工具ART生成800×，Insertion/deletion均为0.1的11组物 种种类、浓度、干扰程度均不同仿真数据，数据内容如表1：
表1
[软件来源ARThttps://www.niehs.nih.gov/research/resources/software/biostatistics/art/]
本发明命名为PGMicroD工具，与现存的五种工具：Bwa、Harp[Kessner D, TurnerT L,Novembre J.Maximum Likelihood Estimation of Frequencies of KnownHaplotypes from Pooled Sequence Data[J].Molecular Biology and Evolution,2013,30(5):1145-1158.]、Kallisto[Bray N L,Pimentel H,Melsted,Páll,et al. Near-optimal probabilistic RNA-seq quantification[J].Nature Biotechnology,2016.]、Karp[Reppell M,Novembre J,Mchardy A C.Using pseudoalignment and base qualityto accurately 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 of Rhizobia[J].PLOS ONE, 2012,7.]对仿真数据中物种类型的类型与浓度检测结果如图4： 其中，ti是物种的估计浓度，τi是物种的真实浓度；
MaxRE＝max|(ti-τi)|，其中，ti是物种的估计浓度，τi是物种的真实浓度。
PGMicroD与其他五种工具对来自血液、尿液、脑脊液的25组临床样本的检测 结果如图5；在仿真数据和真实数据上的作用效果，PGMicroD均表现出较佳状 态。
本实验处列举六种工具在第七组仿真数据上的物种类型与浓度检测结果如 下所示，PGMicroD的结果明显优于其他五种工具。
表2真实物种类型编号与浓度
表3Harp检测的物种类型编号与浓度
物种编号 浓度
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检测的物种类型编号与浓度
表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检测的物种类型编号与浓度
表7 Mothur检测的物种类型编号与浓度
物种编号 浓度
61 0.13903
62 0.124266
63 0.040484
64 0.114136
65 0.038944
66 0.006653
67 0.005883
68 0.005245
69 0.003181
70 0.014302
71 0.026068
72 0.010465
73 0.056945
74 0.010673
75 0.00765
76 0.00065
77 0.004481
78 0.066614
79 6.31E-06
80 0.005037
表8 PGMicroD检测的物种类型编号与浓度
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发 明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明 的保护范围之内。

Claims (6)

1.一种基于高通量测序数据的微生物成分及浓度检测方法，其特征在于，所述基于高通量测序数据的微生物成分及浓度检测方法量化碱基匹配或错配指示度、比对产生的gap与特定可变区指示度信息；建立基于NGS数据的读段-物种库的隶属关系体系；通过设置隶属分数阈值过滤掉潜在的干扰读段，为准确鉴定与估计样本中物种成分与浓度奠定基础；提取物种的覆盖率、比对的gap得分、可变区得分三个特征，使用机器学习方法进行样本中物种成分的鉴定。
2.如权利要求1所述的基于高通量测序数据的微生物成分及浓度检测方法，其特征在于，所述基于高通量测序数据的微生物成分及浓度检测方法具体包括：
第一步，病原菌数据库选取与NGS数据预处理，对于获取到的NGS测序样本，采用FastQC工具分析数据的质量情况；基于测序读段数据及参考序列，采用BWA进行比对，保留比对成功的读段序列；
第二步，构建读段-物种库的隶属关系体系，从比对成功的数据中量化碱基匹配或错配指示度、比对产生的gap与读段来自可变区的指示度等信息，建立测序读段到物种库的隶属关系体系，通过过滤掉比对到物种下的低隶属分数的读段来矫正测序读段的比对结果；
第三步，过滤后的数据，提取物种下比对结果的位点覆盖率coverage、物种比对的gap分数记作gapscore、物种的可变区得分记作HVRScore，用于构建物种成分判别的分类器。
3.如权利要求2所述的基于高通量测序数据的微生物成分及浓度检测方法，其特征在于，所述第二步的物种的参考数据库集合为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的指示度衡量如下：
使用EMBOSS工具获取物种序列的可变区范围，rj比对到fi的可变区指示度量化为wij：
综合序列的比对形式、测序错误与可变区的指示度，将rj比对到fi的隶属分数定义如下：
将所有的测序数据比对到物种库后就可构建出读段-物种库的隶属关系体系；score(rj,fi)表示读段rj隶属于fi的程度。
4.如权利要求2所述的基于高通量测序数据的微生物成分及浓度检测方法，其特征在于，所述第三步的成分鉴定与浓度估计具体包括：
a)从物种的位点覆盖率而言，位点覆盖率越高则该物种存在的可能性越大，定义如下：
b)从空间特征而言，若物种比对后得到的gap集合为G＝{g1,g2,…,gn}，gi表示第i个gap的长度，则物种比对结果下产生gap的得分如下：
c)从物种的生物特性而言，16S共包含9个可变区HVR，可变区下比对到的读段越多，则该物种存在的可能性越大；HVRi表示物种的第i个可变区，令xi为物种比对结果下HVRi产生的gap长度，令yi为物种HVRi的长度，从矫正结果中量化物种的可变区得分如下：
通过构建800×的仿真数据，计算物种库中每个物种的特征向量coverage,gapscore,HVRScore，用作SVM分类器的训练集，实现对病变样本中物种成分的鉴定。
5.一种应用权利要求1～4任意一项所述基于高通量测序数据的微生物成分及浓度检测方法的基因序列数据检测系统。
6.一种应用权利要求1～4任意一项所述基于高通量测序数据的微生物成分及浓度检测方法的信息数据处理终端。
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 true CN110517726A (zh) 2019-11-29
CN110517726B 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)

Cited By (3)

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

Citations (7)

* 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 华中科技大学鄂州工业技术研究院 一种高通量测序的微生物数据处理方法
US20190087538A1 (en) * 2017-09-20 2019-03-21 Regeneron Pharmaceuticals, Inc. Immunotherapy Methods For Patients Whose Tumors Carry A High Passenger Gene Mutation Burden
CN109757103A (zh) * 2016-07-14 2019-05-14 百时美施贵宝公司 针对tim3的抗体及其用途

Patent Citations (7)

* 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及其分析技术和应用
US20190087538A1 (en) * 2017-09-20 2019-03-21 Regeneron Pharmaceuticals, Inc. Immunotherapy Methods For Patients Whose Tumors Carry A High Passenger Gene Mutation Burden
CN109273053A (zh) * 2018-09-27 2019-01-25 华中科技大学鄂州工业技术研究院 一种高通量测序的微生物数据处理方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
XIAO YANG: "Error Correction and Clustering Algorithms for Next Generation Sequencing", 《2011 IEEE INTERNATIONAL SYMPOSIUM ON PARALLEL AND DISTRIBUTED PROCESSING WORKSHOPS AND PHD FORUM》 *
贺碧芳: "噬菌体展示的生物信息学研究与应用", 《中国博士学位论文全文数据库 (基础科学辑)》 *

Cited By (5)

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

Also Published As

Publication number Publication date
CN110517726B (zh) 2023-07-04

Similar Documents

Publication Publication Date Title
Ebrahimie et al. Hierarchical pattern recognition in milking parameters predicts mastitis prevalence
CN110517726A (zh) 一种基于高通量测序数据的微生物成分及浓度检测方法
CN104603283B (zh) 确定异常状态相关生物标志物的方法及系统
Gupta et al. Identification of discriminatory antibiotic resistance genes among environmental resistomes using extremely randomized tree algorithm
Spergser et al. Matrix-assisted laser desorption ionization–time of flight mass spectrometry is a superior diagnostic tool for the identification and differentiation of mycoplasmas isolated from animals
Faron et al. Multicenter evaluation of the Bruker MALDI Biotyper CA system for the identification of clinical aerobic gram-negative bacterial isolates
CA2304876A1 (en) Methods for classifying samples and ascertaining previously unknown classes
Ulrich et al. ReadBouncer: precise and scalable adaptive sampling for nanopore sequencing
CN108064272B (zh) 用于类风湿性关节炎的生物标记物及其用途
CN110111843B (zh) 对核酸序列进行聚类的方法、设备及存储介质
WO2018160899A1 (en) Systems and methods for metagenomic analysis
CN110459264A (zh) 基于梯度增强决策树预测环状rna与疾病相关性的方法
CN115064215B (zh) 一种通过相似度进行菌株溯源及属性鉴定的方法
CN112331268B (zh) 目标物种特有序列的获取方法及目标物种检测方法
US20130045878A1 (en) Process for identification of pathogens
Teng et al. MALDI-TOF MS for identification of Tsukamurella species: Tsukamurella tyrosinosolvens as the predominant species associated with ocular infections
CN115171792A (zh) 一种毒力因子和抗生素抗性基因的混合预测方法
US20180330044A1 (en) Methods Associated With A Database That Stores A Plurality Of Reference Genomes
ES2659268T3 (es) Detección de mezclas en la identificación de microbios por espectrometría de masas
Centeno-Martinez et al. The bovine nasal fungal community and associations with bovine respiratory disease
Zhang et al. MetaDomain: a profile HMM-based protein domain classification tool for short sequences
Uricaru et al. YOC, A new strategy for pairwise alignment of collinear genomes
CN105567831B (zh) 一种食品微生物定性与定量的检测方法
Sun et al. Eliminate false positives in metagenomic profiling based on type IIB restriction sites
Marić et al. Approaches to metagenomic classification and assembly

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