一种病原微生物参考知识库的最大多样性聚类构建方法

摘要

本发明涉及生物技术领域，具有涉及一种病原微生物参考知识库的最大多样性聚类构建方法。该构建采用贪心缩放算法，通过特定方法计算冗余基因组的多序列相似性，去除扩展相似度高的序列，可以保证在去除大量冗余信息的同时，最大程度保留微生物基因组的多样性。实验表明，本发明提供的病原微生物知识库的最大多样性聚类方法准确性高、检测效率高且数据冗余率低，尤其适用于大规模的病原微生物全基因组序列。

说明书

技术领域

本发明涉及生物技术领域，特别涉及一种病原微生物参考知识库的最大多样性聚类构建方法。

背景技术

高通量测序的持续发展导致可用的基因组测序数据空前增长，如此庞大的数据集给数据的存储和计算带来了巨大的挑战。为了从测序样本中精准鉴定病原微生物的种类，需要提前构建临床级病原微生物参考基因组知识库，知识库的质量直接决定了致病原检测结果的精确性，而知识库的大小决定了检测的速度。专家共识建议知识库中病原体的种类应该涵盖细菌、真菌、寄生虫、病毒、支原体或衣原体、分枝杆菌等，涵盖的种类建议超过2万种。建设时不仅需要包括染色体，还需要考虑线粒体、转录组序列及非编码序列等。因此，知识库的完整性和准确性是临床应用检测的关键。

当前病原微生物的知识库建立的方法是采用从国际公共数据库中通过获取、整合、筛选等一系列处理过程进行建立。通常知识库的建立方式有三种，分别是冗余方法、菌株代表序列法和打断重组法。第一，冗余方法是将获取的数据经过序列校正后全部纳入参考知识库，保留了物种的所有信息，可以实现所有物种的准确比对和筛查，但缺点是知识库的数据量庞大且重复序列多，极大降低了检测速度，这对实际临床应用患者往往是不可取的。第二，菌株代表序列法是从所有的序列中选取每个菌种的一条代表序列，其优点是可以实现快速检测。但是，菌株代表序列法也有明显的缺点，由于同种菌株的序列存在变异差异，舍弃大部分的基因组序列会丢失较多的病原微生物信息，造成检测的假阳性非常高。第三，打断重组法是将同一菌种的基因组进行打断后重新组装，可以尽可能保留基因组的亚种变异信息，但是由于病原微生物的种类繁多，从种的水平将基因组打断后重新组装非常困难。同时，在组装过程中打断重组法也会引入两个新的问题：一是由于受组装软件的限制，组装效率受打断后的读长（reads）长度、组装软件参数、基因组GC含量、基因组长度等多个因素影响，导致组装性能低；二是对于通过组装获取的病原微生物基因组的准确性需要进一步验证。

当前病原微生物知识库的构建方法存在的难题较多，主要原因是可用参考基因组冗余度高且数据量大，数据量甚至达到了TB级别，这极大影响了致病菌的检测速度，从而导致感染类疾病不能及时得到精准筛查。急危重症等临床级应用知识库通常需要人工筛选，耗费了巨大的人工成本。随着新的病原微生物基因组不断被发布，病原知识库需要不断地迭代更新，如果更新不及时会影响病原微生物的检测精度，导致漏检现象发生。因此，建立非冗余病原检测的知识库用于病原微生物检测参考知识库，达到非冗余基因组的准确性，提高致病微生物的检测速度和精度，节省计算和存储资源，具有重要的临床价值。

发明内容

有鉴于此，本发明提供一种病原微生物参考知识库的最大多样性聚类构建方法。该方法准确性高、检测效率高且数据冗余率低，是尤其适用于大规模的病原微生物全基因组序列。

本发明提供一种病原微生物参考知识库的最大多样性聚类构建方法，包括：

步骤1：收集病原微生物的全基因组序列，获得所有病原微生物的总冗余基因组R；

步骤2：根据病原微生物的界的分类对所述总冗余基因组R进行分箱，将每一种类的病原微生物的冗余基因组R

（1）将冗余基因组R

（2）将原始冗余集R

（3）计算S

（4）将获得的核心基因组集C

C1中的每条序列对应步骤（2）中的S

（5）将C

依次类推，将每一轮比对结束后获得的核心序列C

步骤3：按照步骤（1）~（5）分别构建得到每类病原微生物的参考知识库Q

一些实施方案中，所述分箱具体为：将同类型病原微生物的基因组序列进行合并，不同种的病原微生物的基因组进行分箱。

一些实施方案中，所述索引为：对所述原始冗余集中的每条序列进行编号，取值范围为：0~(n-1)，其中n为原始冗余集的序列总数，获取每条序列的编号和序列长度。

一些实施方案中，所述并行排序为：根据序列的编号对进程号进行求余计算，对原始冗余集中的每一个冗余知识库按照序列长度进行降序排列。每一个冗余知识库即同类型的原始病原微生物基因组序列集。

一些实施方案中，所述并行排序的判定条件为：seq_num%p_num_ = c_id，其中，seq_num表示序列编号，取值范围为：0~（n-1）之间的整数，其中n为原始冗余集的序列总数；p_num表示进程总数或采用的CPU核数，取值范围为：大于等于1的正整数；c_id为当前进程的进程号，取值范围为：0~（n-1）之间的整数，其中n为原始冗余集的序列总数；%表示取余运算符，即序列编号对进程总数取余后等于当前进程号，则表示该基因组存储在当前进程的内存区。

一些实施方案中，所述核心基因组序列S1为原始冗余集Rm’中长度最长的序列。

一些实施方案中，采用子序列扩展最大匹配的比对策略对子序列集A进行序列聚类。

进一步，所述子序列扩展最大匹配的比对策略选自以下情况之一：

①如果子序列m1、子序列m2被子序列m3所覆盖，则m3作为扩展后的匹配子序列；

②如果子序列m3被子序列m1的5'端和子序列m2的3'端所覆盖，则将m1、m2、m3进行拼接，得到子序列w，最终将子序列w最为扩展后的匹配子序列；

③如果子序列m1的5'端和子序列m2的3'端有重叠部分，且覆盖子序列m3 ，则将m1和m2进行拼接，得到子序列w，最终将子序列w最为扩展后的匹配子序列；

④如果子序列m1和子序列m2没有重叠部分，且不覆盖m3，则采用奖惩得分矩阵对m1、m2和m3分别向两端进行扩展。扩展过程允许空位匹配，如果出现错配或者空位则赋予负值得分，则将子序列m1、m2和m3均作为扩展后的匹配子序列。

一些具体实施例中，本发明病原微生物参考知识库的构建方法，包括：

步骤a、冗余基因组获取：收集并获取病原微生物的所有全基因组序列，即冗余基因组R；

步骤b、分箱操作：按照病原微生物的界的分类对冗余基因组R中的序列按照同类合并、异类分箱的原则进行操作；然后将每一类病原微生物的冗余基因组R

所述参考知识子库Q

b1、索引提取：对冗余基因组R

b2、序列并行排序：根据编号对进程号的求余计算，对冗余基因组R

b3、核心基因组设定：将原始冗余集Rx’中长度最长的基因组作为核心基因组S1；

b4、基因组扩展：采用子序列扩展最大匹配的比对策略对原始冗余集Rx’进行序列聚类：首先，将R

①如果子序列m1、子序列m2被子序列m3所覆盖，则m3作为扩展后的匹配子序列；

②如果子序列m3被子序列m1的5'端和子序列m2的3'端所覆盖，则将m1、m2、m3进行拼接，得到子序列w，最终将子序列w最为扩展后的匹配子序列；

③如果子序列m1的5'端和子序列m2的3'端有重叠部分，且覆盖子序列m3 ，则将m1和m2进行拼接，得到子序列w，最终将子序列w最为扩展后的匹配子序列；

④如果子序列m1和子序列m2没有重叠部分，且不覆盖m3，则采用奖惩得分矩阵对m1、m2和m3分别向两端进行扩展。扩展过程允许空位匹配，如果出现错配或者空位则赋予负值得分，则将子序列m1、m2和m3均作为扩展后的匹配子序列；

b5、扩展序列相似性计算：通过上述扩展策略，得到最大精确匹配集A的扩展最大精确匹配子序列集A'；计算核心基因组序列和剩余序列之间的扩展相似性S

b6、核心基因组标记：将满足上一步骤中扩展相似性阈值的序列被标记为冗余基因组序列，否则标记为核心基因组序列；获得核心基因组集C1和冗余序列集R1；

b7、贪心缩放聚类：将获得的核心基因组集C

C1中的每条序列对应步骤b4中的S

将C

依次类推，将每一轮比对结束后获得的核心序列C

步骤c、独立集整合：将每类病原微生物的参考知识库Q

本发明提供的病原微生物参考知识库的最大多样性聚类构建方法包括11个模块，分别如下：冗余基因组获取：收集病原微生物的全基因组序列；分箱操作：将获取的全基因组进行同类合并，异类分箱；索引提取：抽取基因组编号和序列长度；序列并行排序：采用多进程模式对获取的冗余序列排序；基因组扩展：采用最大匹配子序列扩展策略进行序列聚类；扩展序列相似性计算：计算核心基因组序列和剩余序列之间的扩展相似性；核心基因组设定：将每轮次聚类中最长的序列最为核心基因组S1；贪心缩放聚类：采用贪心缩放算法，将核心基因组作为参考基因组完成所有冗余序列聚类；冗余序列分离：合并核心基因组，得到全基因组参考序列集；独立集整合：合并的独立集，获得病原微生物全基因组参考知识库。该知识库去除了全基因组数据集中的高冗余序列，提高致病微生物的检测速度和精度，降低了知识库冗余率，节省了计算和存储资源。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。

图1为实施例2中流感病毒检测方法准确率对比示意图;

图2为实施例2流感病毒检测的准确率和冗余率示意图，图中采用“#1”表示最长代表序列方法，横坐标所有数字表示不同的序列相似性;

图3为实施例2中流感病毒检测时间示意图，图中采用“#1”表示最长代表序列方法，横坐标所有数字表示不同的序列相似性;

图4为实施例3中大肠杆菌建立方法对比示意图;

图5为实施例3中大肠杆菌检测时间对比示意图;

图6为实施例3中大肠杆菌加速比对比示意图;

图7为实施例3中大肠杆菌文件大小示意图;

图8为实施例3中大肠杆菌序列数示意图;

图9为实施例3中大肠杆菌比对上的reads数示意图。

具体实施方式

本发明公开了病原微生物参考知识库的最大多样性聚类构建方法，本领域技术人员可以借鉴本文内容，适当改进工艺参数实现。特别需要指出的是，所有类似的替换和改动对本领域技术人员来说是显而易见的，它们都被视为包括在本发明。为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的较佳实施例。但是,本发明可以以许多不同的形式来实现,并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。

本发明提供的一种病原微生物参考知识库的构建方法中所用原料及试剂均可由市场购得。

下面结合实施例，进一步阐述本发明：

实施例1 病原微生物参考知识库的构建

构建方法具体如下：

步骤a、冗余基因组获取：收集并获取病原微生物的所有全基因组序列，即冗余基因组R；

步骤b、分箱操作：按照病原微生物的界的分类对冗余基因组R中的序列按照同类合并、异类分箱的原则进行操作；然后将每一类病原微生物的冗余基因组R

所述参考知识子库Q

b1、索引提取：对冗余基因组R

b2、序列并行排序：根据编号对进程号的求余计算，对冗余基因组R

b3、核心基因组设定：将原始冗余集Rx’中长度最长的基因组作为核心基因组S1；

b4、基因组扩展：采用子序列扩展最大匹配的比对策略对原始冗余集Rx’进行序列聚类：首先，将R

①如果子序列m1、子序列m2被子序列m3所覆盖，则m3作为扩展后的匹配子序列；

②如果子序列m3被子序列m1的5'端和子序列m2的3'端所覆盖，则将m1、m2、m3进行拼接，得到子序列w，最终将子序列w最为扩展后的匹配子序列；

③如果子序列m1的5'端和子序列m2的3'端有重叠部分，且覆盖子序列m3 ，则将m1和m2进行拼接，得到子序列w，最终将子序列w最为扩展后的匹配子序列；

④如果子序列m1和子序列m2没有重叠部分，且不覆盖m3，则采用奖惩得分矩阵对m1、m2和m3分别向两端进行扩展。扩展过程允许空位匹配，如果出现错配或者空位则赋予负值得分，则将子序列m1、m2和m3均作为扩展后的匹配子序列；

b5、扩展序列相似性计算：通过上述扩展策略，得到最大精确匹配集A的扩展最大精确匹配子序列集A'；计算核心基因组序列和剩余序列之间的扩展相似性S

b6、核心基因组标记：将满足上一步骤中扩展相似性阈值的序列被标记为冗余基因组序列，否则标记为核心基因组序列；获得核心基因组集C1和冗余序列集R1；

b7、贪心缩放聚类：将获得的核心基因组集C

C1中的每条序列对应步骤b4中的S

将C

依次类推，将每一轮比对结束后获得的核心序列C

步骤c、独立集整合：将每类病原微生物的参考知识库Q

实施例2

一、数据获取

本发明首先从NCBI的Refseq核酸数据库中获取Viral数据集，共得到病毒全基因组序列12,182条。对得到所有的病毒基因组进行筛选，最终得到78条流感病毒全基因组序列（FLU数据集）用于后续的实验验证，其中，最大长度：2368nt，最小长度：838nt，平均长度：1721nt，这里nt表示核苷酸（Nucleotide）。

二、数据集划分

将FLU数据集划分为22个子集，其中根据序列相似性划分的为21个，步长为1%，分别为： 80%、81%、82%、83%、84%、85%、86%、87%、88%、89%、90%、91%、92%、93%、94%、95%、96%、97%、98%、99%、100%，以及只有参考基因组序列的数据子集1个。

三、测试集生成

（1）在测试集的准备中，本发明随机筛选了3条流感病毒基因组，如表1所示。因此在最长代表序列方法中，选取的代表序列为NC_007375.1，长度为2341nt。

（2）对于上述随机筛选的3条序列进行随机打断成75bp长度的reads，共200000条序列，即生成20M的流感病毒短reads集。然后，将该reads集作为流感病毒检测的测试输入数据。

表1 用于生成测试reads的流感病毒基因组列表

四、知识库建立分析与对比

将上述流感病毒短reads集与划分的22个FLU数据子集进行比对，采用的方法是BWA软件的mem短序列比对。为了评估聚类方法的性能，本发明采用多个指标进项衡量，具体包括：准确率、加速比、检测时间、冗余率、数据库大小、全基因组序列数及比对上的reads数。

（1）准确率：流感病毒短reads集中比对到数据集上的reads数的占比；

（2）检测时间：将将高通量测序reads数据比对到所建参考知识库上的CPU执行时间；

（3）数据库大小：菌种序列数据文件的大小；

（4）全基因组序列数：所建知识库中包含的基因组数量；

（5）加速比：以冗余方法的检测时间为基准，除以当前方法的检测时间的比值作为加速比的值，单位为1；

（6）冗余率：每个数据集的文件大小与冗余方法文件大小的比值；

（7）检测时间：流感病毒短reads集比对到数据集上的CPU执行时间；

（8）比对上的reads数：采用BWA软件的mem功能模块比对结果中mapped对应的reads数。

从图1的测试结果看出，在流感病毒测试集上，最长代表序列方法的检测准确率仅为29.80%，而聚类方法的准确率为97.50%，冗余方法为100%，说明聚类方法的准确率远优于最长代表序列方法，与冗余方法接近。

表2为流感病毒知识库建立方法在检测结果的对比，从表中可以看出，聚类方法在99%相似性条件下检测时间为6.29s，全基因组序列数减少了7条，比对上的reads数为195,065条，远高于最长代表序列方法（59,577条）。

表2 流感病毒知识库建立方法测试结果

图2结果显示随着序列相似性的增加检测准确率和冗余率均呈逐渐上升的趋势。序列相似性在80%-99%区间变化时，99%的相似性条件下准确率最高，且冗余率上升迟缓。需要注意的是，虽然最长代表序列方法的冗余率最低（1.70%），但准确率也是最低；在100%序列相似性条件下等同于最长代表序列方法。

对于流感病毒的检测时间评估，从图3可以看出，聚类方法在99%序列相似性的条件下（6.29s）与最长代表序列方法（4.14s）的CPU执行时间均相近，证明本发明的聚类方法的在保持比对时间稳定的前提下，准确率上具有较好的优势。需要注意的是，上述两种方法与冗余方法（4.65s）的CPU执行时间较长，原因是BWA的整体执行时间较短，因此服务器计算任务波动影响较大导致的。

实施例3

一、数据获取

本发明首先从NCBI的Refseq数据库中获取细菌（Bacteria）数据集中的大肠杆菌（Escherichia coli）进行提取，去冗余前的数据量约为6.7GB。Taxonomy ID是562，对得到所有的细菌基因组进行筛选，最终得到4192条流大肠杆菌全基因组序列，最大长度：5,942,969nt，最小长度：763nt，平均长度：1,689,358nt，该数据集将被用于后续的实验验证。

二、数据集划分

考虑到数据量较大，本实施例中划分为3个区间进行验证，分别为： 90%相似性的聚类方法、100%相似性的冗余方法，以及只有参考基因组序列的最长代表序列方法。

三、测试集生成

（1）在测试集的准备中，本发明随机筛选了3条大肠杆菌基因组，如表3所示。因此在最长代表序列方法中，选取的代表序列为NZ_CP027599.1，长度为5,942,969nt。

（2）对于上述随机筛选的3条序列进行随机打断成75bp长度的reads，共200000条序列，即生成20M的大肠杆菌短reads集。然后，将该reads集作为大肠杆菌检测的测试输入数据。

表3 用于生成测试reads的大肠杆菌基因组列表

四、知识库建立分析与对比

从图4的测试结果看出，在大肠杆菌测试集上，最长代表序列方法的检测准确率为95.51%，而聚类方法和冗余方法的准确率均为100%（90%序列相似性条件下），说明聚类方法的准确率远优于最长代表序列方法，与冗余方法相同。在知识库的冗余率上，最小序列方法最低，冗余率为0.08%，但是会丢失部分基因组信息，导致准确率也是最低。聚类方法的冗余率为19.64%，远低于冗余方法，而且准确性与冗余方法保持一致。

图5为大肠杆菌检测时间对比示意图，从图中可以看出聚类方法（90%相似性）的检测时间为75.884s，冗余方法的检测时间为260.864s，聚类方法的检测时间远低于冗余方法。

图6为大肠杆菌加速比对比示意图，从图中可以看出聚类方法的检测速度是冗余方法的3.44倍，最长代表序列方法是冗余方法的48.14倍。

图7为大肠杆菌文件大小示意图，从图中可以看出聚类方法（90%相似性）的知识库大小为1343.13MB，冗余方法的知识库大小为6838.47MB，冗余序列减少了超过80%。

图8为大肠杆菌序列数示意图，从图中可以看出聚类方法（90%相似性）的知识库序列数为1428，冗余方法的知识库序列数位4192，冗余序列减少了2764条。

图9为大肠杆菌比对上的reads数示意图，从图中可以看出根据聚类方法（90%相似性）获取全基因组作为参考知识库时，比对上的reads数为200,000，冗余方法获取全基因组作为参考知识库时，比对上的reads数也为200,000，最长代表序列方法获取全基因组作为参考知识库时，比对上的reads数为191093条，说明聚类方法可以达到与冗余方法相同的水平，远优于最长代表序列方法。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。