CN108334750B - 一种宏基因组数据分析方法及系统 - Google Patents

Abstract

本发明涉及一种宏基因组数据分析方法及系统。本发明所述数据分析方法和系统基于k‑mer算法获得样本的初步物种鉴定结果，并基于所述初步物种鉴定结果，抽取部分或全部支持序列，采用blast算法对初步物种鉴定结果进行验证，判断初步物种鉴定结果是否为报告检出物种。本发明所述方法及系统能够降低假阳性，在短时间内快速、准确地获得样本的报告检出物种，并与多种主流测序平台兼容，适用于第二代测序技术和第三测序技术；本发明所述方法、系统还能够对样本的耐药基因和耐药突变位点进行准确的鉴定并映射到报告检出物种。进一步地，本发明所述系统可用于鉴定病原体微生物，尤其是心内膜炎病原体，克服其难以培养的缺陷。

Description

一种宏基因组数据分析方法及系统

技术领域

本发明涉及生物信息学领域，具体而言，涉及一种宏基因组数据分析方法及系统。

背景技术

宏基因组(Metagenome)，又称群落基因组，是指特定小生境中全部微小生物遗传物质的总和。宏基因组学(metagenmomics)是指，直接应用基因组学技术对小生境中的微生物群落进行研究，且不需要对单一菌株进行分离培养的学科。

不同于以往的微生物学分析手段，宏基因组学分析不需要筛选得到各微生物群落的培养物，而是直接测定样品中所有微生物的核酸序列，以分析微生物群落的生长情况。宏基因组学分析可避免由于环境改变而引起的微生物序列变化所带来的偏差，特别适合用于鉴定在常规条件下难以培养的微生物，例如，极端微生物或某些难以培养的病原体。以导致心内膜炎的病原微生物为例，常规培养不仅耗时长(2～4周)，而且培养成功率不高，阴性率高达30％。对于许多感染性心内膜炎患者而言，不论是血液或者是心脏赘生物样本都不能通过培养给出物种和耐药性鉴定结果，影响治疗方案和复发干预措施的指定。而宏基因组学可以通过直接对患者的生物样本进行测序，绕过培养步骤，获得生物样本的病原体群落信息。

随着测序技术的发展，宏基因组学技术已有长足发展，但其在数据分析方面仍然存在瓶颈，制约其发展和应用，具体表现为：

1)宏基因组高通量检测具有高敏感性，但检测结果中的假阳性过多，特异性差，不能满足特异性要求高的鉴定方法的需求，例如，对病原微生物的临床鉴定。

2)现有的宏基因组测序数据分析方法尚难以在保证鉴定结果准确性的基础上，大幅度加快分析速度，缩短分析时间。

3)现有的宏基因组数据分析平台兼容性差，不能普遍适用于各类测序场景。

4)现有的宏基因组分析技术尚不能将物种鉴定与功能基因分析有机地整合，无法提供更为全面、深度加工的信息分析结果。

有鉴于此，特提出本发明。

发明内容

本发明的目的在于提供宏基因组数据分析方法和系统，以解决上述技术问题中的至少一个。

为了实现本发明的上述目的，特采用以下技术方案。

一种宏基因组数据分析方法，所述方法包括：

A1、将从样本获得的宏基因组测序序列与第一核酸数据库进行比对，基于k-mer算法获得初步物种鉴定结果；

A2、针对所述初步物种鉴定结果中的每个物种，统计支持序列的总数目，如果支持序列的总数目≥预定值，则随机抽取部分支持序列，构成验证序列集，如果支持序列的总数目＜预定值，则抽取全部支持序列，构成验证序列集；基于blast算法，在第二核酸数据库中对所述验证序列集中的序列进行物种鉴定，当验证序列集中50％以上序列的鉴定结果与所述初步物种鉴定结果一致，则所述初步物种鉴定结果通过验证，所述初步物种鉴定结果即为样本的报告检出物种。

本发明所述方法在鉴定样本的物种来源时，创造性地采用“两步法”对宏基因组数据进行分析。第一步，将测序序列与微生物基因组数据库进行比对，基于k-mer算法在短时间内处理大量测序数据，快速地获得初步的物种鉴定结果，该鉴定结果包含较高的假阳性结果，常规的过滤参数设置难以将其剔除。第二步，在初步鉴定结果的基础上，本发明抽取部分支持序列，采用更为准确但速度较慢的blast算法进行验证，能够大大消除k-mer算法缺陷造成的假阳性结果，同时，由于blast比对对象为初步鉴定结果中抽取的部分支持序列，计算量显著下降，能够在短时间内完成数据分析，既达到验证目的又控制计算时间。另外，在支持序列的总数目较少，低于预定值的情况下，为避免验证序列集不能满足数据统计要求，本发明所述方法还进一步限定，当支持序列的总数目小于预定值时抽取全部支持序列，用于构成验证序列集，从而确保有足够数量的支持序列用于blast算法的验证。根据实施例1～3的实验结果可知，本发明所述方法对于第二代测序数据，由于读长较短，会在同属异种的物种中出现一些错误或难以判断的情况，造成少量的假阳性结果，但对于读长更长的新一代测序数据，基本没有假阳性结果产生，能得到十分准确的鉴定结果。

在一些具体的实施方式中，所述第一核酸数据库与所述第二核酸数据库相同或不同，优选地，所述第一核酸数据库与所述第二核酸数据库不同；更优选地，所述第一核酸数据库为微生物基因组数据库，所述第二核酸数据库为NCBI的nt数据库。

在一些具体的实施方式中，为平衡数据的统计学意义和数据量带来的统计学压力，所述预定值设置为既满足验证序列集对统计数据的要求，又满足blast算法对运算速度的要求；优选地，所述预定值设置为150～250，更优选为200。

在一些具体的实施方式中，随机抽取的所述部分支持序列的数目≥预定值，优选地，随机抽取的所述部分支持序列的数目等于预定值。

在一些具体的实施方式中，所述步骤A2还包括在抽取支持序列之前对鉴定结果进行过滤，优选地，所述过滤包括将支持序列数据量占比低于0.2％、长读长支持序列数目低于2条、短读长支持序列数目低于10条的鉴定结果剔除。

在一些具体的实施方式中，所述步骤A2进行blast比对时，E value＜E^-30。

在一些具体的实施方式中，为获取样本中各微生物的相对定量信息，本发明所述方法还包括：A3、计算每个报告检出物种的支持序列的长度总和，经所述报告检出物种的基因组大小标准化后，得到所述报告检出物种富集程度的参考值，以所述参考值为基准，计算每个报告检出物种在所述样本中的相对丰度。

本发明还涉及一种宏基因组数据分析，所述方法包括：B1、基于blast算法，将样本的宏基因组测序序列与耐药基因数据库的核酸序列信息进行比对，获得耐药基因的初步鉴定结果；B2、对所述耐药基因的初步鉴定结果进行过滤，只保留支持序列对所述耐药基因的序列覆盖度＞30％且不集中分布于所述耐药基因的3’或5’端的的鉴定结果，即为报告检出耐药基因。

在一些具体的实施方式中，所述步骤B2在过滤前先将不符合E value＜E^-30的耐药基因鉴定结果剔除。

本发明上述方法不但能够鉴定样品中的微生物物种，还能鉴定样品的耐药性。并且，本发明所述方法对耐药基因数据准确性的兼容性强。对于耐药基因数据库来说，对于序列两端的注释不一定完全准确，故本发明所述方法在设计中对3’和5’端序列的比对过滤均存在一定的兼容度。这种兼容表现在两个方面。一是控制序列对基因的覆盖区域，过滤掉只有3’或5’有序列覆盖的鉴定结果，二是对3’和5’端比对的宽容，即如果中间区域有符合要求(准确度和覆盖度均满足)对3’和5’端的不准确给予宽容处理。

在一些具体的实施方式中，所述方法还包括：B1’、基于blast算法，将宏基因组测序序列与耐药基因数据库的耐药突变位点氨基酸序列进行比对，获得耐药突变位点的初步鉴定结果；B2’、对所述耐药突变位点的初步鉴定结果进行过滤，只保留长读长支持序列数目＞2条或短读长支持序列数目＞3条的鉴定结果，即为报告检出耐药突变位点。

在一些具体的实施方式中，所述方法还包括：B3、整合所述报告检出耐药基因和/或所述报告检出耐药突变位点，获得样本的全局耐药结果。

在一些具体的实施方式中，所述方法还包括执行权利要求1所述步骤A1～A2，或A1～A3，获得样本的报告检出物种，将B3所述全局耐药结果映射到报告检出物种。

本发明上述方法将宏基因组的物种鉴定和耐药基因/耐药突变位点的鉴定进行有机整合。其中，结合物种鉴定结果，本发明所述方法能够将耐药突变基因或耐药突变位点信息映射到报告检出物种，获得具体的耐药菌株，对微生物的培养、临床用药等具有重要的指导意义。

在一些具体的实施方式中，所述方法还包括在进行物种鉴定、耐药基因鉴定和/或耐药突变位点鉴定前，对所述宏基因组测序序列进行质量评估和/或去宿主基因组。

在一些具体的实施方式中，所述质量评估包括：根据测序平均质量分数、质量分数分布以及测序长度中的至少一种评估测序数据的可用性，对可用数据中的低质量部分进行过滤处理；优选地，所述过滤处理包括：对于长读长测序数据，保留长度＞500bp且测序质量分数＞6的测序数据；对于短读长测序数据保留长度＞50bp且质量分数＞20的测序数据。

在一些具体的实施方式中，所述去宿主基因组包括：将所述宏基因组测序序列与宿主基因组序列进行比对，剔除能比对到所述宿主基因组的序列，保留不能比对到所述宿主基因组的序列，用于后续的数据分析。

在一些具体的实施方式中，所述去宿主基因组采用bwa算法进行比对，去短读长测序序列和长读长测序序列均兼容且适宜。

在一些具体的实施方式中，所述质量评估在所述去宿主基因组之前执行。

在一些具体的实施方式中，所述方法为非诊断目的和/或非治疗目的。

在一些具体的实施方式中，所述方法中的宏基因组测序序列为长读长序列或短读长序列。

本发明前述方法对测序数据的可用性进行评估，并将可用数据中的低质量部分过滤，能够提高数据分析的准确性，进一步地，本发明前述方法在进行质量评估时能够针对不同的测序数据采用不同的过滤处理方式，对测序数据具有高度的兼容性。本发明前述方法还包括去宿主基因组步骤，该步骤能够极大程度地降低后续分析的数据量，能够加快数据分析的速度，尤其能够减轻blast算法的运算压力。

本发明还涉及用于前述宏基因组数据分析方法的系统，所述系统包括物种鉴定模块，所述物种鉴定模块包括初步鉴定模块和验证模块，其中：

所述初步鉴定模块包括接收单元和比对单元；

所述接收单元用于接收从样本中获得的宏基因组测序序列；

所述比对单元用于将所述测序序列与第一核酸数据库进行比对，基于k-mer算法获得初步物种鉴定结果；

所述验证模块包括验证序列集单元、物种鉴定单元和验证单元；

其中，所述验证序列集单元：针对每个初步的物种鉴定结果，统计支持序列的总数目，如果支持序列的总数目≥预定值，则随机抽取部分支持序列，构成验证序列集，如果支持序列的总数目＜预定值，则抽取全部支持序列，构成验证序列集；

所述物种鉴定单元用于在第二核酸数据库中基于blast算法对验证序列集进行物种鉴定；

所述验证单元，用于比较初步物种鉴定结果与验证序列集的物种鉴定结果的一致性，如果验证序列集中50％以上序列的鉴定结果与初步物种鉴定结果一致，则判定通过验证，所述初步物种鉴定结果即为报告检出物种；

优选地，所述系统还包括相对丰度计算模块，所述相对丰度计算模块包括参考值计算单元和相对丰度计算单元；

其中，所述参考值计算单元，用于计算每个报告检出物种的支持序列的长度总和，并将所述支持序列的长度总和用所述报告检出物种的基因组大小标准化，获得所述报告检出物种富集程度的参考值；

所述相对丰度计算单元，用于以所述参考值为基准，计算每个报告检出物种在所述样本中的相对丰度。

本发明还涉及用于前述宏基因组数据分析方法的系统，所述系统包括耐药基因鉴定模块，所述耐药基因鉴定模块包括初步鉴定单元和过滤单元；

其中，所述初步鉴定单元，用于基于blast算法将宏基因组测序序列与耐药基因数据库的核酸序列信息进行比对，获得耐药基因的初步鉴定结果；

所述过滤单元，用于过滤所述耐药基因的初步鉴定结果，只保留支持序列对所述耐药基因的序列覆盖度＞30％且不集中分布于所述耐药基因的3’或5’端的宏基因组测序序列的鉴定结果，即为报告检出耐药基因；优选地，在过滤前先将E value＜E^-30的鉴定结果剔除；

优选地，所述系统还包括耐药突变位点鉴定模块，所述耐药突变位点鉴定模块包括初步鉴定单元和过滤单元；

其中，所述初步鉴定单元，用于基于blast算法将所述宏基因组测序序列与耐药基因数据库的耐药突变位点氨基酸序列进行比对，获得耐药突变位点的初步鉴定结果；

所述过滤单元，用于过滤所述耐药突变位点的初步鉴定结果，只保留长读长支持序列数目＞2条或短读长支持序列＞3条的鉴定结果，即为报告检出耐药突变位点；

更优选地，所述系统还包括全局耐药报告模块，所述全局耐药报告模块包括整合单元，所述整合单元用于整合所述报告检出耐药基因和/或所述报告检出耐药突变位点，获得样本的全局耐药结果；

最优选地，所述系统还包括前述物种鉴定模块和/或相对丰度计算模块，所述全局耐药报告模块还包括映射单元，所述映射单元用于将全局耐药结果映射到报告检出物种，所述报告检出物种通过执行所述物种鉴定模块获得。

在一些具体的实施方式中，所述系统还包括质量评估模块和/或去宿主基因组模块；

其中，所述质量评估模块包括评估单元和过滤单元；

所述评估单元，用于根据测序平均质量分数、质量分数分布以及测序长度中的至少一种评估测序数据的可用性；

所述过滤单元，用于过滤可用数据中的低质量部分；优选地，所述过滤处理用于：对于长读长测序数据，保留长度＞500bp且测序质量分数＞6的测序数据；对于短读长测序数据，保留长度＞50bp、质量分数＞20的测序数据；

所述去宿主基因组模块包括比对单元和过滤单元；

其中，所述比对单元用于将所述宏基因组测序序列与宿主基因组序列进行比对；

所述过滤单元用于根据比对结果，剔除能比对到所述宿主基因组的序列，保留不能比对到所述宿主基因组的序列，用于后续的数据分析。

优选地，所述系统先执行质量评估模块，再执行去宿主基因组模块。

优选地，所述系统先执行质量评估模块和去宿主基因组模块，再执行物种鉴定模块、耐药基因鉴定模块和/或耐药突变位点鉴定模块。

优选地，在一些具体的实施方式中，所述分析系统用于鉴定环境微生物或病原微生物，优选地，所述环境微生物为极端微生物，所述病原微生物为心内膜炎病原体；优选地，所述系统接收的宏基因组测序序列获取自第二代测序平台或第三代测序平台，更优选地，所述宏基因组测序序列获取自第三代测序平台。

本发明还涉及前述系统在鉴定环境微生物或病原微生物中的应用，所述环境微生物为极端微生物，所述病原微生物为心内膜炎病原体；优选地，所述系统接收的宏基因组测序序列获取自第二代测序平台或第三代测序平台，更优选地，所述宏基因组测序序列获取自第三代测序平台。

术语定义

本说明书中使用的术语具有以下定义。

术语“第一核酸数据库”是指数据库形式和数据库大小与k-mer算法兼容的核酸数据库。

术语“第二核酸数据库”是指数据库形式和数据库大小与blast算法兼容的核酸数据库。

术语“支持序列”是指物种鉴定或者耐药鉴定过程中，经算法识别或比对后，判定为与参考基因组或耐药数据库参考序列一致的序列，这个序列的集合为某一物种或耐药基因鉴定的总支持序列。

术语“长读长”是指能够产生＞500bp序列的测序平台所产生的数据，例如，由Nanopore或Pacbio等测序平台产生的数据。

术语“短读长”是指由能够产生50-300bp序列的测序平台所产生的数据，例如，由Illumina，454或BGI等测序平台产生的数据。

术语“第二代测序平台或第三代测序平台”是指包含但不限于SMRT，Nanopore，Illumina，ION torrent或BGI等测序平台。

技术效果

与现有技术相比，本发明的有益效果为：

1)本发明所述方法和系统基于k-mer算法和blast算法对宏基因组数据进行分析处理，能够降低假阳性，在短时间内快速、准确地获得样本的报告检出物种；进一步，本发明所述方法和系统还能够计算报告检出物种在样本中的相对丰度，获得相对定量信息；更进一步，本发明所述方法和系统还对测序数据的质量进行评估以及去除宿主基因组，能够进一步提高数据分析的准确性、降低数据的处理量，缩短处理时间。

2)本发明所述方法和系统还能够鉴定耐药基因和耐药突变位点，获得样本的全局耐药性，并结合物种鉴定结果，映射到报告检出物种，获得具体的耐药微生物菌株；进一步地，本发明所述方法和系统在鉴定耐药基因时对耐药基因数据准确性的兼容性强，减少耐药基因数据库序列两端的注释不准确的缺陷。

3)本发明所述方法和系统通过算法的选取、流程以及参数的设计可以适用于各类主流测序平台(如SMRT，Nanopore，Illumina，ION torrent或BGI)，应用场景更为普遍，并充分发挥各平台的特点和优势。

4)本发明所述方法和系统解决了难培养微生物，例如心内膜炎病原体的检测问题，检测全面，快速，可覆盖细菌、真菌、少见病原菌、单一病原体感染和混合病原体感染等各种类型。同时，本发明所述方法和系统还包括心内膜炎病原体的全局耐药基因位点检测和耐药株的鉴定结果，可以为临床的治疗提供帮助，预防手术后病人再次感染。

具体实施方式

下面将结合实施例对本发明的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购买获得的常规产品。

实施例1基于Nanopore测序平台的心脏赘生物样本的宏基因组检测与数据分析

收集7例临床培养阴性的感染性心内膜炎患者瓣膜置换手术取出的心脏赘生物样本A1～A7，-80℃冰箱进行保存。

所得样本按照以下流程提取核酸：从冰箱取出赘生物样本先在室温放置30分钟，然后使用灭菌剪刀将赘生物样本剪碎，使用天根试剂盒(TIANamp Micro DNA kit)按照说明书进行核酸提取。

提取后的核酸样本按照以下流程进行建库和上机测序，建库方案选取牛津纳米孔公司提供的1D Native barcoding protocol：

1)使用g-TUBE(Covaris)将1.2μg核酸样本在5000转/分钟的条件下打断1分钟，获得片段化的DNA；

2)核酸的末端修复：在45μl片段化的DNA中加入3μl Ultra II End-prep enzymemix(NEB)、7μl Ultra II End-prep reaction buffer(NEB)和5μl nuclease-free water，0.2ml PCR管中进行混匀，20℃反应5分钟过后65℃再反应5分钟；

3)加barcode：每个样本取末端修复后的核酸500ng，各加入2.5μl NativeBarcode和25μl Blunt/TA Ligase Master Mix，混匀后21℃反应30分钟；

4)将上一步加好barcode的所有样本共取700ng，加入20μl Barcode Adapter Mix(BAM)和10μl Quick T4 DNA Ligase，混匀后室温反应10分钟；

5)构建好的文库上机(Nanopore MinION测序仪)进行测序。

7个样本中，样本A1与A2平行建库两次，有2个重复，则7个样本产生9套测序数据。测序数据的分析流程如下：

1)将Nanopore MinION测序产生数据通过开源的basecalling算法Albacore，将电信号转换为碱基信号，得到fastq格式的序列信息。随后根据添加的barcode序列拆分成属于不同样本的序列集合，同时去除接头序列。

2)测序质量鉴定。计算统计数据质量分数(quality score)，对于Nanopore的长读长数据，过滤掉长度＜500bp、平均测序质量＜6的序列。

3)去宿主序列。通过bwa开源算法，将过滤后的序列比对到人类基因组(基因组版本HG38)，保留比对失败的序列，进入下一步分析。

4)采用“两步法”对去宿主后序列进行病原体物种鉴定。在第一步中，序列被比对到构建好的包含5591种细菌以及236种真菌的微生物基因组数据库中，基因参考序列来源于NCBI的refGENE数据库。Centrifuge软件在该数据库集合中加以鉴定，并过滤掉支持序列数目≤2，富集比例＜0.2％的鉴定病原体。过滤后的结果进入第二步的验证流程。在第二步的验证中，使用的验证数据库为NCBI的NT数据库。用开源算法blastn将需要验证鉴定结果的支持序列在该数据库做比对。对于每个待鉴定物种，随机筛选200条支持序列作为验证集(不足200条的，将全部支持序列作为验证集)，如验证集中超过50％的序列都可以比对到目标物种上(满足E value<E^-30)，则通过验证。经过两次比对后的鉴定结果为最终报告检出物种。计算每个物种支持序列的长度总和，再用病原体基因组大小标准化序列长度和，得到病原体富集程度的参考数值，以该参考值为基准，计算每个鉴定物种在所述样本中的相对丰度比例。

5)耐药基因及耐药位点突变的检测。将去宿主后得到的非宿主序列数据用开源算法blast比对到耐药基因位点数据库CARD中的耐药基因核酸序列信息以及突变位点氨基酸信息，以E value<E^-30为比对结果的过滤标准。过滤后的比对结果中，针对每条耐药基因及耐药位点所在基因的鉴定结果，统计测序序列对目标基因的覆盖度，保留覆盖度＞30％且覆盖区域不仅仅集中在5’及3’端的结果。对于耐药突变位点的鉴定，除所在基因的比对结果符合过滤标准外，对于每个氨基酸突变位点需要有对应一致的核酸突变位点检出，且与注释信息突变一致的支持序列>2条。整合过滤后的耐药基因及耐药突变位点检出结果，根据数据库注释信息报告全局的耐药结果。再结合上一步物种鉴定信息中对序列所属物种的鉴定信息，如全局耐药基因中>2条序列对应到同一物种，则将该耐药检测结果对应到耐药株上。

数据分析结果：实施例1所述测序数据的基本信息如表1所示，病原体鉴定结果如表2所示，耐药性检测结果如表3所示。

根据表1所示结果可知，对于Nanopore平台产生的数据，测长较长，本数据集中平均长度在1k左右。

根据表2所示结果可知，7个样本均有病原体检出，平行重复数据检出物种一致。检出病原体均是临床证明与感染性心内膜炎高度相关的物种，特别是其中Coxiellaburnetii与Bartonella quintana是两种感染性心内膜炎高度相关的少见病原菌。每套数据的数据量均在100-1000M的数据水平，数据量不大，但仍可经过该方法得到准确的数据检出。

实施例1中A4，A5及A7样本的覆盖度符合耐药检出的需求，实际分析发现A5和A7样本有耐药基因检出，检出结果如表3所示。根据耐药基因检测结果，对应到药物层面，A5的Streptococcus oralis可能对链霉素、大环内酯物抗生素，以及林酰胺具有耐药性，而A7的Streptococcus anginosus可能对四环素、链霉素、大环内酯物抗生素，林酰胺具有耐药性。

表1

表2

表3

实施例2基于Ilumina测序平台的心脏赘生物样本的宏基因组检测与数据分析

以实施例1中的A1～A2为样本，提取基因组核酸并构建文库，采用Illumina HiseqPE150进行测序，所得测序数据通过去除接头和高N比例的序列后，得到fastq格式的序列信息。对各个样本分别进行如下的数据分析：

1)将Ilumina测序产生的fastq格式的数据去除接头和高N比例的序列后进入下一步的质量评估分析。

2)测序质量鉴定。该数据的建库读长为150，过滤掉长度＜100bp、平均测序质量＜25的序列，该数据前10个碱基GC比例较为异常则切除每条序列的前10个碱基。

3)去宿主序列。通过比对到人类基因组(基因组版本HG38)，将比对失败的序列保留，进入下一步分析。比对软件使用bwa开源算法；

4)采用“两步法”对去宿主后序列进行病原体物种鉴定。在第一步中，序列被比对到构建好的包含5591种细菌以及236种真菌的微生物基因组数据库中，基因参考序列来源于NCBI的refGENE数据库。Centrifuge软件在该数据库集合中加以鉴定，并过滤掉支持序列数目≤10，富集比例＜0.2％的鉴定病原体。过滤后的结果进入第二步的验证流程。在第二步的验证中，使用的验证数据库为NCBI的NT数据库。用开源算法blastn将需要验证鉴定结果的支持序列在该数据库做比对。对于每个待鉴定物种，随机筛选200条支持序列作为验证集(不足200条的，使用全部支持序列作为验证集)，如验证集中超过50％的序列都可以比对到目标物种上(满足E value<E^-30)，则通过验证。经过两次比对后的鉴定结果为最终报告检出物种。计算每个物种支持序列的长度总和，再用病原体基因组大小标准化序列长度和，得到病原体富集程度的参考数值，以该参考值为基准，计算每个鉴定物种在样本中的相对丰度比例。

5)耐药基因及耐药突变位点的检测。将去宿主后得到的非宿主序列数据用开源算法blast比对到耐药基因位点数据库CARD中的耐药基因核酸序列信息以及突变位点氨基酸信息，以E value<E^-30为比对结果的过滤标准。过滤后的比对结果中，针对每条耐药基因及耐药突变位点所在基因的鉴定结果，统计测序序列对目标基因覆盖度，保留覆盖度＞30％且覆盖区域不仅仅集中在5’及3’端的结果，且支持的序列数需≥3。对于耐药突变位点的鉴定，除所在基因的比对结果符合过滤标准外，对于每个氨基酸突变位点需要有对应一致的核酸突变位点检出，且与注释信息突变一致的支持序列>5条。整合过滤后的耐药基因及耐药位点检出结果，根据数据库注释信息报告全局的耐药结果。再结合上一步物种鉴定信息中对序列所属物种的鉴定信息，如全局耐药基因中>2条序列对应到同一物种，则将该耐药检测结果对应到耐药株上。

数据分析结果：样本A1和A2的测序数据量分别为55G和54G，将此数据集用于病原体物种鉴定，其结果如表4所示，耐药性检测结果如表5所示。

根据表4所示结果可知，A1和A2两个样本虽然病原体滴度非常低，但由于数据量的优势，对检出病原菌的基因组覆盖度均超过50％，深度为3倍左右。同时，由于较高的基因组覆盖度，两样本均有较好的耐药检出(参见表5、6)：由于读长较短，对于同源基因的耐药鉴定分辨率较低(参见表5)，但是由于同源基因指示的耐药表型一致，故对耐药预测结果不会造成影响(参见表6)。

表4

表5

表6

耐药药物类型	A1	A2
			tetracycline	tetM/S/O	tetM/S/O
macrolide	tetM/S/O	tetM/S/O
			lincosamide	ErmB	ErmB,RlmA(II)
streptogramin	ErmB	ErmB
			fluoroquinolone	-	patB,pmrA

实施例3基于BGI测序平台的后心脏赘生物样本的耐药基因检测

以实施例1中的A1～A2为样本，提取基因组核酸并构建文库，采用BGI测序平台进行测序，将BGI测序产生数据针对各个样本分别进行如下的数据分析：

1)将BGI测序产生的fastq格式的数据去除接头和高N比例的序列后进入下一步的质量评估分析。

2)测序质量鉴定。该数据的建库读长为150，过滤掉长度＜100bp、平均测序质量＜25的序列。

3)去宿主序列。通过比对到人类基因组(基因组版本HG38)，将比对失败的序列保留，进入下一步分析。

4)采用“两步法”对测序序列进行病原体物种鉴定。在第一步中，序列被比对到构建好的包含5591种细菌以及236种真菌的微生物基因组数据库中，基因参考序列来源于NCBI的refGENE数据库。Centrifuge软件在该数据库集合中加以鉴定，并过滤掉支持序列数目≤10，富集比例＜0.2％的鉴定病原体。过滤后的结果进入第二步的验证流程。在第二步的验证中，使用的验证数据库为NCBI的NT数据库。用开源算法blastn将需要验证鉴定结果的支持序列在该数据库做比对。对于每个待鉴定物种，随机筛选200条支持序列作为验证集(不足200条的，使用全部支持序列作为验证集)，如验证集中超过50％的序列都可以比对到目标物种上(满足E value<E^-30)，则通过验证。经过两次比对后的鉴定结果为最终报告检出物种。计算每个物种支持序列的长度总和，再用病原体基因组大小标准化序列长度和，得到病原体富集程度的参考数值，以该参考值为基准，计算每个鉴定物种在所有鉴定结果中计算相对丰度比例。两步法的策略用第一步保证了快速的初步检测，而第二步保证了最终报告结果的准确性；

5)耐药基因检测。将去宿主后得到的非宿主序列数据用开源算法blast比对到耐药基因位点数据库CARD中的耐药基因核酸序列信息以及突变位点氨基酸信息，以E value<E^-30为比对结果的过滤标准。过滤后的比对结果中，对每条耐药基因及耐药位点所在基因的鉴定结果统计测序序列对目标基因覆盖度，保留覆盖度＞30％且覆盖区域不仅仅集中在5’及3’端的结果，且支持的序列数需≥3。对于耐药位点位点鉴定，除所在基因的比对结果符合过滤标准外，对于每个氨基酸突变位点需要有对应一致的核酸突变位点检出，且与注释信息突变一致且的支持序列>5条。报告整合过滤后的耐药基因及耐药位点检出结果，根据数据库注释信息报告全局的耐药结果。再结合上一步物种鉴定信息中对序列所属物种的鉴定信息，如全局耐药基因中>2条序列对应到同一物种，则将该耐药检测结果对应到耐药株上。

数据分析结果：样本A1～A7的测序数据的物种鉴定结果如表7所示，耐药鉴定结果如表8所示；不同平台下，样本A1～A7的耐药鉴定结果如表9所示。

根据表7所示检测结果可知，样本A1～A7在实施例3中所得的物种鉴定结果与实施例1完全一致。

根据表8～9所示结果可知，采用本发明的方法对不同数据集进行耐药基因和位点发现，虽然由于数据量和目的物种覆盖度的影响，耐药结果的检出有不稳定的现象，但不同数据集中相同样本的检出仍然十分一致。

表7

表8

表9

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

Claims (26)

1.一种宏基因组数据分析方法，其特征在于，所述方法包括：

A1、将样本的宏基因组测序序列与第一核酸数据库进行比对，基于k-mer算法获得初步物种鉴定结果；

A2、针对所述初步物种鉴定结果中的每个物种，统计支持序列的总数目，如果支持序列的总数目≥预定值，则随机抽取部分支持序列，构成验证序列集，如果支持序列的总数目＜预定值，则抽取全部支持序列，构成验证序列集；基于blast算法，在第二核酸数据库中对所述验证序列集中的序列进行物种鉴定，当验证序列集中50%以上序列的鉴定结果与所述初步物种鉴定结果一致，则所述初步物种鉴定结果通过验证，所述初步物种鉴定结果即为样本的报告检出物种；

A3、计算每个报告检出物种的支持序列的长度总和，经所述报告检出物种的基因组大小标准化后，得到所述报告检出物种富集程度的参考值，以所述参考值为基准，计算每个报告检出物种在所述样本中的相对丰度。

2.根据权利要求1所述的宏基因组数据分析方法，其特征在于，所述方法还包括：

B1、基于blast算法，将样本的宏基因组测序序列与耐药基因数据库的核酸序列信息进行比对，获得耐药基因的初步鉴定结果；

B2、对所述耐药基因的初步鉴定结果进行过滤，只保留支持序列对所述耐药基因的序列覆盖度＞30%且不集中分布于所述耐药基因的3’或5’端的鉴定结果，即为报告检出耐药基因。

3.根据权利要求2所述的宏基因组数据分析方法，其特征在于，所述方法还包括：

B1’、基于blast算法，将样本的宏基因组测序序列与耐药基因数据库的耐药突变位点氨基酸序列进行比对，获得耐药突变位点的初步鉴定结果；

B2’、对所述耐药突变位点的初步鉴定结果进行过滤，只保留长读长支持序列的数目＞2条或短读长支持序列的数目＞3的鉴定结果，即为报告检出耐药突变位点。

4.根据权利要求2所述的宏基因组数据分析方法，其特征在于，所述方法还包括：B3、整合所述报告检出耐药基因，获得样本的全局耐药结果。

5.根据权利要求3所述的宏基因组数据分析方法，其特征在于，所述方法还包括：B3、整合所述报告检出耐药基因和所述报告检出耐药突变位点，获得样本的全局耐药结果。

6.根据权利要求4或5所述的宏基因组数据分析方法，其特征在于，所述方法还包括B4、根据所述样本的报告检出物种，将步骤B3所述全局耐药结果映射到报告检出物种。

7.根据权利要求1所述的宏基因组数据分析方法，其特征在于，所述方法还包括在进行物种鉴定前，对所述宏基因组测序序列进行质量评估和去宿主基因组。

8.根据权利要求2所述的宏基因组数据分析方法，其特征在于，所述方法还包括在进行物种鉴定和耐药基因鉴定前，对所述宏基因组测序序列进行质量评估和去宿主基因组。

9.根据权利要求3所述的宏基因组数据分析方法，其特征在于，所述方法还包括在进行物种鉴定、耐药基因鉴定和耐药突变位点鉴定前，对所述宏基因组测序序列进行质量评估和去宿主基因组。

10.根据权利要求7~9任一项所述的宏基因组数据分析方法，其特征在于，所述质量评估包括：根据测序平均质量分数、质量分数分布以及测序长度中的至少一种评估测序数据的可用性，对可用数据中的低质量部分进行过滤处理。

11.根据权利要求10所述的宏基因组数据分析方法，其特征在于，所述过滤处理包括：对于长读长测序数据，保留长度＞500 bp且测序质量分数＞6的测序数据；对于短读长测序数据保留长度＞50bp且质量分数＞20的测序数据。

12.根据权利要求7~9任一项所述的宏基因组数据分析方法，其特征在于，所述去宿主基因组包括：将所述宏基因组测序序列与宿主基因组序列进行比对，剔除能比对到所述宿主基因组的序列，保留不能比对到所述宿主基因组的序列，用于后续的数据分析。

13.根据权利要求7~9任一项所述的宏基因组数据分析方法，其特征在于，所述质量评估在所述去宿主基因组之前执行。

14.一种宏基因组数据分析系统，其特征在于，所述系统包括物种鉴定模块，所述物种鉴定模块包括初步鉴定模块和验证模块，其中：

所述初步鉴定模块包括接收单元和比对单元；

所述初步鉴定模块的所述接收单元用于接收从样本中获得的宏基因组测序序列；

所述初步鉴定模块的所述比对单元用于将所述测序序列与第一核酸数据库进行比对，基于k-mer算法获得初步物种鉴定结果；

所述验证模块包括验证序列集单元、物种鉴定单元和验证单元；

其中，所述验证序列集单元用于：针对每个初步的物种鉴定结果，统计支持序列的总数目，如果支持序列的总数目≥预定值，则随机抽取部分支持序列，构成验证序列集，如果支持序列的总数目＜预定值，则抽取全部支持序列，构成验证序列集；

所述物种鉴定单元用于在第二核酸数据库中基于blast算法对验证序列集进行物种鉴定；

所述验证单元用于比较初步物种鉴定结果与验证序列集的物种鉴定结果的一致性，如果验证序列集中50%以上序列的鉴定结果与初步物种鉴定结果一致，则判定通过验证，所述初步物种鉴定结果即为报告检出物种；

所述系统还包括相对丰度计算模块，所述相对丰度计算模块包括参考值计算单元和相对丰度计算单元；

其中，所述参考值计算单元，用于计算每个报告检出物种的支持序列的长度总和，并将所述支持序列的长度总和用所述报告检出物种的基因组大小标准化，获得所述报告检出物种富集程度的参考值；

所述相对丰度计算单元，用于以所述参考值为基准，计算每个报告检出物种在所述样本中的相对丰度。

15.根据权利要求14所述的宏基因组数据分析系统，其特征在于，所述系统还包括耐药基因鉴定模块，所述耐药基因鉴定模块包括初步鉴定单元和过滤单元；

其中，所述耐药基因鉴定模块的所述初步鉴定单元，用于基于blast算法将样本的宏基因组测序序列与耐药基因数据库的核酸序列信息进行比对，获得耐药基因的初步鉴定结果；

所述耐药基因鉴定模块的所述过滤单元，用于过滤所述耐药基因的初步鉴定结果，只保留支持序列对所述耐药基因的序列覆盖度＞30%且不集中分布于所述耐药基因的3’或5’端的宏基因组测序序列的鉴定结果，即为报告检出耐药基因。

16.根据权利要求15所述的宏基因组数据分析系统，其特征在于，所述系统还包括耐药突变位点鉴定模块，所述耐药突变位点鉴定模块包括初步鉴定单元和过滤单元；

其中，所述耐药突变位点鉴定模块的所述初步鉴定单元，用于基于blast算法将所述样本的宏基因组测序序列与耐药基因数据库的耐药突变位点氨基酸序列进行比对，获得耐药突变位点的初步鉴定结果；

所述耐药突变位点鉴定模块的所述过滤单元，用于过滤所述耐药突变位点的初步鉴定结果，只保留长读长支持序列数目＞2条或短读长支持序列＞3条的鉴定结果，即为报告检出耐药突变位点。

17.根据权利要求16所述的宏基因组数据分析系统，其特征在于，所述系统还包括全局耐药报告模块，所述全局耐药报告模块包括整合单元，所述整合单元用于整合所述报告检出耐药基因和所述报告检出耐药突变位点，获得样本的全局耐药结果。

18.根据权利要求17所述的宏基因组数据分析系统，其特征在于，所述全局耐药报告模块还包括映射单元，所述映射单元用于将全局耐药结果映射到报告检出物种，所述报告检出物种通过执行所述物种鉴定模块获得。

19.根据权利要求14~18任一项所述的宏基因组数据分析系统，其特征在于，所述系统还包括质量评估模块和去宿主基因组模块；

其中，所述质量评估模块包括评估单元和过滤单元；

所述质量评估模块的所述评估单元，用于根据测序平均质量分数、质量分数分布以及测序长度中的至少一种评估测序数据的可用性；

所述质量评估模块的所述过滤单元，用于过滤可用数据中的低质量部分；

所述去宿主基因组模块包括比对单元和过滤单元；

其中，所述去宿主基因组模块的所述比对单元用于将所述宏基因组测序序列与宿主基因组序列进行比对；

所述去宿主基因组模块的所述过滤单元用于根据比对结果，剔除能比对到所述宿主基因组的序列，保留不能比对到所述宿主基因组的序列，用于后续的数据分析。

20.根据权利要求19所述的宏基因组数据分析系统，其特征在于，所述质量评估模块的所述过滤单元用于：对于长读长测序数据，保留长度＞500 bp且测序质量分数＞6的测序数据；对于短读长测序数据保留长度＞50bp且质量分数＞20的测序数据。

21.根据权利要求19所述的宏基因组数据分析系统，其特征在于，所述系统先执行质量评估模块，再执行去宿主基因组模块。

22.根据权利要求14~18任一项所述的宏基因组数据分析系统，其特征在于，所述系统用于鉴定环境微生物。

23.根据权利要求22所述的宏基因组数据分析系统，其特征在于，所述环境微生物为极端微生物。

24.根据权利要求14~18任一项所述的宏基因组数据分析系统，其特征在于，所述系统用于鉴定病原微生物。

25.根据权利要求24所述的宏基因组数据分析系统，其特征在于，所述病原微生物为心内膜炎病原微生物。

26.根据权利要求14~18任一项所述的宏基因组数据分析系统，其特征在于，所述系统接收的宏基因组测序序列获取自第二代测序平台或第三代测序平台。