病原微生物基因组数据库及其建立方法与流程

本发明涉及宏基因组技术领域，特别是涉及一种病原微生物基因组数据库及其建立方法。

背景技术：

病原宏基因(转录)组测序是一种应用于临床病原感染检测的新兴技术，具有检测病原种类广、灵敏度高、准确度高、时效快等优点，正在逐渐成为临床一线的病原感染检测手段。而病原宏基因组技术核心之一是病原微生物基因组数据库，数据库的质量直接影响着病原宏基因组检测的病原微生物种类数、准确度及分析性能等。

病原微生物基因组数据库是由多种物种的基因组构成的，大部分都收集于ncbi等公共数据库，经过筛选、加工构建而成。在公共数据库中，每个物种可能测有多个不同的菌株基因组，如何选择、加工这些菌株基因组来代表该物种是一个难题。

目前大多数宏基因组分析的流程都只是随机或通过聚类等分析后选取其中一株作为该物种的代表基因组，或不加选择地将全部菌株纳入到数据库中，这两种方法都各有优缺点。只选取一株作为代表的方法的优点是构建的数据库数据量比较小，分析时效快，但是缺点也比较突出，因为微生物基因组进化较快，同一个物种的菌株间的基因组存在差异，不同地区、不同时间，甚至是同时分离的不同菌株测序的基因组就可能存在差异，部分进化较快的物种菌株间的基因组差异可达3％以上。

目前病原宏基因组在临床实际应用中，大多数病原体的检出序列数在几十到几百之间，甚至是个位数，基因组覆盖率在0.1％以下，而菌株间基因组的突变率已高达3％以上，甚至含有特有序列，因此对于只选取一株作为物种代表基因组的方法，0.1％覆盖率的检出难以覆盖突变区域或特有序列区，常常会造成漏检而得出假阴的检测结果。

而将同一物种所有菌株基因组纳入数据库的方法则可以有效避免这种情况的漏检，但这方法的缺点也尤为明显。一方面，将所有菌株基因组纳入以后，数据库数据量变得很大，导致分析时间变得很长，甚至需要一天以上，这对临床应用的时效要求是难以接受的，往往早一个小时就可能多救治一个病患，同时对计算服务器或集群的资源需求也极大增加，分析成本大大增加；另一方面，来源于公共数据库的菌株基因组的测序质量参差不齐，部分菌株含有污染序列，甚至含有分类错误的菌株，如果不加以筛选过滤，就容易导致假阳的结果产生，给临床诊治带来极大的困扰。

技术实现要素：

基于此，有必要针对上述问题，提供一种病原微生物基因组数据库及其建立方法，采用该方法得到的基因组数据库，既有准确率高的优点，又具有分析时效短，节约成本的优势。

一种病原微生物基因组数据库的建立方法，包括以下步骤：

数据获取：获取病原微生物基因组数据；

菌株基因组筛选：按照预定筛选规则选取物种菌株基因组；

去除质粒序列：去除上述得到的菌株基因组中存在的质粒序列；

过滤：按照预定过滤规则，去除标注信息有误、染色体组装不完整，以及分类错误的菌株，得到该物种的参考菌株基因组；

构建融合基因组：将上述参考菌株基因组打断，得到序列集，将上述序列集中的序列通过比较去除冗余，再对非冗余数据集进行组装得到重新组装的基因组，并过滤低质量的contigs，根据长度将contigs从大到小重新拼接，得到该物种的融合基因组；

组库：重复上述步骤，得到预定物种的融合基因组，汇总，即得病原微生物基因组数据库。

上述数据库的建立方法，通过下载微生物基因组，筛选高质量基因组，去除质粒序列，去除标注、分类错误基因组，过滤基因组低质量或污染片段，物种菌株间基因组去冗余后融合，最后将所有物种融合基因组进行组库得到病原微生物基因组数据库。

可以理解的，上述去除质粒序列步骤中，可根据具体情况调整，对于古菌、细菌等有质粒序列的基因组，进行去除质粒序列处理，没有质粒的基因组则不需要去除。

在其中一个实施例中，所述数据获取步骤中，所述病原微生物基因组数据来源于patric数据库和/或ncbi的refseq和genbank数据库中分类为archaea、bacteria、fungi、protozoa和viral的基因组数据。上述数据库中的病原微生物基因组数据具有较高的数据质量，上述数据既包括基因组序列，也包括描述该序列的说明性文件。

在其中一个实施例中，所述菌株基因组筛选步骤中，所述筛选规则为：

如物种具有若干个菌株基因组数据，则选取组装完成且组装质量高的菌株基因组数据；

如物种仅有单个菌株基因组数据，则直接选取该菌株基因组数据。

具体的，在上述数据库中，组装完成指组装完成度为“completegenome”，质量高指质量为“good”或有文献支持的菌株基因组。

在其中一个实施例中，所述过滤步骤中，所述过滤规则包括勘误过滤规则，所述勘误过滤规则为：

如物种为多菌株基因组的物种，统计各菌株基因组的chromosome和contig数量，若contig数量大于chromosome数量，剩余contig无信息标注的，则舍弃该菌株基因组；

如物种为单菌株基因组的物种，计算该菌株基因组的contig数及每个contig的长度，进而获得该单菌株基因的n90，舍弃n90以外的contigs，或者，舍弃长度小于设定值的contigs。优选地，病毒中，该设定值为300；细菌、古菌、真菌、寄生虫中，该设定值为1000；

上述n90以外的contigs具体指：通过组装软件直接拼接成的序列片段称作contig，按照contig长度从大到小排序，并从大到小进行累加，当累加值大于或等于该物种基因组长度的90％时所加上的那个contig的长度定义为n90，小于该长度的contigs即为n90以外的contigs。

上述勘误过滤规则用于去除标注信息有误、染色体组装不完整的数据。

在其中一个实施例中，所述构建融合基因组步骤中，具体包括以下步骤：

基因组打断：将所述参考菌株基因组打断为长度为n且移步步长为k的序列集；例如，可截取基因组序列第1个碱基至第n个碱基之间的序列为第一条序列，然后起始位置定位至第k个碱基，截取第k至k+n位碱基间的序列为第二条序列，依次类推，对于截取的第n条序列，其位置为(n-1)k至(n-1)k+n；

序列集去冗余：以上述序列集中的一条序列为参考序列集，将其余序列分别与该参考序列集比较，如100％匹配，则将该条序列舍弃，若非100％匹配则加入到参考集中，重复上述过程，历遍整个序列集，得到非冗余序列集；

非冗余序列集组装：对上述非冗余序列集进行组装，组装完成后统计contigs的长度，指定长度l，过滤长度小于l的contigs，得到重新组装的基因组；

基因组contigs拼接：根据长度将contigs从大到小排序，contigs间使用m个连续n连接，融合为完整的基因组，得到的该物种的融合基因组。

可以理解的，上述n指在数据库中以符号“n”表示序列间隙和不同菌株的连接，也可根据需求，选用除“a、t、c、g”外等其它不会产生混淆的符号代替。

在其中一个实施例中，所述n为30～10000的自然数；所述k为1～n的自然数。

在其中一个实施例中，所述n选自：100～500。

在其中一个实施例中，所述非冗余序列集组装步骤中，所述长度l为100～10000。优选100-1000，更优选300-500。

在其中一个实施例中，所述聚类序列集坐标回溯步骤中，m为大于4的自然数。优选的，m选自10～100。

本发明还公开了上述的建立方法得到的病原微生物基因组数据库。

上述病原微生物基因组数据库，既整合了物种所有可靠菌株的差异序列，保留了丰富的物种菌株信息，能有效避免病原宏基因(转录)组检测结果的假阳性、假阴性，提高检测准确度；同时也去除了物种内的冗余序列，大大减少了数据库的数据量，减少了分析计算资源的需求，缩短了分析时间，降低了分析成本。

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

本发明的一种病原微生物基因组数据库的建立方法，通过对病原微生物基因组数据的筛选、过滤和融合，去除了冗余序列，得到的病原微生物基因组数据库，既整合了物种所有可靠菌株的差异序列，保留了丰富的物种菌株信息，能有效避免病原宏基因(转录)组检测结果的假阳性、假阴性，提高检测准确度；同时也去除了物种内的冗余序列，大大减少了数据库的数据量，减少了分析计算资源的需求，缩短了分析时间，降低了分析成本。

附图说明

图1为实施例2中5株模拟测序集数据比对到三种方法构建的数据库的比对准确率；

图2为实施例2中5株模拟测序集数据比对到三种方法构建的数据库的准确率箱型图；

图3为实施例2中5株模拟测序集数据比对到三种方法构建的数据库分析时间；

图4为实施例2中5株模拟测序集数据比对到三种方法构建的数据库分析时间箱型图。

具体实施方式

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

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

实施例1

一种病原微生物基因组数据库，通过以下方法建立：

一、数据获取

从patric下载细菌基因组数据，具体如下：

patric全称为美国病原微生物资源整合中心，网站上收录了目前已知的绝大部分细菌类病原微生物基因组数据，可从其ftp服务器上下载所有收录的细菌基因组数据。

在patric网站ftp服务器(ftp://ftp.patricbrc.org/)上面下载了所有分类为archaea、bacteria的基因组数据，并下载了对应的基因组信息统计文件patric_genome.txt，共包含227577株基因组信息。

二、菌株基因组筛选

根据文件表头信息，选取“public”列为“true”，“genomestatus”列为“complete”，“genomequality”列为“good”的基因组，筛选后得到13537株基因组。

以下通过肺炎克雷伯菌(klebsiellapneumoniae)为例，进行物种内的菌株筛选与处理说明。

三、去除质粒序列

去除上述得到的菌株基因组中存在的质粒序列，具体为：

经上述初步筛选后，肺炎克雷伯菌共有289株菌株基因组，根据基因组fasta文件里面的序列名称标注，带有“plasmid”或“plasmid”关键词即可确认该序列为质粒序列，将该质粒序列去除，得到289株不含质粒序列的基因组。

四、过滤

按照预定过滤规则，去除标注信息有误、染色体组装不完整，以及分类错误的菌株，具体为：

1、去除标注信息有误、染色体组装不完整的基因组。

统计上述289株不含质粒序列的菌株基因组的contig数，肺炎克雷伯菌的染色体(chromosome)为1个，所以将所有contig数大于1的菌株都舍弃，得到110株只含有一个contig的基因组。

五、构建融合基因组

以上述得到的110株肺炎克雷伯菌优质参考菌株基因组为例，构建肺炎克雷伯菌的融合基因组。

1、基因组打断

将110株肺炎克雷伯菌优质菌株基因组按照长度为300bp，步长为1bp打断，即n＝300，k＝1，得到序列集。

2、序列集去冗余

将上述序列集进行去冗余处理，第一条序列为参考集，从第二序列开始历遍整个序列集，去相似度为100％的序列，得到非冗余序列集。

3、非冗余序列集组装

将非冗余序列集转换为fastq格式文件，使用spades软件对fastq文件进行组装，组装完成后统计所有contigs长度，过滤长度小于100的contigs，即l＝100，得到重新组装的基因组。

4、基因组contigs拼接

根据长度将上述重新组装的基因组contigs按照从大到小进行排序，contigs间用10个“n”碱基进行连接，连接得到的即为肺炎克雷伯菌的融合基因组，可用于微生物基因组组建。

六、组库

重复上述步骤，得到上述下载13537株基因组所覆盖的物种的融合基因组，汇总，即得病原微生物基因组数据库。

实施例2

为了评价上述实施例1所构建的肺炎克雷伯菌的融合基因组效果，对肺炎克雷伯菌的未经处理的所有菌株基因组、肺炎克雷伯菌的ncbi参考菌株基因组、上述肺炎克雷伯菌的融合基因组进行准确度和分析时效等进行分析比较。

一、数据量评估比较

首先在ncbigenome数据库查询下载肺炎克雷伯的referencegenome，用于构建“一个物种只选取一株菌株基因组”方法(传统方法一)的数据库(肺炎克雷伯菌参考菌株基因组)，在ncbiassembly数据库查询肺炎克雷伯菌的所有组装基因组，仅refseq数据库上就有7280条组装结果，总碱基数据量超过20g，一个物种的数据量就如此之大，普通计算服务器已经难以处理，为了顺利分析，只选取了组装状态为“completegenome”的菌株，共347株，从中随机选取5株出来作为临床模拟菌株，生成长度为75bp，深度为2x的模拟测序数据集，用于评价测试。

剩余342株菌株基因组用于构建“一个物种选取所有菌株基因”(传统方法二)方法的数据库(肺炎克雷伯菌所有菌株基因组)。

上述传统方法一构建的数据库数据量大小为5.5m，传统方法二所构建的数据库数据量大小为1.8g，本专利实施例1方法所构建的数据库(肺炎克雷伯菌融合基因组)数据量大小为54m。

相对于传统方法一，本专利方法构建得到的数据库数据量上仅增加了8.8倍，却包含了肺炎克雷伯菌几乎所有的真实基因组序列，而传统方法二的数据是本专利方法的33倍，且两者包含的真实基因组序列接近。在数据库数据量大小方面的比较已初步看出本专利方法的优势。

二、数据分析效果比较

使用bwa软件mem模块将模拟测序数据集分别比对到上述构建的三种数据库中，统计分析时间及比对的准确度。

1、准确率。

结果如图1和图2所示，图1为5株模拟测序集数据比对到三种方法构建的数据库的准确率示意图，图2为5株模拟测序集数据比对到三种方法构建的数据库的准确率的波动情况示意图。

从图中可以看出，5株模拟测序数据集中，传统方法一的准确率平均值为81.81％，传统方法二的准确率平均值为98.88％，本专利方法的准确率平均值为98.73％。即本专利方法的准确率与传统方法二几乎一致，且显著高于传统方法一。

2、分析时间。

在分析时间方面，如图3、4所示，图3为5株模拟测序集数据比对到三种方法构建的数据库分析时间示意图，图4为5株模拟测序集数据比对到三种方法构建的数据库分析时间波动情况示意图。

从图中可以看出，方法一平均分析时间为1.59s，方法二平均分析时间为22.98s，本专利方法平均分析时间为2.53s。本专利方法时间上是方法一的1.6倍，是方法二的0.1倍。综上可以看出本专利方法具有高准确率、分析资源需求较低、分析时间较短等优点，显著优于两种传统的数据库构建方法。

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

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