CN109273053A - 一种高通量测序的微生物数据处理方法 - Google Patents

Abstract

本发明公开了一种高通量测序的微生物数据处理方法，其中，所述方法包括：高通量测序的微生物16sRNA读段进行重叠群组装、分箱，以q‑PCR标记微生物重叠群，使所述微生物重叠群包含标记基因，去除含有标记基因的生物重叠群，获得高质量微生物宏基因组测序数据。本发明通过序列聚类等方法鉴定的去除来源于污染物的序列，得到更为高纯度的微生物宏基因组测序数据，保证基于微生物宏转录组测序数据的基因表达结果更为准确。本发明以微生物宏基因组测序数据作为研究对象，基于生物信息学思路，提高微生物宏基因组测序数据的质量。

Description

一种高通量测序的微生物数据处理方法

技术领域

本发明涉及一种高通量测序的微生物数据处理方法，属于高通量测序质量控制领域。

背景技术

下一代测序技术(NGS)又称高通量测序，以高输出量和高解析度为主要特色，能一次并行对几十万到几百万条DNA分子进行序列读取，在提供丰富的遗传学信息的同时，还可以大大降低测序费用，缩短测序时间的测序技术。由于高通量测序技术数据处理量大，处理内容繁杂，因此对于测序质量的控制、污染源的确定与排除成为了一个重要的研究课题。测序质量的影响因素是多方面的，常见的影响因素多数为操作中的误差，跨越日期和组处理数据的批次效应的主要来源已被确定为实验性的，如来自DNA提取试剂盒，PCR批次或测序仪器，而不是生物学。这一问题在“多物种”NGS数据处理中尤为突出，一旦污染和测序，将读数与目标和污染物分开并不是一项简单的任务，即使污染物可以轻松识别。在大多数情况下，目标和污染都没有完整的基因组，这使得分配读数非常具有挑战性。因此迫切需要利用目标物种的有限信息去除环境微生物污染物。

目前已经研究并提出了一些基于相似性或组成信息的分箱程序。然而，互补的上游和下游加工方法需要与重叠群装箱结合以达到更高的灵敏度和特异性。一种有前途的污染物鉴定和过滤解决方案是宏基因组方法，它促进了污染微生物基因组的分类学和功能分析。已经针对基于不同计算方法的可能污染物的分析进行了一些软件研发：SourceTracker，其应用贝叶斯推断方法来估计微生物污染的组成和丰度；DeconSeq，可以通过长读取对齐处理可能来自人类的污染，目前采用上述两种软件并基于重叠群聚类方法，可以成功地区分读数与目标物种和污染物。然而，读取分配的假阳性率仍然很高，并且没有考虑潜在有价值的信息，例如多个样本(具有相似污染物)中某些目标物种的丰度相关性。因此，迫切需要对当前流程进行高级优化。

发明内容

针对现有技术存在的上述问题，本发明的目的是获得一种高通量测序的微生物数据处理方法。

为实现上述发明目的，本发明采用的高通量测序的微生物数据处理方法的技术方案如下：

所述方法包括：高通量测序的微生物16sRNA读段进行重叠群组装、分箱，以q-PCR标记微生物重叠群，使所述微生物重叠群包含标记基因，去除含有标记基因的生物重叠群，获得高质量微生物宏基因组测序数据。

优选的，所述数据处理方法包括如下步骤：

a)通过已公开序列建立模拟数据集，通过高通量测序数据建立真实宏基因数据集；

b)对数据集内数据进行质量控制，去掉低质量的碱基和读段，提取16sRNA基因；

c)采用Parallel-Meta pipeline(version 2.0)软件生成分类学概况，通过HMM从原始测序数据中提取16s rRNA序列，以Greengene数据库搜索序列以确定物种的来源，16sRNA数据搜索物种数量；

d)采用VELVET、MEGAHIT软件进行重叠群组装，使用velvet命令从velveth获得的k-mers构建de Bruijn图并提取重叠群，使用MEGAHIT对de Bruijn图进行汇编；

e)以q-PCR方法标记基因，将含有标记基因的重叠群鉴定为目标物质，对目标物种聚类鉴定，获得去除污染的高质量微生物宏基因组测序数据。

优选的，模拟数据集(simulated)收集现有数据库的目标和污染物种的高通量读段模拟真实情况。模拟数据集的目标和污染物种的reads来自NeSSM。在这项研究中，我们假设每个样本中只有一种目标细菌种类，本研究中使用的目标细菌种类包括几种模式生物：如枯草芽孢杆菌，金黄色葡萄球菌，大肠杆菌。

优选的，真实宏基因数据集包括测序结果集(ab initio)和体内对照集(invivo)。对于ab initio数据集，将污染物种的真实测序数据与来自样本的真实宏基因组序列混合，混合不同的浓度梯度。

优选的，真实宏基因数据集包括高通量测序数据集和体内对照集。

其中，步骤b)具体为：

对测序数据进行质量控制，通过Trimmomatic-0.36修剪实际测序数据以去掉低质量的碱基和读段，切断3个前导/拖尾碱基，丢弃长度小于100bp的读段。所述低质量的碱基为Qvalue_average<35或Qvalue<20的碱基。

优选的，步骤c)具体为：为了增加序列的特异性，应用两个组装器来组装来自总读段中的重叠群。选择的组装工具之一是Velvet，它可以处理从头基因组装配和短序列读取比对。对于Velvet，使用velveth命令构建数据集作为准备工作，并使用velvet命令从velveth获得的k-mers构建de Bruijn图并提取重叠群。对于所有分析，k设置为12。其他参数设置为自动或默认。另一个是MEGAHIT，这是最近通过简洁的de Bruijn图设计的复杂宏基因组学的汇编程序。值得一提的是，通过使用这两个工具，本质上都考虑了丰度信息。

对于模拟的宏基因组，在两个装配工上进行装配以比较它们的性能。提取并比较基本装配统计数据。由于MEGAHIT通过模拟数据分析显示优于Velvet，因此仅使用MEGAHIT处理ab initio和体内数据集。

优选的，模拟数据库的污染物种选自：酿酒酵母、枯草芽孢杆菌、金黄色葡萄球菌和/或大肠杆菌。

优选的，当采用酿酒酵母处理数据时，添加真核生物18S rRNA数据库。

优选的，枯草芽孢杆菌的标记包括：ftsZ、lytF、nsrR、spo0A、ygxB、yjbH和/或yjbI。

优选的，大肠杆菌的标记包括：acpP、casA、cof、dxs、fabB、fabF、leuO、tesA和/或uidA。

本发明是一种基于微生物组学和生物信息学的微生物宏基因组测序数据的数据处理方法。本发明获得待分析的原始微生物宏基因组测序数据；通过生物信息学方法分析可能的污染物定性和定量；通过序列聚类等方法鉴定的去除来源于污染物的序列，得到更为高纯度的微生物宏基因组测序数据，保证基于微生物宏转录组测序数据的基因表达结果更为准确。本发明以微生物宏基因组测序数据作为研究对象，基于生物信息学思路，提高微生物宏转录组测序数据的质量。

与现有技术相比，本发明以微生物宏基因组测序数据为研究对象，基于生物信息学思路，利用比对组装聚类等技术对其进行研究，相比传统方法而言，该方法有以下优点：

(1)聚类纯度高，重叠群分箱方法可以解决来自多个物种的单个高浓度和纯净的目标簇；

(2)灵敏度和特异性好，目标物种的灵敏度和特异性值读取MEGAHIT处理数据的分配都较高。测试发现四个样本的敏感度值平均为92.7％，而目标重叠群和读数的敏感性值均更高，可以达到100％；

(3)覆盖率高，读段和重叠群分布集中；

(4)效率高，QC-Blind在具有4,000个配对末端读数的数据集上的运行时间在单个处理器中小于12小时。

附图说明

图1为本发明的方法流程图。

具体实施方式

下面结合实施例对本发明提供的一种高通量测序的微生物数据处理方法作进一步详细、完整地说明。下面描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

下述实施例中的实验方法，如无特殊说明，均为常规方法。下述实施例中所用的实验材料如无特殊说明，均为市场购买得到。

本实施例以人类样本唾液内的微生物高通量测序结果为例进行质量控制，具体操作步骤如下：

一、模拟和真实数据集

1.有关模拟和真实宏基因组数据集的信息。

本实施例选择了三种类型的宏基因数据集：simulated(模拟数据集)、ab initio(测序结果集)和in vivo(体内对照集)；数据集具体信息如下表1所示。

表1

表1中，“Simu_BS_5％”表示枯草芽孢杆菌是目标物种，并且来自该目标物种的读数构成该样品中所有读数的5％，其他条目内容类似。对于ab initio数据集，样本名称的定义类似。所有物种的参考基因组均从NCBI Microbial Genomes网站下载。

2.宏基因组数据准备

对于计算机模拟数据集，目标和污染物种的读段(reads)来自NeSSM。在这项研究中，我们假设每个样本中只有一种目标细菌种类，本研究中使用的目标细菌种类包括几种模式生物：枯草芽孢杆菌，金黄色葡萄球菌，大肠杆菌。将它们的读段与人口腔微生物群落中5或10个代表性物种的基因组产生的读段混合(称为HOB(5/10))，其被用作可能的人类污染物。来自目标物种的读段的梯度比例设定为5％，35％，65％，95％。此外，我们将酿酒酵母与枯草芽孢杆菌和10种口腔细菌结合起来，模拟真核污染的特殊情况(数据集D)。

3.ab initio数据准备

对于ab initio数据集，我们将枯草芽孢杆菌的真实测序数据与来自人类唾液样本的真实宏基因组序列混合(表1)，其中来自目标物种(枯草芽孢杆菌)的读数的相对比例设定为35％，65％，95％对于不同的数据集。样品命名为AB_BS_35％，AB_BS_65％，AB_BS_95％。

4.in vivo数据准备

在该研究中使用的体内数据集是来自如下制备的真实群落样品的宏基因组(非16srRNA)数据集：在活化后，将枯草芽孢杆菌168培养过夜直至OD600值在0.6和0.8之间。将所有枯草芽孢杆菌以12000转/分钟(12114g)离心进行以下实验。从三个健康成人中收集新鲜唾液，不含饮用水或在采集样品前约30分钟漱口。然后在离心前将200ul新鲜唾液加入枯草芽孢杆菌培养物中。此示例已命名为Real_BS。

二、数据分析

质量控制过程如图1所示。首先，通过Trimmomatic-0.36修剪实际测序数据以去除低质量(如Qvalue_average<35或Qvalue<20)碱基和读数。如果质量得分低于质量阈值，则切断3个前导/拖尾碱基。长度太短(默认值<100bp)的读段丢弃掉。然后提取16s rRNA基因用于物种鉴定和定量。最后，按顺序进行读段收集、重叠群分箱和目标基因对比。

1.识别目标物种和污染物种

分类学概况由Parallel-Meta pipeline(version2.0)生成。通过HMM从原始测序数据中提取16s rRNA序列，并且针对Greengene数据库搜索这些序列以确定物种的来源。产生的物种总数作为重叠群分箱步骤的输入，为聚类提供了更好的准确性。另外，当用酿酒酵母(Saccharomyces cerevisiae)处理数据集时，添加真核生物18S rRNA数据库以供参考。对于未在16s RNA或18S rRNA数据库中记录信息的未知污染物，选择聚类簇的数目将有困难，但这种方法足够现实，可以很好地用于目标识别和污染过滤。

2.重叠群组装

为了增加序列的特异性，应用两个组装器组装来自community读段的重叠群。选择的组装工具之一是Velvet，它可以处理从头基因组装配和短序列读取比对。对于Velvet，我们使用velveth命令构建数据集作为准备工作，并使用velvet命令从velveth获得的k-mers构建de Bruijn图并提取重叠群。对于所有分析，k设置为12。其他参数设置为自动或默认。另一个是MEGAHIT，这是最近通过简洁的de Bruijn图设计的复杂宏基因组学的汇编程序。值得一提的是，通过使用这两个工具，本质上都考虑了丰度信息。

对于模拟的宏基因组，在两个装配工上进行装配以比较它们的性能。提取并比较基本装配统计数据。由于MEGAHIT通过模拟数据分析显示优于Velvet，因此仅使用MEGAHIT处理ab initio和in vivo数据集。

3.重叠群分箱

重叠群(conting)分箱是该方法中最重要的步骤之一。在所有现有的分箱算法中，选择了CONCOCT，因为CONCOCT在重叠群组装中考虑了多个样品的序列组成和覆盖率，并且它可以处理单个样品和多个样品，这可能适合用于批量处理受污染的样品。在CONCOCT处理后，重叠群将聚集成物种水平的群体。

4.标记基因选择和比对

利用标记基因进行目标物种聚类鉴定适用于没有完整或部分参考基因组但只有少数标记基因的目标物种。基因越独特，识别就越灵敏和特异。通过咨询文献，选择q-PCR标记ftsZ，lytF，nsrR，spo0A，ygxB，yjbH，yjbI用于枯草芽孢杆菌，acpP，casA，cof，dxs，fabB，fabF，leuO，tesA，uidA用于大肠杆菌。

然后，将含有目标物种的标记基因的重叠群鉴定为属于目标物种(定义为targetcontigs)。基于这些分配，通过BOWTIE2将原始读段映射到target contigs。然后可以评估每个步骤中的总读段和目标读段的统计数据。对于ab initio和in vivo数据集，只有目标读段或重叠群分类，因为将每种污染物读段分类到其来源物种是不切实际的，特别是其中许多还没有对其全基因组测序。

三、结果评估

3.1在读段组装和重叠群聚类方面

这种重叠群分箱方法可以解决来自多个物种的单个高浓度和纯净的目标簇。考虑到在模拟基因组上的读取映射期间产生的可能的伪像，可以合理预期该方法实际上对于真实数据集实际上更好。

3.2在灵敏度特异性方面

一般而言，目标物种的灵敏度和特异性值读取MEGAHIT处理数据的分配都较高。四个样本的敏感度值平均为92.7％，而目标重叠群和读数的敏感性值均更高：Simu_BS_5％，Simu_BS_35％和Simu_BS_65％的100％分配特异性，表明目标群中的目标信息可以是保留很少的污染物。然而，Velvet处理数据的灵敏度和特异性评估在具有5％目标读数的数据集中极低(34.3％，与MEGAHIT中的93.5％相比)，这引起了对Velvet处理严重污染数据的能力的质疑。在重叠群水平的Velvet敏感性也不乐观。结合起来，对目标物种读取分配的敏感性和特异性的评估显示了在QC-Blind方法中使用MEGAHIT的优越性。因此，在以下分析中，在QC-Blind方法中采用MEGAHIT作为默认值。

3.3在效率方面

通常，QC-Blind在具有4,000个配对末端读数的数据集上的运行时间在单个处理器中小于12小时，这取决于测序质量和重叠群数量。最大比例的时间用于重叠群装箱。由于聚类算法的时间复杂性，由于重叠群长度的截止值较低，运行时间明显增加，而且重叠群数量增加。因此，应该考虑合理的截止，改进的聚类算法以及确定多个处理器的利用率以达到更好的效率。

最后有必要在此说明的是：以上实施例只用于对本发明的技术方案作进一步详细地说明，不能理解为对本发明保护范围的限制，本领域的技术人员根据本发明的上述内容作出的一些非本质的改进和调整均属于本发明的保护范围。

Claims (10)

1.一种高通量测序的微生物数据处理方法，其特征在于，所述方法包括：高通量测序的微生物16sRNA读段进行重叠群组装、分箱，以q-PCR标记微生物重叠群，使所述微生物重叠群包含标记基因，去除含有标记基因的生物重叠群，获得高质量微生物宏基因组测序数据。

2.根据权利要求1所述的高通量测序的微生物数据处理方法，其特征在于，所述数据处理方法包括如下步骤：

a)通过已公开序列建立模拟数据集，通过高通量测序数据建立真实宏基因组数据集；

b)对数据集内数据进行质量控制，去掉低质量的碱基和读段，提取16s rRNA基因；

c)生成分类学概况，通过HMM从原始测序数据中提取16s rRNA序列，以Greengene数据库搜索序列以确定物种的来源，16sRNA数据搜索物种数量；

d)进行重叠群组装，使用velvet命令从velveth获得的k-mers构建de Bruijn图并提取重叠群，使用MEGAHIT对de Bruijn图进行汇编；

e)以q-PCR方法标记基因，将含有标记基因的重叠群鉴定为目标物质，对目标物种聚类鉴定，获得去除污染的高质量微生物宏基因组测序数据。

3.根据权利要求1所述的高通量测序的微生物数据处理方法，其特征在于，模拟数据集收集现有数据库的目标和污染物种的高通量读段模拟真实情况。

4.根据权利要求1所述的高通量测序的微生物数据处理方法，其特征在于，真实宏基因数据集包括高通量测序数据集和体内对照集。

5.根据权利要求1所述的高通量测序的微生物数据处理方法，其特征在于，步骤b)具体为：

对测序数据进行质量控制，通过Trimmomatic-0.36修剪实际测序数据以去掉低质量的碱基和读段，切断3个前导/拖尾碱基，丢弃长度小于100bp的读段。

6.根据权利要求1所述的高通量测序的微生物数据处理方法，其特征在于，模拟数据库的污染物种选自：酿酒酵母、枯草芽孢杆菌、金黄色葡萄球菌和/或大肠杆菌。

7.根据权利要求6所述的高通量测序的微生物数据处理方法，其特征在于，当采用酿酒酵母处理数据时，添加真核生物18S rRNA数据库。

8.根据权利要求6所述的高通量测序的微生物数据处理方法，其特征在于，枯草芽孢杆菌的标记包括：ftsZ、lytF、nsrR、spo0A、ygxB、yjbH和/或yjbI。

9.根据权利要求6所述的高通量测序的微生物数据处理方法，其特征在于，大肠杆菌的标记包括：acpP、casA、cof、dxs、fabB、fabF、leuO、tesA和/或uidA。

10.根据权利要求1所述的高通量测序的微生物数据处理方法，其特征在于：真实宏基因数据集包括测序结果集和体内对照集。