CN115064215B - 一种通过相似度进行菌株溯源及属性鉴定的方法 - Google Patents

一种通过相似度进行菌株溯源及属性鉴定的方法 Download PDF

Info

Publication number
CN115064215B
CN115064215B CN202210991537.3A CN202210991537A CN115064215B CN 115064215 B CN115064215 B CN 115064215B CN 202210991537 A CN202210991537 A CN 202210991537A CN 115064215 B CN115064215 B CN 115064215B
Authority
CN
China
Prior art keywords
sequence
cluster
similarity
genome
strains
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
CN202210991537.3A
Other languages
English (en)
Other versions
CN115064215A (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.)
Peking University Peoples Hospital
Original Assignee
Peking University Peoples Hospital
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 Peking University Peoples Hospital filed Critical Peking University Peoples Hospital
Priority to CN202210991537.3A priority Critical patent/CN115064215B/zh
Publication of CN115064215A publication Critical patent/CN115064215A/zh
Application granted granted Critical
Publication of CN115064215B publication Critical patent/CN115064215B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10Sequence alignment; Homology search
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/30Detection of binding sites or motifs
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00ICT programming tools or database systems specially adapted for bioinformatics

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Biophysics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Bioethics (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Chemical & Material Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Molecular Biology (AREA)
  • Data Mining & Analysis (AREA)
  • Genetics & Genomics (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提供一种通过相似度进行菌株溯源及属性鉴定的方法,其包含如下步骤:a)将样本中菌株的测序序列或组装序列和聚类泛基因组数据库中每个簇的泛基因组进行序列比对;b)根据比对上序列的覆盖度筛选比对结果;c)根据筛选的比对结果计算簇相似度,选取簇相似度最大的L个簇作为候选簇;d)从每个所述候选簇中选取株相似度最大的M个原始菌株作为候选菌株;e)从所述候选菌株中,选取株相似度最大的N个原始菌株,作为来源菌株;f)根据所述来源菌株获取所述样本中菌株的来源信息及属性信息。采用此方法可快速、准确地鉴定出病原微生物的属性和来源。

Description

一种通过相似度进行菌株溯源及属性鉴定的方法
技术领域
本发明属于生信分析领域,具体涉及一种通过相似度进行菌株溯源及属性鉴定的方法。
背景技术
传统药敏检测方法如纸片扩散法和肉汤稀释法,耗时长,且需要先将细菌分离培养为纯培养物。传统的溯源或同源性分析往往依赖于从头开始的系统发育分析,每得到一株新纯培养物,对菌株序列进行一次新的系统发育分析,观察该序列在系统发育树上的位置。此类方法由于需要如细菌的分离培养、基因组提取、基因组组装、基因功能注释、核心基因组分析、系统发育树分析等步骤,所需时间长,且比对结果严重依赖于各实验室或医院纳入的菌株范围。
宏基因组学规避了对样本的微生物的分离培养,可进行样本直接检测,报告所有检出序列的病原体,其检测速度之快为危急感染患者的临床诊断提供新保障。虽然能通过宏基因组能得到所有检出序列的病原体,但是无法得知该病原体的来源,基因组属性等。这可能会错过院内暴发或者社区暴发,导致更多感染性事件的发生。无法得知基因组属性还可能导致用药错误或药物滥用。
有鉴于此,提出本发明。
发明内容
本发明提供一种通过聚类泛基因组和相似度快速鉴定菌株属性及其来源的方法,采用此方法可快速、准确地鉴定出病原微生物的属性和来源,具体包括以下实施方式:
实施方式1. 一种通过相似度进行菌株溯源及属性鉴定的方法,其特征在于,包含如下步骤:
a)序列比对:将样本中菌株的测序序列或组装序列和聚类泛基因组数据库中每个簇的泛基因组进行序列比对,其中所述聚类泛基因组数据库中包含泛基因组序列-原始序列-原始菌株关联关系,和任选的属性信息;
b)比对结果筛选:对a)得到的序列比对结果,根据比对上序列的覆盖度筛选比对结果;
c)计算相似度,筛选候选簇:根据b)筛选的比对结果计算样本中菌株与某簇的簇相似度,选取簇相似度最大的L个簇作为样本中菌株所在的候选簇,其中,L为正整数;
d)筛选候选菌株:从每个所述候选簇中选取株相似度最大的M个原始菌株作为与样本中菌株相似的候选菌株,其中,M为正整数,所述株相似度用以衡量样本中菌株与所述聚类泛基因组数据库中关联的原始菌株的相似性;
e)确定来源菌株:从所述候选菌株中,选取株相似度最大的N个原始菌株,作为样本中菌株的来源菌株,其中,N为正整数;
f)获取来源信息及属性信息:根据所述来源菌株获取所述样本中菌株的来源信息及属性信息。
实施方式2. 根据实施方式1所述的方法,其特征在于,所述覆盖度为某条测序序列或组装序列与该簇泛基因组某条序列比对上的长度占所述该簇泛基因组某条序列的长度的百分比。
实施方式3. 根据实施方式1所述的方法,其特征在于,所述覆盖度为某条测序序列或组装序列与该簇泛基因组某条序列比对上的长度占所述某条测序序列或组装序列的长度的百分比。
实施方式4. 根据实施方式1所述的方法,其特征在于,所述簇相似度为比对上某簇泛基因组所有序列的相似度的平均值。
实施方式5. 根据实施方式1所述的方法,其特征在于,所述筛选候选簇如下进行:
所述L选取5个至10个,或者
选取簇相似度达到最大簇相似度的99%的簇作为样本中菌株所在的候选簇。
实施方式6. 根据实施方式1所述的方法,其特征在于,所述株相似度为比对上某簇泛基因组中可以关联到某原始菌株的所有序列相似度之和。
实施方式7. 根据实施方式1所述的方法,其特征在于,所述聚类泛基因组数据库按以下步骤构建:
步骤1)序列相似性计算:计算微生物菌株基因组序列的相似性;
步骤2)序列聚类:根据序列的相似性计算结果进行聚类,将相似序列聚为同一簇;
步骤3)构建聚类泛基因组数据库:对于聚类之后的每一簇,构建该簇所有菌株的泛基因组,作为该簇菌株的序列特征,每簇泛基因组序列包含序列ID以及与原始菌株对应信息,进而构建成物种的聚类泛基因组数据库;其中,所述的微生物菌株基因组序列包括收集或自建的菌株序列数据及属性信息,和/或收集的公共数据库中菌株序列数据及属性信息。
实施方式8. 根据实施方式7所述的方法,其特征在于,所述的微生物菌株基因组序列为二代组装或三代组装或二代和三代混合组装的基因组;所述公共数据库采用NCBI中的Refseq数据库。
实施方式9. 根据实施方式7所述的方法,其特征在于,所述方法还包括:
根据微生物菌株基因组序列及属性信息,整理所有菌株、序列及其属性信息,构建菌株-序列-属性信息表;
根据构建的物种聚类泛基因组数据库和菌株-序列-属性信息表构建聚类泛基因组的序列ID-菌株-序列-属性信息表。
实施方式10. 根据实施方式7所述的方法,其特征在于,所述的属性信息包括菌株ID、菌株名称、菌株分类、序列ID、收集地点、提交地点、收集时间、提交时间、基因组信息、耐药基因、MLST分型、KL分型、药敏结果及分类属性信息中的一个或多个。
实施方式11. 根据实施方式7所述的方法,其特征在于,
所述的步骤1)中,取ANI距离作为序列相似性的衡量标准,或通过提取核心基因组构建系统进化树或bac120基因集构建系统进化树,根据进化距离作为序列相似性的衡量标准;
所述的步骤2)中,序列聚类采用层次聚类法。
实施方式12. 根据实施方式7所述的方法,其特征在于,所述方法还包括:构建耐药和毒力基因数据库:从耐药数据库和毒力基因数据库下载序列及相关数据,作为本地耐药和毒力基因数据库。
实施方式13. 根据实施方式7至12任一项所述的方法,其特征在于,所述方法还包括:收集微生物菌株基因组序列实时更新数据库。
实施方式14. 根据实施方式1所述方法,其特征在于,所述测序序列包括质控测序序列,其通过如下步骤获取:
数据预处理:将样本中菌株的测序数据剔除接头、低质量序列和长度过短序列,得到预处理后的数据;
去宿主处理:将获取的所述预处理后的数据与宿主基因组比对,去除比对上宿主基因的序列,从而获得非宿主序列,即为所述质控测序序列;
所述测序序列组装后即获得所述组装序列。
实施方式15. 根据实施方式1所述的方法,其特征在于,所述样本中菌株的来源信息及属性信息包括耐药基因和/或毒力基因数据库比对结果,所述耐药基因和/或毒力基因数据库比对结果通过如下方法获得:
将样本中菌株的测序序列或组装获取的序列与耐药基因和/或毒力基因数据库进行序列比对,序列比对结果中根据比对上序列的基因相似度进行筛选,判断样本中菌株携带的耐药基因和/或毒力基因。
实施方式16. 根据实施方式1所述的方法,其特征在于,所述样本中菌株包括一种或多种病原微生物。
实施方式17.一种鉴定病原微生物暴发的方法,其特征在于,对特定区域内在不同时间获取的样本中菌株的测序序列或组装序列,按实施方式1至16中任一项所述的方法进行菌株溯源和属性鉴定,从而鉴定病原微生物暴发。
实施方式18. 一种电子设备,其特征在于,包括:存储器、与所述存储器连接的处理器,及存储在所述存储器上并可在所述处理器上运行的计算机程序,所述处理器运行所述计算机程序时执行实施方式1-17任一项所述的方法。
本发明提供一种通过聚类泛基因组和相似度快速鉴定菌株属性及其来源的方法,与现有技术相比,本发明提供的技术方案能在宏基因组学报告出病原体的半小时或一小时内,报告出与该病原体相近的病原体的基因组属性、分型和来源,药敏结果等基本信息,同时报告出该病原体所拥有的耐药基因和/或毒力基因,从而能根据之前的药敏结果及时精确用药,对症治疗,拯救危及患者,减少抗生素的滥用。另外,本申请的技术方案还能够通过报告该病原体的属性和可能来源,进而提示是否发生院内或社区暴发,尽早阻断传播链,降低院内或社区传播风险。本发明的方法适用于一代、二代、三代测序数据,尤其适用于三代测序,无论是纯病原微生物或者是三代宏基因组测序得到的数据。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案,下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图是本发明的一些实施方式,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明的基础流程图。
具体实施方式
下面将结合附图对本发明的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
以下术语或定义仅仅是为了帮助理解本发明而提供。这些定义不应被理解为具有小于本领域技术人员所理解的范围。
除非在下文中另有定义,本发明具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解,但仍然阐述以下定义以更好地解释本发明。
如本发明中所使用,术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的,且不排除其它未列举的元素或方法步骤。术语“由…组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案,这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。
在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”,“所述”,包括该名词的复数形式。
本发明中的术语“大约”、“大体”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10%,优选±5%。
此外,说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类,是用于区分相似的元素,不是描述顺序或时间次序必须的。应理解,如此应用的术语在适当的环境下可互换,并且本发明描述的实施方案能以不同于本发明描述或举例说明的其它顺序实施。
本申请中,如无相反说明,“测序序列”和“测序read”具有相同的含义,“组装序列”和“拼接序列”具有相同的含义。
本发明一方面提供一种通过相似度进行菌株溯源及属性鉴定的方法,其特征在于,包含如下步骤:
a)序列比对:将样本中菌株的测序序列或组装序列和聚类泛基因组数据库中每个簇的泛基因组进行序列比对,其中所述聚类泛基因组数据库中包含泛基因组序列-原始序列-原始菌株关联关系,和任选的属性信息;
b)比对结果筛选:对a)得到的序列比对结果,根据比对上序列的覆盖度筛选比对结果;
c)计算相似度,筛选候选簇:根据b)筛选的比对结果计算样本中菌株与某簇的簇相似度,选取簇相似度最大的L个簇作为样本中菌株所在的候选簇,其中,L为正整数;
d)筛选候选菌株:从每个所述候选簇中选取株相似度最大的M个原始菌株作为与样本中菌株相似的候选菌株,其中,M为正整数,所述株相似度用以衡量样本中菌株与所述聚类泛基因组数据库中关联的原始菌株的相似性;
e)确定来源菌株:从所述候选菌株中,选取株相似度最大的N个原始菌株,作为样本中菌株的来源菌株,其中,N为正整数;
f)获取来源信息及属性信息:根据所述来源菌株获取所述样本中菌株的来源信息及属性信息。从而完成菌株溯源及属性鉴定。
本申请所述的属性信息包括菌株ID、菌株名称、菌株分类、序列ID、收集地点、提交地点、收集时间、提交时间、基因组信息、耐药基因、MLST分型、KL分型、药敏结果及分类属性信息中的一个或多个。
在一些实施方式中,所述覆盖度为某条测序序列或组装序列与该簇泛基因组某条 序列比对上的长度占所述该簇泛基因组某条序列的长度的百分比。某条序列的覆盖度如下计算:
式中, 为该read或组装序列与泛基因组序列seq比对上的碱基长 度;
为该read或组装序列比对上的泛基因组序列seq的碱基长度。
优选地,选取>80%的比对结果。
在一些实施方式中,所述覆盖度为某条测序序列或组装序列与该簇泛基因组某条 序列比对上的长度占所述某条测序序列或组装序列的长度的百分比。某条序列的所述覆盖 度如下计算:
式中,为该read或组装序列与泛基因组序列seq比对上的碱基长度;
为该read的碱基长度。
在一些实施方式中,所述簇相似度为比对上某簇泛基因组所有序列的相似度的平 均值。所述簇相似度所述簇相似度按以下方法计算:
其中, 为某条测序read或拼接序列与该簇泛基因组某条序列seq比对 上的长度占该seq的覆盖度;
为某条测序read或组装序列与该簇泛基因组某条序列seq的一 致性;
为某条测序read或组装序列与该簇泛基因组某条序列seq的相似度;
为与该簇泛基因组比对上的第i个测序read或第i个组装序列;
为比对上该簇泛基因组的某条序列seq的read或组装序列的个数;
为与泛基因组某序列seq的相似度,即所有比对上该簇泛基因组某序列seq 的测序read或组装序列与该序列的相似度的最大值;
为该簇泛基因组中比上的所有seq的数目;
为聚类泛基因组数据库中簇的总数。
在一些实施方式中,所述筛选候选簇如下进行:所述L选取5个至10个,或者选取簇相似度达到最大簇相似度的99%的簇作为样本中菌株所在的候选簇。
在一些实施方式中,所述株相似度为比对上某簇泛基因组中可以关联到某原始菌 株的所有序列相似度之和。优选地,所述株相似度如下计算:
式中,为与泛基因组某序列seq的相似度,即所有比对上该簇泛基因组某序 列seq的测序read或组装序列与该序列的相似度的最大值;
为比对上的簇内泛基因组序列的个数。
优选地,选取株相似度最大的Top5的菌株作为候选菌株。
在一些实施方式中,所述聚类泛基因组数据库按以下步骤构建:
步骤1)序列相似性计算:计算微生物菌株基因组序列的相似性;
步骤2)序列聚类:根据序列的相似性计算结果进行聚类,将相似序列聚为同一簇;
步骤3)构建聚类泛基因组数据库:对于聚类之后的每一簇,构建该簇所有菌株的泛基因组,作为该簇菌株的序列特征,每簇泛基因组序列包含序列ID以及与原始菌株对应信息,进而构建成物种的聚类泛基因组数据库;其中,所述的微生物菌株基因组序列包括收集或自建的菌株序列数据及属性信息,和/或收集的公共数据库中菌株序列数据及属性信息。
在一些实施方式中,所述聚类泛基因组数据库的构建方法如下:
1)构建菌株代表性基因组序列库:收集多个地区不同医院,尤其是同一医院或地区短时间内的暴发菌株并测序,作为自收集数据库;收集公共数据库(如Refseq)中不同地区或国家上传的序列数据及相关信息;2)构建序列/菌株属性信息表:根据所述的菌株代表性基因组序列库及相关信息,整理所有序列/菌株的属性信息,优先的为菌株ID、序列ID、收集或提交地点、收集时间、提交时间、基因组信息、耐药基因、MLST分型、KL分型、药敏结果及分类等属性信息,构建序列/菌株属性信息表;3)序列相似性计算:计算微生物菌株基因组序列的相似性,即计算所有步骤1)中菌株基因组序列的相似性,优选的,取ANI距离作为序列相似度的衡量标准,也可以用提取核心基因组构建系统进化树或bac120基因集构建系统进化树计算序列之间的相似性;4)序列聚类:根据序列的相似性计算结果进行聚类,将相似序列聚为同一簇,优选的采用层次聚类法中的最大距离法聚类;5)构建聚类泛基因组数据库:对于聚类之后的每一簇,不重复地提取该簇所有菌株的泛基因组,作为该簇菌株的序列特征,每簇泛基因组序列片段包含序列ID以及与原始菌株对应信息,进而构建成物种的聚类泛基因组数据库;6)构建属性信息和泛基因组关联信息:根据微生物菌株基因组序列及属性信息,整理所有菌株、序列及其属性信息,构建菌株-序列-属性信息表;根据构建的物种聚类泛基因组数据库和菌株-序列-属性信息表构建聚类泛基因组的序列ID-菌株-序列-属性信息表,进而构建物种的聚类泛基因组及属性数据库。7)构建耐药和毒力基因数据库:从耐药数据库和毒力基因数据库下载序列及相关数据,作为本地耐药和毒力基因数据库。
在一些实施方式中,所述测序序列包括质控测序序列,其通过如下步骤获取:
数据预处理:将样本中菌株的测序数据剔除接头、低质量序列和长度过短序列,得到预处理后的数据;
去宿主处理:将获取的所述预处理后的数据与宿主基因组比对,去除比对上宿主基因的序列,从而获得非宿主序列,即为所述质控测序序列;
所述测序序列组装后即获得所述组装序列。
在一些实施方式中,所述样本中菌株的来源信息及属性信息包括耐药基因和/或 毒力基因数据库比对结果,所述耐药基因和/或毒力基因数据库比对结果通过如下方法获 得:将样本中菌株的测序序列或组装获取的序列与耐药基因和/或毒力基因数据库进行序 列比对,序列比对结果中根据比对上序列的基因相似度进行筛选,判断样本中菌株 携带的耐药基因和/或毒力基因。在一些优选的实施方式中,将下机的fastq数据与耐药基 因和/或毒力基因数据库相比较,以基因相似度判定数据中是否含有该基因,所述 基因相似度如下计算:
其中为该read或组装序列比对上长度占耐药基因和毒力基因数据库 中基因的覆盖度;
为该read或组装序列与耐药基因和毒力基因数据库中基因比对上 的碱基长度;
为该read或组装序列比对上的耐药基因和毒力基因数据库中基因的 碱基长度;
为某条测序read或组装序列比对上的耐药基因和毒力基因数据 库中基因的一致性。
优选地,选取基因相似度>90%的基因,综合判断样本中菌株携带的耐药基 因和/或毒力基因。
在一些实施方式中,所述样本中菌株包括一种或多种病原微生物。
本申请另一方面提供一种鉴定病原微生物暴发的方法,其特征在于,对特定区域内在不同时间获取的样本中菌株的测序序列或组装序列,按前述任一方法进行菌株溯源和属性鉴定,从而鉴定病原微生物暴发。进一步地,选择同一医院或地区短时间内的暴发菌株测序后组装基因组序列及时按前述方法处理实时更新数据库。
本申请还提供一种电子设备,其特征在于,包括:存储器、与所述存储器连接的处理器,及存储在所述存储器上并可在所述处理器上运行的计算机程序,所述处理器运行所述计算机程序时执行前述任一方法。
本发明还提供一种计算机存储介质,所述计算机存储介质存储有计算机程序,所述计算机程序包括程序指令,所述程序指令当被处理器执行时,执行如上任一项所述的方法。所述计算机存储介质还存储有本专利构建的暴发数据库,当程序指令被执行时,读取数据库中的数据,进行后续分析。
在一些实施方式中,所述的菌株溯源及属性鉴定的方法包括如图1所示:
1)数据预处理:将测序数据剔除接头、低质量序列和长度过短序列read,测序数据 优选的指一代、二代、三代、四代测序数据;优选的,二代数据先拼接,得到组装数据;2)去宿 主处理:将上述获取测序或组装序列与宿主基因组比对,去除比对上宿主的序列,从而获得 非宿主序列;3)序列比对:将2)获取的序列和聚类泛基因组数据库中的每个簇的泛基因组 进行序列比对;4)比对结果筛选:将3)对比的结果中根据对上序列的覆盖度筛选比对结果,优选的>80%。5)计算相似度,筛选候选簇:根据4)筛选的 比对序列结果计算该菌株与某簇的簇相似度,选取簇相似度最大的L (L≥1,L为自然数)个簇作为该菌株所在的候选簇,优选的候选簇取Top5或簇相似度大于最 大相似度 * 99%内的簇。6)筛选候选菌株:将每个候选簇中选取株 相似度最大的M(M≥1,M为自然数)个作为与样本相似的候选菌株,即候选簇中取最大的对应的菌株,优选的选取株相似度最大的Top5的菌株。7)确定 相似菌株:依据候选序列根据病原微生物数据库中关联的泛基因组序列-原始序列信息关 系,获取最相似的原始菌株信息,进而获取该样本可能的来源、MLST分型、KL分型、耐药基因 等信息;8)综合判断样本中菌株携带的耐药基因和毒力基因:将下机的fastq数据与耐药基 因和毒力基因数据库相比较,以基因相似度判定数据中是否含有该基因;将2)获取 的序列read或组装序列与耐药基因和毒力基因数据库中的每个基因进行序列比对。对比的 结果中根据对上序列的基因相似度筛选比对结果,优选的候选>90%的基 因。9)将7)获取到的这些序列的属性和来源以及8)获取到的耐药毒力基因数据库比对的结 果作为该微生物基因组的属性和可能来源。
进一步的,步骤3)中选定的结果可以改变限制,例如相似度可以设定为大于85%,或者90%,或者95%,或者98%。
本申请还公开了一种宏基因组微生物鉴定系统,所述系统包括如下模块:1)数据 预处理:将测序数据剔除接头、低质量序列和长度过短序列read,测序数据优选的指一代、 二代、三代、四代测序数据;优选的,二代数据先拼接,得到组装数据;2)去宿主处理:将上属 获取测序或组装序列与宿主基因组比对,去除比对上宿主的序列,从而获得非宿主序列;3) 序列比对:将2)获取的序列和聚类泛基因组数据库中的每个簇的泛基因组进行序列比对; 4)比对结果筛选:将3)对比的结果中根据对上序列的覆盖度筛选 比对结果,优选的,>80%。5)计算相似度,筛选候选簇:根据4)筛选的比对序列结 果计算该菌株与某簇的相似度,选取相似度最大的N(N≥1,N为自然 数)个簇作为该菌株所在的候选簇,优选的候选簇取Top5或相似度大于最大相似度 *99%内的簇。6)筛选候选菌株:将每个候选簇中选取相似度最大的N(N≥1,N为自然数)个作为与样本相似的候选菌株,即候选簇中取最大的对应的菌株,优选的选取相似度最大的Top5的菌株。7)确定相 似菌株:依据候选序列根据病原微生物数据库中关联的泛基因组序列-原始序列信息关系, 获取最相似的原始菌株信息,进而获取该样本可能的来源、MLST分型、KL分型、耐药基因等 信息;8)综合判断样本中菌株携带的耐药基因和毒力基因:将下机的fastq数据与耐药基因 和毒力基因数据库相比较,以相似度判定数据中是否含有该基因;将2)获取的序列 read或组装序列与耐药基因和毒力基因数据库中的每个基因进行序列比对。对比的结果中 根据对上序列的相似度筛选比对结果,优选的候选>90%的基因。9)将7)获 取到的这些序列的属性和来源以及8)获取到的耐药毒力基因数据库比对的结果作为该微 生物基因组的属性和可能来源。
下面为具体的实施例。
实施例
实施例1三代测序样本检测流程构建
1.数据预处理和去宿主处理: 对于mNGS测序1h后的下机样本,首先进行预处理:质控和人源数据库比对进行去人源,得到去人源之后的fastq.gz文件,其内部是非人源的reads。对于纯微生物数据,首先拼接下机进行预处理后的纯三代数据,得到组装数据fasta文件。
2. 序列比对与比对结果筛选:对于微生物基因组所得的序列和去人源的mNGS序 列,先分别通过blast软件和聚类泛基因组数据库中的40个cluster的泛基因组比较。对每 一个比较的基因, 优选的,根据对上序列的覆盖度以>80%为界限。覆盖度计算如 下:
其中为该read或组装序列比对上长度占seq的覆盖度;
为该read或组装序列与泛基因组序列seq比对上的碱基长度;
为该read或组装序列比对上的泛基因组序列seq的碱基长度。
3. 计算相似度:根据2. 中筛选的比对序列结果计算该菌株与某簇的簇相似度,选取簇相似度S最大的3个簇作为该菌株所在的候选簇。簇相似度计 算如下:
其中为fastq或fasta数据中某条测序read或组装序列与该簇泛基因组 某条序列seq比对上的长度占seq的覆盖度;
为fastq或fasta数据中某条测序read或组装序列与该簇泛基 因组某条序列seq的一致性;
为fastq或fasta数据中某条测序read或组装序列与该簇泛基因组某条序 列seq的相似度;
为与该簇泛基因组比对上的第i个测序read或第i个组装序列;
为比对上该簇泛基因组的某条序列seq的read或组装序列的个数;
为泛基因组某序列seq的相似度,即所有比对上该簇泛基因组某序列seq的 测序read或组装序列与该序列的相似度的最大值;
为该簇泛基因组中比上的所有seq的数目。
4. 筛选候选菌株:从3个候选簇中选取株相似度最大的5个作为与样本 相似的候选菌株,即候选簇中取最大的对应的菌株。株相似度计算如 下:
为比对上的簇内泛基因组序列的个数。
5.确定来源菌株并获取来源信息和属性信息:按照株相似度排序,由大 到小,输出聚类泛基因组数据库中汇总的前5株菌的来源,如收集医院、时间等;基因组信 息,如ST-KL型,毒力基因和耐药基因的等基本情况。
6. 对所有微生物基因组和去人源的mNGS三代测序样本,优选的,利用blast软件, 将去人源的mNGS数据或组装数据fasta与聚类泛基因组数据库中的耐药基因和毒力基因数 据库相比较,以相似度判定数据中是否含有该基因;将1)获取的序列read或组装序 列与耐药基因和毒力基因数据库中的每个基因进行序列比对。对比的结果中根据对上序列 的基因相似度筛选比对结果,优选的候选>90%的基因。相似度计算 如下:
其中为fastq或fasta数据中某条测序read或组装序列比对上长度占 耐药基因和毒力基因数据库中基因的覆盖度;
为fastq或fasta数据中某条测序read或组装序列与耐药基因和毒 力基因数据库中基因比对上的碱基长度;
为fastq或fasta数据中某条测序read或组装序列比对上的耐药基因和 毒力基因数据库中基因的碱基长度。
7. 对所有微生物基因组和去人源的宏基因组三代测序样本,在进行2-5步的同时,同时进行第6步。
实施例2 微生物基因组三代测序样本检测
1. 通过三代nanopore测序得到微生物基因组序列Sample1_kpn.fastq.gz,Sample2_kpn.fastq.gz。
将Sample1和Sample2送二代测序,以二代测序的信息作为本发明流程的验证信息。Sample1和Sample2进行二代测序和拼接之后基本的基因组信息,部分药敏信息以及来源见下表1:
表1 样本的基本信息
2. 将Sample1和Sample2利用unicycler软件进行三代测序数据单独拼接,得到拼接Sample1.fasta,Sample2.fasta。
将Sample1.fasta和Sample2.fasta与40个cluster的泛基因组进行比较,得到,按从大到小排序后,前5个的结果见下表2:
表2 部分的情况
3. 以Sample1和Sample2的最大的两个为例,分别与聚类泛基因组数 据库相对比。按从大到小排序后,前5个的结果见下表3:
表3 Cluster1-Sample1部分的情况
表4 Cluster4-Sample2部分的情况
表5 Cluster26-Sample1部分的情况
表6 Cluster29-Sample2部分的情况
注:cluster29内有两株菌。
4. 由3可见Sample1和Sample2在分别最大的两个cluster中,的大小相差很大。比较各个cluster的,将各个cluster的合 并,选取最大的前五个序列作为与待测病原微生物来源相似和基因组属性相似 的菌株输出,此时能通过数据库中菌株初步判定该待测病原微生物各方面的属性以及来 源。
5. 与耐药基因数据库和毒力基因数据库相比较,表7-10为相似度 > 95% 的部分耐药基因和毒力基因展示。由于耐药基因和毒力基因有多个亚型,此时取最 大的亚型为该基因亚型。
表7 Sample1部分耐药基因大于95%的情况
表8 Sample1部分重要毒力基因大于95%的情况
表9 Sample2部分耐药基因大于95%的情况
表10 Sample2部分重要毒力基因大于95%的情况
实施例3 mNGS三代测序样本检测
1. 通过naonpore测序mNGS样本1小时,进行预处理和去人源后,得到Sample3_mNGS.fastq.gz和Sample4_mNGS.fastq.gz。样本的基本信息见下表11:
表11 样本的基本信息
2. 将Sample3和Sample4与40个cluster的泛基因组进行对比,得到,按 从大到小排序后,前5个的结果见下表12:
表12 部分的情况
3. 以Sample3和Sample4的最大的两个所在的cluster为例,与聚类泛 基因组数据库相对比。按从大到小排序后,前5个结果见下表13-16:
表13 Cluster1-Sample3部分的情况
表14 Cluster2-Sample3部分的情况
表15 Cluster3-Sample4部分的情况
表16 Cluster23-Sample4部分的情况
4. 将各个cluster的合并,选取最大的前五个作为与它们自 身的基因组属性以及来源相似的菌株输出,此时能通过数据库中菌株初步判定该菌各方面 的属性以及来源。需注意,Sample4前5个的KL型均不太相同,这是一类无法判定 的情况,原因是mNGS样本对菌株的测序深度不够。我们对Sample4所有测完的数据进行ANI 分析,与它相近的前5个序列有3个是KL19,一个是KL28,一个是KL146。这可能是mNGS建库的 时候菌株浓度不够,无法更加细节的区分KL型。
5. 与耐药基因数据库和毒力基因数据库相比较,表17-19为 > 95%的部 分耐药基因和毒力基因展示。由于耐药基因和毒力基因有多个亚型,此时取最大的 那一亚型为该基因亚型。其中与Sample3比上的毒力基因中大于95%的基因无特殊 需要列出的基因。
表17 Sample3部分耐药基因大于95%的情况
表18 Sample4部分重要耐药基因的情况
表19 Sample4部分重要毒力基因的情况
实施例4数据分析效果比较
在菌株属性分析鉴定领域,通常使用系统进化树分析菌株在进化树上的位置,或使用fastANI计算与其它菌株的相似度,但一次比较的菌株数量过大时会耗费过多时间。本发明方法将从分析用时和准确度方面与直接将序列数据和聚类泛基因组数据库中所有菌相比较。
1.分析鉴定用时分析
将上述2个病原体样本和2个mNGS样本采用本发明方法和fastANI方法和聚类泛基因组数据库比较分析,在同一台服务器并保证相同的CPU下分析。
分析时间如下表20所示,本发明方法用时短于fastANI方法。且随着待测序列数据量的增加,本发明方法的用时优势更明显。
表20 本发明方法和fastANI方法分析时间(s)
2.准确率分析
对2个病原体样本和2个mNGS样本分别采用本发明和fastANI方法基于聚类泛基因组数据库进行分析,统计输出结果。本发明方法在精确度方面与fastANI相当。具体统计结果如下表21所示。
表21 本发明方法和fastANI方法在各个方面的精确度(ST-KL型、耐药和毒力基因存在情况、菌株来源)
前述对本发明的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本发明限定为所公开的精确形式,并且很显然,根据上述教导,可以进行很多改变和变化。对示例性实施例进行选择和描述的目的在于解释本发明的特定原理及其实际应用,从而使得本领域的技术人员能够实现并利用本发明的各种不同的示例性实施方案以及各种不同的选择和改变。本发明的范围意在由权利要求书及其等同形式所限定。

Claims (9)

1.一种通过相似度进行菌株溯源及属性鉴定的方法,其特征在于,包含如下步骤:
a)序列比对:将样本中菌株的测序序列或组装序列和聚类泛基因组数据库中每个簇的泛基因组进行序列比对,其中所述聚类泛基因组数据库中包含泛基因组序列-原始序列-原始菌株关联关系,和属性信息;
b)比对结果筛选:对a)得到的序列比对结果,根据比对上序列的覆盖度筛选比对结果;
c)计算相似度,筛选候选簇:根据b)筛选的比对结果计算样本中菌株与某簇的簇相似度,选取簇相似度最大的L个簇作为样本中菌株所在的候选簇,其中,L为正整数;
d)筛选候选菌株:从每个所述候选簇中选取株相似度最大的M个原始菌株作为与样本中菌株相似的候选菌株,其中,M为正整数,所述株相似度用以衡量样本中菌株与所述聚类泛基因组数据库中关联的原始菌株的相似性;
e)确定来源菌株:从所述候选菌株中,选取株相似度最大的N个原始菌株,作为样本中菌株的来源菌株,其中,N为正整数;
f)获取来源信息及属性信息:根据所述来源菌株获取所述样本中菌株的来源信息及属性信息;
其中,
所述属性信息包括MLST分型、KL分型、收集地点、收集时间中的一个或多个;
步骤c)中,所述簇相似度为比对上某簇泛基因组所有序列的相似度的平均值,所述簇相似度Scluster按以下方法计算:
Sread=Covseq*identityread
Sseq=max0≤i≤n(Sread,i),
其中,Covseq为某条测序read或拼接序列与该簇泛基因组某条序列seq比对上的长度占该seq的覆盖度;
identityread为某条测序read或组装序列与该簇泛基因组某条序列seq的一致性;
Sread为某条测序read或组装序列与该簇泛基因组某条序列seq的相似度;
i为与该簇泛基因组比对上的第i个测序read或第i个组装序列;
n为比对上该簇泛基因组的某条序列seq的read或组装序列的个数;
Sseq为与泛基因组某序列seq的相似度,即所有比对上该簇泛基因组某序列seq的测序read或组装序列与该序列的相似度的最大值;
m为该簇泛基因组中比对上的所有seq的数目;
所述覆盖度为某条测序序列或组装序列与该簇泛基因组某条序列比对上的长度占所述该簇泛基因组某条序列的长度的百分比,
所述属性信息还包括耐药基因和药敏结果。
2.根据权利要求1所述的方法,其特征在于,所述筛选候选簇如下进行:
所述L选取5个至10个,或者
选取簇相似度达到最大簇相似度的99%的簇作为样本中菌株所在的候选簇。
3.根据权利要求1所述的方法,其特征在于,所述株相似度为比对上某簇泛基因组中可以关联到某原始菌株的所有序列相似度之和。
4.根据权利要求1所述的方法,其特征在于,所述聚类泛基因组数据库按以下步骤构建:步骤1)序列相似性计算:计算微生物菌株基因组序列的相似性;
步骤2)序列聚类:根据序列的相似性计算结果进行聚类,将相似序列聚为同一簇;
步骤3)构建聚类泛基因组数据库:对于聚类之后的每一簇,构建该簇所有菌株的泛基因组,作为该簇菌株的序列特征,每簇泛基因组序列包含序列ID以及与原始菌株对应信息,进而构建成物种的聚类泛基因组数据库;
其中,所述的微生物菌株基因组序列包括收集或自建的菌株序列数据及属性信息,和/或收集的公共数据库中菌株序列数据及属性信息。
5.根据权利要求1所述方法,其特征在于,所述测序序列包括质控测序序列,其通过如下步骤获取:
数据预处理:将样本中菌株的测序数据剔除接头、低质量序列和长度过短序列,得到预处理后的数据;
去宿主处理:将获取的所述预处理后的数据与宿主基因组比对,去除比对上宿主基因的序列,从而获得非宿主序列,即为所述质控测序序列;
所述测序序列组装后即获得所述组装序列。
6.根据权利要求1所述的方法,其特征在于,所述样本中菌株的来源信息及属性信息包括耐药基因和/或毒力基因数据库比对结果,所述耐药基因和/或毒力基因数据库比对结果通过如下方法获得:
将样本中菌株的测序序列或组装获取的序列与耐药基因和/或毒力基因数据库进行序列比对,序列比对结果中根据比对上序列的基因相似度进行筛选,判断样本中菌株携带的耐药基因和/或毒力基因。
7.根据权利要求1所述的方法,其特征在于,所述样本中菌株包括一种或多种病原微生物。
8.一种鉴定病原微生物暴发的方法,其特征在于,对特定区域内在不同时间获取的样本中菌株的测序序列或组装序列,按权利要求1至7中任一项所述的方法进行菌株溯源和属性鉴定,从而鉴定病原微生物暴发。
9.一种电子设备,其特征在于,包括:存储器、与所述存储器连接的处理器,及存储在所述存储器上并可在所述处理器上运行的计算机程序,所述处理器运行所述计算机程序时执行权利要求1至8任一项所述的方法。
CN202210991537.3A 2022-08-18 2022-08-18 一种通过相似度进行菌株溯源及属性鉴定的方法 Active CN115064215B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210991537.3A CN115064215B (zh) 2022-08-18 2022-08-18 一种通过相似度进行菌株溯源及属性鉴定的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210991537.3A CN115064215B (zh) 2022-08-18 2022-08-18 一种通过相似度进行菌株溯源及属性鉴定的方法

Publications (2)

Publication Number Publication Date
CN115064215A CN115064215A (zh) 2022-09-16
CN115064215B true CN115064215B (zh) 2023-10-24

Family

ID=83207811

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210991537.3A Active CN115064215B (zh) 2022-08-18 2022-08-18 一种通过相似度进行菌株溯源及属性鉴定的方法

Country Status (1)

Country Link
CN (1) CN115064215B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115527612B (zh) * 2022-10-28 2023-11-14 四川天瓴创新科技集团有限公司 基于数值特征表达的基因组二四代融合组装方法及系统
CN117153248B (zh) * 2023-09-05 2024-05-07 天津极智基因科技有限公司 一种基于泛基因组的基因区变异检测及可视化方法、系统
CN117037912B (zh) * 2023-09-13 2024-06-18 青岛极智医学检验实验室有限公司 一种泛基因组的构建方法、终端设备及存储介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103627800A (zh) * 2013-11-14 2014-03-12 浙江天科高新技术发展有限公司 环境微生物快速检测方法
CN106886689A (zh) * 2015-12-15 2017-06-23 浙江大学 一种病原微生物基因组快速分析方法及系统
CN110004239A (zh) * 2019-04-03 2019-07-12 河海大学 一种基于微生物溯源解析水体污染源的方法
CN111276185A (zh) * 2020-02-18 2020-06-12 上海桑格信息技术有限公司 一种基于二代高通量测序的微生物鉴定分析系统及装置
CN111916151A (zh) * 2020-07-21 2020-11-10 深圳海关动植物检验检疫技术中心 一种苜蓿黄萎病菌的溯源检测方法及应用
WO2021154561A1 (en) * 2020-01-31 2021-08-05 Becton, Dickinson And Company Methods and systems for classifying fluorescent flow cytometer data
CN113744807A (zh) * 2021-11-03 2021-12-03 微岩医学科技(北京)有限公司 一种基于宏基因组学的病原微生物检测方法及装置
CN114420212A (zh) * 2022-01-27 2022-04-29 上海序祯达生物科技有限公司 一种大肠杆菌菌株鉴定方法和系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103627800A (zh) * 2013-11-14 2014-03-12 浙江天科高新技术发展有限公司 环境微生物快速检测方法
CN106886689A (zh) * 2015-12-15 2017-06-23 浙江大学 一种病原微生物基因组快速分析方法及系统
CN110004239A (zh) * 2019-04-03 2019-07-12 河海大学 一种基于微生物溯源解析水体污染源的方法
WO2021154561A1 (en) * 2020-01-31 2021-08-05 Becton, Dickinson And Company Methods and systems for classifying fluorescent flow cytometer data
CN111276185A (zh) * 2020-02-18 2020-06-12 上海桑格信息技术有限公司 一种基于二代高通量测序的微生物鉴定分析系统及装置
CN111916151A (zh) * 2020-07-21 2020-11-10 深圳海关动植物检验检疫技术中心 一种苜蓿黄萎病菌的溯源检测方法及应用
CN113744807A (zh) * 2021-11-03 2021-12-03 微岩医学科技(北京)有限公司 一种基于宏基因组学的病原微生物检测方法及装置
CN114420212A (zh) * 2022-01-27 2022-04-29 上海序祯达生物科技有限公司 一种大肠杆菌菌株鉴定方法和系统

Also Published As

Publication number Publication date
CN115064215A (zh) 2022-09-16

Similar Documents

Publication Publication Date Title
CN115064215B (zh) 一种通过相似度进行菌株溯源及属性鉴定的方法
Bickhart et al. Generating lineage-resolved, complete metagenome-assembled genomes from complex microbial communities
CN111462821B (zh) 病原微生物分析鉴定系统及应用
Pammi et al. Molecular assays for the diagnosis of sepsis in neonates
CN109686439B (zh) 遗传病基因检测的数据分析方法、系统及存储介质
Quicke et al. Utility of the DNA barcoding gene fragment for parasitic wasp phylogeny (Hymenoptera: Ichneumonoidea): data release and new measure of taxonomic congruence
CN111009286A (zh) 对宿主样本进行微生物分析的方法和装置
CN105740650B (zh) 一种快速准确鉴定高通量基因组数据污染源的方法
US20230141128A1 (en) Molecular technology for predicting a phenotypic trait of a bacterium from its genome
CN112687344B (zh) 一种基于宏基因组的人腺病毒分子分型和溯源方法及系统
US20200294628A1 (en) Creation or use of anchor-based data structures for sample-derived characteristic determination
CN1385702A (zh) 提供临床诊断服务的方法
Ames et al. Using populations of human and microbial genomes for organism detection in metagenomes
Sanderson et al. High precision Neisseria gonorrhoeae variant and antimicrobial resistance calling from metagenomic Nanopore sequencing
Kearse et al. The Geneious 6.0. 3 read mapper
CN115312183A (zh) 医学检验报告智能解读方法及系统
CN106951710B (zh) 基于特权信息学习支持向量机的cap数据系统及方法
CN113793647A (zh) 一种基于二代测序宏基因组数据分析装置及方法
CN115083527A (zh) 一种聚类泛基因组数据库构建方法
CN115938491B (zh) 一种用于临床病原诊断的高质量细菌基因组数据库构建方法及系统
CN114496089B (zh) 一种病原微生物鉴定方法
CN111310792A (zh) 一种基于决策树的药敏实验结果识别方法与系统
AU2023261122A1 (en) Construction method for model for analyzing variation detection result
Walter et al. Genomic variant identification methods alter Mycobacterium tuberculosis transmission inference
Czech et al. Scalable methods for post-processing, visualizing, and analyzing phylogenetic placements

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