CN114496085B - 病原微生物宏基因组生信分析参考品及其制备方法和应用 - Google Patents

Abstract

本发明涉及一种病原微生物宏基因组生信分析参考品及其制备方法和应用，属于基因检测技术领域。该方法包括以下步骤：建立丰度分布模型：收集临床样本的宏基因组检测数据，建立自变量为测序序列数目因变量为相对丰度的高斯回归模型；标准化高通量测序数据生成：获取参考基因组序列，模拟生成每种微生物物种预定读长和预定测序错误率的高通量序列数据；Gamma‑泊松分布模型：以Gamma‑泊松分布模型拟合临床样本的宏基因组检测数据；参考品制备：以Gamma‑泊松分布模型随机产生一组模拟样本序列数据，并从标准化高通量测序数据中随机挑选相同数目的测序数据，即得。采用该方法得到的生信分析参考品，可全面地评估生物信息分析流程的灵敏度、特异度、召回率和准确性。

Description

病原微生物宏基因组生信分析参考品及其制备方法和应用

技术领域

本发明涉及基因检测技术领域，特别是涉及一种病原微生物宏基因组生信分析参考品及其制备方法和应用。

背景技术

病原微生物宏基因组技术不依赖于传统的微生物培养，能够快速、无偏差检测临床样本中多种微生物，该项技术解决了临床极大的病原学诊断问题，越来越多的医院将此项技术落地和转化。

但该检测方法由于整体流程复杂，目前仍然缺乏规范化的评价准则和适用的参考品。尤其对于宏基因组学中引起最广泛讨论的生物信息学分析流程，需要一套具有代表性和科学性的生物信息学分析性能确认参考品及其制备的相应技术，仍处于缺失的状态。

为此，国内外已有科学家进行了一些研究，得到了一些制备宏基因组模拟数据的软件和算法。但由于病原宏基因组技术检测出的结果包含临床样本中的致病原、人体微生态定植菌群和试剂工程菌等复杂的成分，仅通过简单地模拟是无法实现真正对病原宏基因组生物信息流程的性能确认。

发明内容

基于此，有必要针对上述问题，提供一种病原微生物宏基因组生信分析参考品的制备方法，采用该方法制备得到的生信分析参考品，可全面地评估生物信息分析流程的灵敏度、特异度、召回率和准确性。

一种病原微生物宏基因组生信分析参考品的制备方法，包括以下步骤：

建立丰度分布模型：收集临床样本的宏基因组检测数据，按照样本类型，针对各微生物物种，分别建立自变量为测序序列数目、因变量为相对丰度的高斯回归模型，得到各微生物物种的丰度高斯分布模型，并获得每种微生物的均值和方差；

标准化高通量测序数据生成：获取上述各微生物物种的参考基因组序列，模拟生成每种微生物物种预定读长和预定测序错误率的高通量序列数据，即为标准化高通量测序数据；

Gamma-泊松分布模型：以Gamma-泊松分布模型方式拟合所述临床样本的宏基因组检测数据，获得模型参数，并建立Gamma-泊松分布模型；

参考品制备：以上述Gamma-泊松分布模型随机产生一组模拟样本序列数据，并根据所述丰度高斯分布模型生成每一种微生物的序列数，再从所述标准化高通量测序数据中随机挑选相同数目的测序数据，将每种微生物挑选出的序列合并，即得所述生信分析参考品。

上述病原微生物宏基因组生信分析参考品的制备方法，首先收集临床宏基因组检测数据（肺泡灌洗液、痰液、血液、脑脊液及其他体液等各种不同类型的生物样本）。利用高斯回归模型统计和建立不同微生物的序列及相对丰度分布模型。同时通过全基因组的进化分析，选择最近源的参考基因组。之后再根据临床样本情况，选择gamma-泊松分布模型进行拟合后建立模型，用于随机挑选微生物参考基因组上产生的高通量测序序列。并根据高斯分布模型和已经建立的微生物序列Gamma-泊松分布模型，随机选择该生物的序列数，从而随机产生已知结果的宏基因组检测数据。该数据可作为生信分析参考品，用于全面地评估生物信息分析流程的灵敏度、特异度、召回率和准确性。

本领域技术人员知晓，上述宏基因组检测数据，也同时表示宏转录组检测数据，其具体是根据病原微生物的具体类型而定。

在其中一个实施例中，所述样本类型为：肺泡灌洗液、痰液、血浆和脑脊液。上述类型的样本，已基本可以涵盖现行宏基因组检测样本类型。

在其中一个实施例中，所述微生物物种包括：致病病原微生物、人体微生态定植菌群和试剂工程菌。可以理解的，本领域技术人员知晓，所述致病病原微生物的具体种类可根据ABX指南和临床微生物手册等所收录的临床致病病原体物种列表选取，也可根据临床共识、论文等进行调整。例如：人体微生态定植菌群包括Propionibacterium_acnes， Comamonas_terrigena，Delftia_acidovorans，Staphylococcus_epidermidis， Veillonella_parvula，Streptococcus_parasanguinis（副溶血链球菌），Streptococcus_ mitis，Streptococcus_oralis，Corynebacterium_striatum，Haemophilus_ parainfluenzae等等；试剂工程菌包括Acinetobacter_guillouiae（桂林不动杆菌）， Acinetobacter_johnsonii（约氏不动杆菌），Sphingomonas_melonis，Comamonas_ terrigena，Moraxella_osloensis，Sphingobium_yanoikuyae，Malassezia_restricta， Pseudomonas_fulva，Pseudomonas_stutzeri，Delftia_tsuruhatensis等等。

在其中一个实施例中，所述预定读长为：50bp，75bp，150bp和300bp；所述标准化高通量测序数据的深度为10X。以上述方式制备得到的参考品，读长上具有广泛地代表性，目前这几种读长的测序策略均为病原宏基因和其他宏基因组学研究常用的测序策略，能够很好地验证现有及未来开发的病原宏基因组生信分析流程；所选择的深度能够很好地保证微生物全基因组都被模拟生成的序列覆盖等优势。

在其中一个实施例中，所述各微生物物种的参考基因组序列通过以下方法选取：计算微生物平均基因核苷酸相似度，计算同一个物种内每一个基因组与其他基因组的相似度，并通过PAM聚类算法计算相似度矩阵的中心点，以所述中心点所属的物种基因组作为参考基因组。

可以理解的，上述相似度的计算可通过本领域常规软件，如FastANI（v1.0）等软件实现。通过上述方式得到的参考基因组序列，更具有代表性。

在其中一个实施例中，所述Gamma-泊松分布模型参数包括：每一个微生物物种分别在肺泡灌洗液、痰液、血浆和脑脊液中的均值、方差及根据Gamma-泊松分布模型拟合的离散度。

在其中一个实施例中，在所述参考品制备步骤之后，还包括指定病原体参考品制备步骤，所述指定病原体参考品制备步骤中，获取指定病原体的参考基因组序列，模拟生成系列梯度数量的高通量测序数据，并将其合并至所述生信分析参考品中，即得指定病原体参考品。

可以理解的，对于系列梯度数量的具体设计，可根据评估要求进行选择，如100,000条、10,000条、1,000条、100条、10条的高通量测序数据等。

在其中一个实施例中，统计每个临床样本出现的最大致病病原微生物种类数量，所述指定病原体的数量上限，不大于单个样本中出现的病原体种类数量均值加3倍标准差。将指定病原体的数量限定在上述范围，根据3σ定理，该范围具有保证所模拟产生的数据所包括的物种数符合99.4%真实世界样本的优势。

本发明公开了上述的病原微生物宏基因组生信分析参考品的制备方法制备得到的病原微生物宏基因组生信分析参考品。

本发明公开了上述的病原微生物宏基因组生信分析参考品在用于评估宏基因组检测生物信息分析流程的灵敏度、特异度、召回率和/或准确性中的应用。

与现有技术相比，本发明具有以下有益效果：

本发明的一种病原微生物宏基因组生信分析参考品的制备方法，通过大样本的数据收集，针对不同样本类型的临床样本、不同类型的微生物占比分布特征（致病原、微生态和试剂工程菌）进行数学建模，最终筛选得到合适的微生物参考基因组进行性能确认参考品数据生成。所得的生信分析参考品可作为临床挑战性实验性能确认标准品使用，可全面地评估生物信息分析流程的灵敏度、特异度、召回率和准确性。

附图说明

图1为肺泡灌洗液样本DNA流程致病病原体、呼吸道微生态和试剂工程菌序列分布；

图2为肺泡灌洗液样本RNA流程致病病原体、呼吸道微生态和试剂工程菌序列分布；

图3为肺泡灌洗液所有物种Gamma-泊松分布模型的重要参数及拟合曲线；

图4为软件Kraken2和Bracken宏基因组物种注释方法的检出性能示意图；

图5为100例肺泡灌洗液模拟数据和10例真实数据β多样性结果展示；

图6为100例肺泡灌洗液模拟数据和10例真实数据Top10物种相对丰度结果展示。

具体实施方式

为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的较佳实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。本文所使用的术语“和／或”包括一个或多个相关的所列项目的任意的和所有的组合。

以下实施例所用方法，如非特别说明，均为常规方法可实现。

实施例1

一种病原微生物宏基因组生信分析参考品，通过以下方法制备得到：

1、建立丰度分布模型。

收集19074例临床样本的宏基因组/转录组检测数据，按照样本类型分为四类，分别为肺泡灌洗液、痰液、血浆和脑脊液。

针对每种样本类型的每种微生物物种，按照常规服从高斯分布随机数的方式分别建立自变量为测序序列数目，因变量为相对丰度的高斯回归模型，得到各微生物物种的丰度高斯分布模型，并获得每种微生物的均值和方差。

上述微生物物种包括致病病原微生物、人体微生态定植菌群和试剂工程菌等，例如，致病病原微生物包括Candida_albicans，Pseudomonas_aeruginosa（铜绿假单胞菌）， Staphylococcus_aureus，Acinetobacter_baumannii，Moraxella_catarrhalis， Klebsiella_pneumoniae（肺炎克雷伯菌），Legionella_pneumophila，Pneumocystis_ jirovecii，Tropheryma_whipplei，Streptococcus_pneumoniae，Human gammaherpesvirus 4（人γ疱疹病毒4），Aspergillus fumigatus，influenza A virus，Sphingomonas_ paucimobilis（少动鞘氨醇单胞菌）等等。人体微生态定植菌群包括Propionibacterium_ acnes，Comamonas_terrigena，Delftia_acidovorans，Staphylococcus_epidermidis（表皮 葡萄球菌），Veillonella_parvula（小韦荣球菌），Streptococcus_parasanguinis， Streptococcus_mitis，Streptococcus_oralis，Corynebacterium_striatum， Haemophilus_parainfluenzae，Rothhia_mucilaginosa（粘液罗氏菌），Prevotella_ melaninogenica（产黑素普雷沃菌），Prevotella_pallens（苍白普雷沃菌）等等；试剂工程菌包括Acinetobacter_guillouiae，Acinetobacter_johnsonii，Sphingomonas_melonis （鞘氨醇单胞菌），Comamonas_terrigena，Meyerozyma_guilliermondii（季也蒙毕赤酵母 菌），Moraxella_osloensis（奥斯陆莫拉菌），Sphingobium_yanoikuyae，Malassezia_ restricta（限制马拉色菌），Pseudomonas_fulva，Pseudomonas_stutzeri，Delftia_ tsuruhatensis，Fusobacterium periodonticum（牙周梭杆菌），Methylobacterium_ radiotolerans（耐辐射甲基杆菌）等等。

2、标准化高通量测序数据生成。

从国家生物技术信息中心(National Center for Biotechnology Information,简称NCBI)上下载上述步骤1中对应的各种临床样本类型不同类型的微生物（常见致病原、微生态和试剂工程菌）的完整基因组序列，通过FastANI（v1.0）软件计算微生物平均基因核苷酸相似度，计算同一个物种内每一个基因组与其他基因组的相似度，并通过PAM聚类算法计算相似度矩阵的中心点，以所述中心点所属的物种基因组作为参考基因组。

再利用单个物种高通量测序数据模拟软件Wgsim，根据参考基因组模拟生成每种微生物不同读长和测序错误率的高通量序列数据，每种物种生成10X深度的数据序列，读长分别为50bp，75bp，150bp和300bp，即得标准化高通量测序数据，备用。

3、Gamma-泊松分布模型。

以Gamma-泊松分布模型方式拟合上述19074例临床样本的宏基因组检测数据，获得Gamma-泊松分布模型参数，包括每一个微生物物种分别在肺泡灌洗液、痰液、血浆和脑脊液中的均值、方差及根据Gamma-泊松分布模型拟合的离散度，并建立Gamma-泊松分布模型。

其中建模方法如下：根据gamma-泊松分布来说，其参数是NB(r,p)的r和p，通过计算条件期望E(Y|X)与条件方差Var(Y|X)，并通过参数重构获得每个物种单独的gamma-泊松分布概率密度函数及其重要参数NB(r,p)，最后使用最大似然估计算法对条件期望E(Y|X)与条件方差Var(Y|X)进行估算。

例如，肺泡灌洗液样本DNA流程致病病原体、呼吸道微生态和试剂工程菌序列gamma-泊松分布模型如图1所示，肺泡灌洗液样本RNA流程致病病原体、呼吸道微生态和试剂工程菌序列gamma-泊松分布模型如图2所示（图中仅展示了部分微生物）。

建模结果如图3所示，图3即为肺泡灌洗液中每个物种的条件期望和条件方差的拟合关系。

4、参考品制备。

第一步，以上述步骤3得到的Gamma-泊松分布模型随机产生一组模拟样本序列数据，即产生每个物种在某一样本类型（如肺泡灌洗液）中的序列数。也即采用Gamma-泊松分布模型拟合、模拟真实环境中临床样本的序列数分布情况。

第二步，由于测序过程是基于随机抽样的过程，其序列的检测为服从高斯分布的的随机过程，因此，依据前述步骤1的丰度高斯分布模型，从前述步骤2已生成的该物种的高通量测序模拟数据（如鲍曼不动杆菌75bp，10X深度）中，随机抽取Gamma-泊松分布模型给出的序列数（即所需抽取的次数），再从所述标准化高通量测序数据中随机挑选相同数目的测序数据，将每种微生物挑选出的序列合并，即得所述生信分析参考品。

实施例2

一种病原微生物宏基因组生信分析指定病原体参考品，参照实施例1的方法制备，其中在步骤3中，统计每个临床样本出现的最大致病病原微生物种类数量。

在步骤4之后，对于输入明确指定的病原体，如指定病原体分别为细菌（Klebsiella pneumoniae，肺炎克雷伯菌）、真菌（Aspergillus fumigatus，烟曲霉）、DNA病毒（Human gammaherpesvirus 4，人γ疱疹病毒4）和RNA病毒（Influenza A virus H1N1，甲型流感病毒），为了更全面分析该病原体检出准确性，模拟了四组数据来分析，根据指定病原体的参考基因组序列，利用wgsim软件模拟生成系列梯度（100,000条、10,000条、1,000条、100条、10条）数量的高通量测序数据，并将这5组数据合并至所述生信分析参考品中，即得指定病原体参考品，可用于指定病原体的干实验性能确认。

所有指定的物种，都按照该原则来操作。而对于指定病原体的数量限制，根据步骤3临床大样本数据统计中，单个样本发现常见病原体数量的均值+3倍标准差（取整）的范围来限定。

实施例3

利用实施例2的生信分析指定病原体参考品评估公开软件Kraken2和Bracken宏基因组物种注释方法的检出性能。

本实施例中，指定病原体分别为细菌（Klebsiella pneumoniae，肺炎克雷伯菌）、真菌（Aspergillus fumigatus，烟曲霉）、DNA病毒（Human gammaherpesvirus 4，人γ疱疹病毒4）和RNA病毒（Influenza A virus H1N1，甲型流感病毒），按照实施例2的方法，得到四组模拟数据作为参考品。

对这四组数据利用Kraken2和Bracken软件进行宏基因组物种注释，根据注释得到的结果和模拟数据中物种真实序列数对软件Kraken2和Bracken宏基因组物种注释方法的检出性能进行评估。

结果如图4所示，图4A为四种病原的注释性能评估情况，图中，横坐标为模拟数据中该物种真实存在的序列数目，纵坐标为Kraken2和Bracken软件注释出的序列数，图中百分比为注释率。

此外，根据是否注释正确，将每条序列归属为真阳性、真阴性、假阳性、假阴性，并计算准确率(Accuracy)、召回率(Recall)以及精确率(Precision)用以评估流程的性能。

结果如图4B所示，可以看到流程精确率基本为100%,准确率和召回率在80%左右，评估指标都非常高。

结合图4A，可以发现此流程注释存在部分假阴性和个别假阳性的结果，这可能是由于Kraken2数据库收纳物种数目较多，由于物种的相似性和序列保守性导致部分序列比对出错导致，整体上来说该流程注释较为准确。

实施例4

本实施例对实施例1得到的生信分析参考品进行验证。

1、方法

使用本算法针对每种样本类型模拟生成100例模拟参考验证数据，每种样本类型设置10例真实临床样本用于验证。

使用Kraken2+Bracken宏基因组生信分析软件进行100例模拟和10例真实临床样本的物种注释分析。

通过计算β多样性和Top10物种相对丰度是否存在显著差异，判定算法的可靠性。

2、结果

结果如图5-6所示，图5为100例肺泡灌洗液模拟数据和10例真实数据β多样性结果展示，图6为100例肺泡灌洗液模拟数据和10例真实数据Top10物种相对丰度结果展示。

结果显示，肺泡灌洗液真实样本和模拟样本的β多样性和Top10物种多样性均无显著差异，体现了本算法的可靠性。

以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种病原微生物宏基因组生信分析参考品的制备方法，其特征在于，包括以下步骤：

建立丰度分布模型：收集临床样本的宏基因组检测数据，按照样本类型，针对各微生物物种，分别建立自变量为测序序列数目、因变量为相对丰度的高斯回归模型，得到各微生物物种的丰度高斯分布模型，并获得每种微生物的均值和方差；

标准化高通量测序数据生成：获取上述各微生物物种的参考基因组序列，模拟生成每种微生物物种预定读长和预定测序错误率的高通量序列数据，即为标准化高通量测序数据；

Gamma-泊松分布模型：以Gamma-泊松分布模型方式拟合所述临床样本的宏基因组检测数据，获得模型参数，并建立Gamma-泊松分布模型；

参考品制备：以上述Gamma-泊松分布模型随机产生一组模拟样本序列数据，并根据所述丰度高斯分布模型生成每一种微生物的序列数，再从所述标准化高通量测序数据中随机挑选相同数目的测序数据，将每种微生物挑选出的序列合并，即得所述生信分析参考品。

2.根据权利要求1所述的病原微生物宏基因组生信分析参考品的制备方法，其特征在于，所述样本类型为：肺泡灌洗液、痰液、血浆和脑脊液。

3.根据权利要求1所述的病原微生物宏基因组生信分析参考品的制备方法，其特征在于，所述微生物物种包括：致病病原微生物、人体微生态定植菌群和试剂工程菌。

4.根据权利要求1所述的病原微生物宏基因组生信分析参考品的制备方法，其特征在于，所述预定读长为：50bp，75bp，150bp或300bp；所述标准化高通量测序数据的深度为10X。

5.根据权利要求1所述的病原微生物宏基因组生信分析参考品的制备方法，其特征在于，所述各微生物物种的参考基因组序列通过以下方法选取：计算微生物平均基因核苷酸相似度，计算同一个物种内每一个基因组与其他基因组的相似度，并通过PAM聚类算法计算相似度矩阵的中心点，以所述中心点所属的物种基因组作为参考基因组。

6.根据权利要求1所述的病原微生物宏基因组生信分析参考品的制备方法，其特征在于，所述Gamma-泊松分布模型参数包括：每一个微生物物种分别在肺泡灌洗液、痰液、血浆和脑脊液中的均值、方差及根据Gamma-泊松分布模型拟合的离散度。

7.根据权利要求1-6任一项所述的病原微生物宏基因组生信分析参考品的制备方法，其特征在于，在所述参考品制备步骤之后，还包括指定病原体参考品制备步骤，所述指定病原体参考品制备步骤中，获取指定病原体的参考基因组序列，模拟生成系列梯度数量的高通量测序数据，并将其合并至所述生信分析参考品中，即得指定病原体参考品。

8.根据权利要求7所述的病原微生物宏基因组生信分析参考品的制备方法，其特征在于，统计每个临床样本出现的最大致病病原微生物种类数量，所述指定病原体的数量上限，不大于单个样本中出现的病原体种类数量均值加3倍标准差。

9.权利要求1-8任一项所述的病原微生物宏基因组生信分析参考品的制备方法制备得到的病原微生物宏基因组生信分析参考品。

10.权利要求9所述的病原微生物宏基因组生信分析参考品在用于评估宏基因组检测生物信息分析流程的灵敏度、特异度、召回率和/或准确性中的应用。

Images