一种微生物基因组数据库构建方法及其应用与流程

本发明涉及生物信息学领域，特别是涉及一种微生物基因组数据库构建方法及其在微生物鉴定中的应用。

背景技术：

宏基因组测序(metagenomicsnextgenerationsequencing,mngs)不依赖于传统的微生物培养，直接对临床样本中的核酸进行高通量测序，能够快速、客观的检测临床样本中的多种病原微生物(包括病毒、细菌、真菌、寄生虫)。随着mngs技术平台的完善和临床研究的增多，mngs在临床上的运用也越来越广泛。宏基因组测序分析中两个重要的部分，一部分是微生物基因组数据库的构建，一部分是比对结果的分析筛选，而微生物基因组数据库的构建方式又影响和决定着比对结果筛选的准确度和特异度。

微生物基因组数据库即收录了多种微生物基因组，一般都是通过ncbi公共数据库的基因组数据的下载筛选加工完成。在公共数据库中，每个物种可能有多个不同的菌株基因组，一般会选择其中一株作为该物种的代表基因组，但由于微生物菌株间基因组的突变率较高，而宏基因组测序的覆盖率较低，这时会造成假阴性现象。而如果将全部菌株纳入到数据库中，则会使数据库的量变的很大，需要更多的计算和存储资源，且会延长分析时间。此外，相同科下各属内微生物各物种的基因组序列之间，存在很多共有序列。近缘物种间共有序列比例则会更高，在仅有某一物种出现时，由于共有序列的比对，可能会导致判断同时出现另一物种的假阳性现象。

因此，构建一种既可以覆盖物种内大部分菌株序列，又可以覆盖物种内或物种间共有序列的数据库是一种优选方案。

鉴于此，提出本发明。

技术实现要素：

本发明的目的是寻求一种既可以覆盖物种内大部分菌株序列，又可以覆盖物种内或物种间共有序列的数据库及其构建方法。本发明同时寻求通过对reads比对结果的分析，获得准确的物种鉴定结果。

为实现上述目的，本发明提出如下技术方案：

本发明首先提供了一种微生物基因组数据库构建方法，包括如下步骤：

1)数据获取：获取不同微生物物种的基因组数据；

2)共有序列集和特异序列集鉴定：将不同物种基因组序列分别打断，得到各物种短片段序列集；将所有物种基因组序列合并，一起进行基因组打断，得到所有物种的短片段序列集；将所有物种短片段序列集与各物种短片段序列集进行序列比较，仅在一个物种短片段序列集中出现的序列为该一个物种的特异序列，出现在大于一个物种序列集的序列为多物种共有序列；

3)数据集合并制备数据库：将上述多物种共有序列集和各物种特异序列集分别标签化处理，并汇总得到数据库。

在一些实施方式中，步骤1)中还包括质粒序列去除。

在一些实施方式中，步骤2)中所述打断是通过切kmer方式打断；

在一些实施方式中，所述打断为长度n移步步长为k的打断方式，其中n取值其中n取值为76～1000，k取值为1～n的自然数，优选的n取值为76～300，k取值为1。

进一步的，所述方法还包括如下步骤：

4)构建比对分值稀疏矩阵：对其中一个物种的特异序列与其它任一物种的特异序列进行两两比对，获得两两物种序列间的比对分值，存储为稀疏矩阵。

在一些具体的实施方式中，所述构建比对分值稀疏矩阵的方法为：使用makeblastdb对其中一个物种的特异序列进行建库，使用其它任一物种的特异序列与此数据库进行比对，将两两物种序列编号，以及两物种间的bitscore值存储为稀疏矩阵。由此，构建存储两两物种间的特异序列，以及特异序列间bitscore值的多个稀疏矩阵。

本发明还提供一种微生物基因组数据库，其特征在于：由上述方法制备。

本发明还提供一种微生物基因组数据库，所述数据库包含不同物种基因序列，其中所述不同物种基因序列为打断后标签化的物种共有序列和标签化的物种特异序列；所述标签化为数据库中存储序列来源信息的字段或文件；

在一些优选的实施方式中，所述打断为长度n移步步长为k的打断方式，更优选的使用jellyfish软件进行打断。

在一些优选的实施方式中，所述打断后还包括去冗余步骤。

进一步的，所述数据库还包括比对分值稀疏矩阵，所述比对分值稀疏矩阵为存储物种特异序列间比对分值的矩阵或表格。

本发明还提供一种微生物鉴定方法，包括如下步骤：

1)数据比对：将测序下机reads比对至上述微生物基因组数据库，得到初步比对结果；

2)比对结果筛选:在各reads的比对结果中，以最大比对分值为基准，将比对分值低于最大比对分值的90-95％的比对结果去除；优选的，将比对分值低于最大比对分值的95％的比对结果去除。

3)比对结果统计。

进一步的，所述比对结果统计步骤为：每条reads的比对结果有如下几种可能，分别按照如下原则，对reads所属物种进行鉴定，

本发明还提供一种上述数据库在微生物物种鉴定中的应用本发明有益的技术效果：

1)本发明通过区分物种内以及物种间共有和特异序列，并加标签方式构建数据库，可以预先获得物种的基因组序列间的相似或特异区段，为后续鉴定提供参考信息。

2)本发明对于比对到两物种的共有序列的reads，按照唯一比对到特异序列的reads的比例进行拆分，相比以往鉴定方法，如每个reads按照1/2拆分给各个物种，出现假阳性的可能性更低。

3)对于物种基因组与数据库中序列有差异的物种有较好的兼容性。本发明使用物种特异序列间的bitscore值作为特异序列的物种区分指标，如果出现了超出此数据库判断范围的序列，则会进行拆分处理。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为实施例2和3中各方法分析模拟测序数据集得到的大肠杆菌reads条数结果。

图2为实施例2和3中各方法分析模拟测序数据集得到的志贺氏菌reads条数结果。

具体实施方式

下面将结合实施例对本发明的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围，并且所述实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

部分术语定义

除非在下文中另有定义，本发明具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本发明。

如本发明中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由…组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。

本发明中的术语“大约”、“大体”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10％，优选±5％。

在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。

此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本发明描述的实施方案能以不同于本发明描述或举例说明的其它顺序实施。

本发明的种微生物基因组数据库构建方法，包括如下步骤：

1)数据获取：获取不同微生物物种的基因组数据；2)共有序列集和特异序列集鉴定：将不同物种基因组序列分别打断，得到各物种短片段序列集；将所有物种基因组序列合并，一起进行基因组打断，得到所有物种的短片段序列集；将所有物种短片段序列集与各物种短片段序列集进行序列比较，仅在一个物种短片段序列集中出现的序列为该一个物种的特异序列，出现在大于一个物种序列集的序列为多物种共有序列；3)数据集合并制备数据库：将上述多物种共有序列集和各物种特异序列集分别标签化处理，并汇总得到数据库。优选的，为避免质粒序列对鉴定影响，步骤1)中还包括质粒序列去除，去除上述基因组中存在的质粒序列，得到去质粒后基因组序列。可以理解，上述去除质粒序列步骤中，可根据具体情况调整，对于有质粒序列的基因组，进行去除质粒序列处理，没有质粒的基因组则不需要去除。步骤2)中所述打断是将所述菌株基因组打断为长度为n且移步步长为k的序列集；例如，可截取基因组序列第1个碱基至第n个碱基之间的序列为第一条序列，然后起始位置定位至第k个碱基，截取第k至k+n位碱基间的序列为第二条序列，依次类推，对于截取的第n条序列，其位置为(n-1)k至(n-1)k+n。打断可以是多种方式，优选的通过切kmer方式打断；根据实际需要，可以设置不同的打断长度和步长，在一些实施方式中，所述打断为长度n移步步长为k的打断方式，其中n取值为76～1000，k取值为1～n的自然数，优选的n取值为76～300，k取值为1。

进一步的，所述方法还包括如下步骤：4)构建比对分值稀疏矩阵：对其中一个物种的特异序列与其它任一物种的特异序列进行两两比对，获得两两物种序列间的比对分值，存储为稀疏矩阵。优选的，所述构建比对分值稀疏矩阵的方法为：使用makeblastdb对其中一个物种的特异序列进行建库，使用其它任一物种的特异序列与此数据库进行比对，将两两物种序列编号，以及两物种间的bitscore值存储为稀疏矩阵。由此，构建存储两两物种间的特异序列，以及特异序列间bitscore值的多个稀疏矩阵。

本发明的微生物基因组数据库，包含不同物种基因序列，其中所述不同物种基因序列为打断后标签化的物种共有序列和标签化的物种特异序列；所述标签化为数据库中存储序列来源信息的字段或文件；在一些优选的实施方式中，所述打断为长度n移步步长为k的打断方式；此外，所述数据库还包括比对分值稀疏矩阵，所述比对分值稀疏矩阵为存储物种特异序列间比对分值的矩阵或表格。

本发明的微生物鉴定方法，包括如下步骤：

1)数据比对：将测序下机reads比对至上述微生物基因组数据库，得到初步比对结果；2)比对结果筛选:在各reads的比对结果中，以最大比对分值为基准，将比对分值低于最大比对分值的90-95％(优选95％)的比对结果去除；3)比对结果统计，所述比对结果统计步骤为：每条reads的比对结果有如下几种可能，分别按照如下原则，对reads所属物种进行鉴定，

以下描述仅是为了帮助理解本发明而提供。这些描述不应被理解为具有小于本领域技术人员所理解的范围。

实施例1方法建立

一、微生物基因组数据库构建方法：

1)数据获取：获取微生物物种的代表基因组数据，每一个物种的每一个菌株可能有多个基因组序列，比如当从ncbi中获得基因组数据时，优先筛选出refseqcategory标记为“referencegenome”和“representativegenome”的基因组序列作为该物种该菌株的基因组序列；如果没有“referencegenome”或“representativegenome”的基因组，则选择标记为“na”的基因组作为基因组序列。

2)质粒序列去除：为避免质粒序列对鉴定影响，去除上述基因组中存在的质粒序列，得到去质粒后基因组序列。

3)鉴定共有序列集和特异序列集：分别将各物种的上述去质粒后基因组打断，得到各物种的短片段序列集；将所有物种的上述去质粒基因组进行合并，然后进行基因组打断，得到所有物种的短片段序列集；将所有物种短片段序列集与各物种的短片段序列集进行序列比较，仅在一个物种短片段序列集中出现的序列为此物种特异序列，出现在大于一个物种序列集的序列为多物种共有序列。

4)数据集合并，构建数据库：将上述多物种共有序列集和各物种特异序列集分别标签化处理，并汇总得到数据库。

进一步的：

5)构建比对分值稀疏矩阵：使用makeblastdb对其中一个物种的特异序列进行建库，使用其它任一物种的特异序列与此数据库进行比对，将两两物种序列编号，以及两物种间的bitscore值存储为稀疏矩阵。由此，构建存储两两物种间的特异序列，以及特异序列间bitscore值的多个稀疏矩阵。

通过上述方法，即获得了微生物基因组数据库。

二、微生物鉴定方法：

1)测序reads数据获取以及过滤：获取高通量测序得到的测序reads数据，并依次进行低质量序列过滤，宿主序列过滤。

2)数据比对：将上述过滤后的reads比对至构建的微生物基因组数据库，得到初步比对结果，比对软件优选blastn。

3)比对结果筛选:在各reads的比对结果中，以最大bitscore值为基准，将比对分值低于最大比对分值的95％的比对结果去除。

4)比对结果统计：每条reads的比对结果有如下几种可能，分别按照如下原则，对reads所属物种进行鉴定：

按照上述方案将reads结果进行整理汇总，即为测序数据的鉴定结果。

实施例2大肠杆菌和志贺氏菌数据构建

以下通过大肠杆菌(escherichiacoli)和志贺氏菌(shigellaflexneri)物种为例，进行数据库构建。

1数据获取：

从ncbi上下载微生物基因组序列，获得大肠杆菌2个菌株的基因组序列gcf_000008865.2和gcf_003697165.2，以及志贺氏菌2个菌株的基因组序列，gcf_000006925.2和gcf_007197595.1。

2去除质粒序列：根据基因组序列文件里面的序列名称，将带有plasmid(plasmid)的序列去除。

3鉴定共有和特异序列集：将大肠杆菌的2个基因组序列合并，然后使用jellyfish按照长度76bp，步长1bp打断，jellyfish包含去冗余过程，得到序列集合1；将志贺氏菌的2个基因组序列合并，然后使用jellyfish按照长度76bp，步长1bp打断，得到集合2；将4个物种的基因组合并，然后使用jellyfish按照长度76bp，步长1bp打断，得到集合3；将序列集合1，集合2，集合3进行比较，仅出现在集合1中，未出现在集合2中的序列为大肠杆菌特异序列；仅出现在集合2中，未出现在集合1中的序列为志贺氏菌特异序列；在序列集合1和2中均出现的序列为共有序列。

获得共有序列2660763条，获得大肠杆菌特异序列6459145条，获得志贺氏菌特异序列4356643条。

4构建数据库：使用blast软件命令makeblastdb对序列集合3进行blast比对库构建。同时向序列集合3的每个序列的序列名称中添加新的字段，该字段根据序列来源分别填充为共有序列，大肠杆菌特异，志贺氏菌特异。

5构建比对分值稀疏矩阵：将序列集合1进行blast比对库构建，将序列集合2中的序列与此库进行比对，获得两数据库中两两序列的bitscore值，将两两序列的的编号以及bitscore值构建稀疏矩阵进行存储。

6对构建数据库的4个基因组序列，抽取模拟数据进行测试。

每个基因组序列随机抽取1000条76bp长度的序列模拟为测序下机reads，获得4个模拟数据；

按照一定比例组合抽取不同物种的1000条76bp长度的序列，模拟为测序下机reads，生成2个模拟数据；

对gcf_000008865.2去质粒后序列，按照变异程度5％的比例进行碱基的随机替换，然后抽取1000条76bp长度的序列，模拟为测序下机reads，作为变异菌株模拟数据。

7对模拟数据使用blastn软件，与构建好的数据库进行比对。

8比对结果筛选:各read的比对结果中，获得比对结果中最大的比对分值，去除本read其它比对结果比对分值低于最大比对分值95％的结果。

9比对结果统计，按照物种分配原则，对比对结果中各reads的物种归属进行统计：

如果此reads的比对到各物种特异序列的bitscore值，均高于两物种的特异序列间相互比对的bitscore值，则此reads分配给两个物种，各1/2，此处理适用的情况举例如下：序列a为大肠杆菌的特异序列，序列b为志贺氏菌的特序列，假定a与b比对的bitscore值为126；此时测序获得序列c，序列c与a的比对的bitscore值为130，c与b比对的bitscore值为135，即均高于126，则可以理解为目前构建的物种的数据库对c序列不具有物种区分度，c序列相对于a序列与b序列更相似，c序列相对于b序列与a序列更相似，则c序列判定为各物种1/2。

共有序列的物种划分：计算reads结果中，唯一比对到各物种特异序列的比例，以此比例拆分分配为共有序列的reads，即为各物种reads。此处理适用的情况如下：1000条reads中，比对结果中唯一比对到大肠杆菌的是100条，唯一比对到志贺氏菌的是200条，比对到共有序列库的序列是400条；则400条按照1:2划分给两个物种，即大肠杆菌133.33条，志贺氏菌266.67条。

实施例3与常规建库和筛选方法进行比较

按照常规方法构建数据库，即使用4个下载的大肠杆菌和志贺氏菌物种的基因组序列，各去除质粒序列后，将序列合并在一起，作为此4种细菌的微生物基因组参考数据库。

对此数据库，使用blast软件命令makeblastdb构建比对库。使用上述7个模拟数据，与此比对库进行blast比对，对于比对结果进行筛选。

第一种方法，使用每条reads的所有比对结果中，均为同一个物种时，即为此reads的来源物种；如果单条reads的所有比对结果不是同一个物种时，reads舍弃。结果如下：

第二种方法，每条reads仅保留bitscore值最高的比对结果，此结果来源的物种即为此reads的来源物种。结果如下：

以上三种方法鉴定大肠杆菌和志贺氏菌的结果请见图1和图2。鉴定结果表明，本发明方法与传统的建库筛选方法相比，在物种的鉴定上有较为明显的优势，可以更为真实的反映模拟数据中的物种情况；即使物种出现了部分的变异突变，本发明方法的鉴定结果也更为接近模拟数据物种的真实情况。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，但本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。