CN111933218A - 一种优化的宏基因组binning分析微生物群落的方法 - Google Patents

Abstract

本发明公开了一种优化的宏基因组binning分析微生物群落的方法，包括对测序数据进行过滤，得到高质量测序数据，然后根据样本的来源及测序数据量的大小，选择不同的组装策略得到contigs，接着进行基因数据分析。与现有技术相比，本发明是一种针对微生物群落同时开展群落整体的生物信息分析和不依赖分离培养的“单菌”基因组分析。在宏基因组层面，提供了更贴合样本特征、测序数据量的高效优质组装算法，并包含丰富全面的信息分析内容，个性新颖的可视化。实现了宏基因组分析由群落到单菌的质变，方案包含可提高准确性的数据校正，包含全面完善的bin信息汇总，有利于更方便高效地筛选到有价值的目标bin，还包含系统完善的目标bin后续分析的挖掘思路。

Description

一种优化的宏基因组binning分析微生物群落的方法

技术领域

本发明涉及高通量测序技术和测序数据生物信息技术分析领域，特别涉及一种优化的宏基因组binning分析微生物群落的方法。

背景技术

群落微生物多种多样，如空气微生物、肠道微生物、土壤微生物等，随着高通量测序技术的普及，用于测序群落微生物多样性及功能的宏基因组学也随着兴起，宏基因组测序可获得一个群落内所有微生物的物种DNA信息。

宏基因组分箱(binning)是将宏基因组测序得到的混合了不同生物的序列或序列组装得到的contigs按物种分开归类的过程。。传统的单物种全基因组序列都是经纯培养之后，再进行全基因组de novo测序才获得的，但是环境中存在着大量的不可培养微生物，宏基因组分箱技术有助于获得不可培养微生物的全基因组序列，获得新物种的基因组序列和功能，预测未知物种的培养方法等等。

宏基因组binning技术的应用前景广泛，开始被更多的人关注和重视，也逐渐出现了很多宏基因组binning工具。现有的分箱工具主要如下表：

不同的分箱工具各有优势，默认设置下均需要消耗大量的计算资源，处理效率较为低下。

对于16S、ITS等主要的微生物群落研究方法，样本测序数据量一般为15M左右，而宏基因组测序数据量通常10G甚至100G。受数据量的影响，宏基因组分析普遍存在资源消耗较大、分析流程复杂、分析难度大、分析内容简略等问题。

虽然目前已有专利公开宏基因组测序及数据分析方法，但主要是分析整个环境微生物群落基因组信息的建库测序、数据质控等比较基础的方面。对于binning分析研究不够深入，主要存在如下不足：

宏基因组方面：

1)组装策略单一。未充分考虑样本类型(如自然环境样本、共生菌样本)、样本测序数据量(如6G、10G、20G，100G等)、生物学重复数量等，直接采用软件默认的组装参数，未进行测试评估，导致组装没有针对性、运算效率低、占用较多运行资源等弊端；

2)注释采用软件默认参数，未进行对比优化且注释策略单一。现有流程是基于基因进行物种注释，在相同测序深度、质控条件阈值下，此种方法非常依赖组装软件类型、组装参数、基因预测软件、预测参数的选择，存在较大波动，由于基因数一般上百万，故耗时耗资源；

3)分析结果片面。现有流程未充分对基因、功能层级、各类物种水平统一进行统计检验，未对各分类水平的物种均开展beta多样性分析。

binning方面，主要存在如下不足：

1)binning流程不严谨。虽然部分专利已有基于宏基因组的binning分析流程，但无耗时耗资源的数据校正部分，而这部分将直接影响bin的聚类、bin鉴定的准确性；

3)难以提供详细的bin基因注释。由于bin数量极多、每个bin基因数量繁多，软件不会输出每个bin每个基因的功能注释，并且软件默认结果没有基因所有7个层级的物种注释。但这部分信息对于用户使用分析有极高的参考价值，不仅可以提供基于bin的功能锁定目标bin的分析策略，还为用户建立了样本特有的所有高质量bin物种与基因的对应关系，可作为宏基因组、宏转录组等分析的个性化参考基因组(现有数据库中不存在的基因组，具有极强的样本针对性)，为客户批量开展bin的多组学分析提供了个性化的策略，可大大提高研究的针对性和研究深度。

3)bin物种注释结果参考价值极低。现有流程默认输出结果为每个物种的属水平物种或NA(属水平未知)，一方面，只展示属水平，过于抽象，用户没有直观的认识，无法获得有效的参考信息，需要手动依次去数据库查看其他分类层级的注释，而一般一个样本获得的高质量bin都可能上百，故工作量也较大；另一方面，仅停留在属水平，不够精细，bin在理论上是一个潜在的单菌基因组，即可作为具体的物种，但现有流程却未提供种水平的注释，故可能妨碍客户高效鉴定同属不同种的新物种；此外，属水平未知，不代表上级分类层级注释未知，很多特别新的物种，确实存在只能注释到目、科等分类层级，故这部分注释为NA的bin由于信息缺失导致无法进行科学的后续评估判断、甚至可能错过了发现目标新物种的可能；并且以主流客户的生信分析水平，无法实现对所有NA的bin重新分析鉴定；

4)可视化形式缺乏。仅提供了最简单的表格输出结果，未考虑用户在获得bin以后高效的信息提取，汇总需求。

5)现有流程分析缺乏对bin的后续数据挖掘体现。目前binning分析后即将bin作为了一个潜在的基因组而中断分析，或仅提供了bin重组装后的基因组序列，均未提供基于目标bin的数据挖掘思路和结果。而binning分析仅仅是中介，bin仅仅是后续分析的起点，真正的数据挖掘还未开始。由于bin分析是从群落总DNA跨至“基因组”，是跳入完全不同的分析领域，故对用户后期深入分析十分不友好，对用户生物学基础要求过高。

对现有的宏基因组binning技术进行进一步的改进，减少分析占用的资源，同时提高其分析结果的可用性，具有非常重要的意义。

发明内容

本发明的目的在于克服现有技术的至少一个不足，提供一种优化的宏基因组binning分析微生物群落的方法，特别是一种优化的宏基因组binning分析微生物群落中单菌株的方法。。

发明人在研究过程中意外发现，样本类型、样本测序数据量等对宏基因组组装的影响，经过样本类型对比、默认参数调整等整理出更有针对性、选择性更多、效率更高的宏基因组分析流程。

本发明所采取的技术方案是：

本发明的第一个方面，提供：

一种微生物群宏基因组的分析方法，包括：

测序数据过滤：对宏基因组测序所得原始数据进行过滤，并统计过滤前后碱基组成、reads数量、碱基质量，得到高质量测序数据；

数据组装：根据样本来源和测序数据量的大小，选择不同的数据组装策略，得到组装结果contigs，优选的，组装完成后过滤掉长度在500bp以下的短序列，其中：

策略1：若样本为非自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略2：若样本为非自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略3：若样本为自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

策略4：若样本为自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

基因数据分析：将S2步骤获得的数据结果进行基因数据分析，分析内容选自基因预测、基因丰度统计、核心基因和泛基因分析；

在一些实例中，所述基因预测的操作包括：利用MetaGeneMark软件对contigs进行基因预测，获得基因序列，然后采用CD-HIT软件，将基因序列相似度大于95％，序列比对区域大于90％的基因聚类为一个cluster，95％identity、90％coverage，选取最长的基因作为每个cluster的代表序列，称为Unigene，所有代表序列即为获得的非冗余基因集合Unigene catalogue。

在一些实例中，所述基因丰度统计的操作包括：利用bowtie2将高质量数据的reads重新比对Unigene，计算各样本中，每个基因的reads比对数目；过滤掉在各个样品中reads支持数目都≤2的基因，获得用于后续步骤分析的基因集合；基于比对上的reads数目及基因长度出发，按公式

计算得到各基因在各样品中的丰度信息，式中，r为比对上某基因k的reads数目，L为基因k的长度；

在一些实例中，所述核心基因和泛基因分析的操作包括：基于各样品的基因丰度表，从所有样本中，随机抽取1个样本，统计基因数目，多次随机抽取后，获得1个样本时的基因数目分布盒型图；再从所有样本中随机抽取2个样本，统计两样本基因交集和并集的基因数量，多次随机抽取后，获得2个样本时的基因数目分布盒型图；重复抽取足够样本，最终获得不同数目样品组合时的基因数目分布，得到核心基因和泛基因稀释曲线。

在一些实例中，重复抽取不少于6个，8个，或10个样本，最终获得核心基因和泛基因稀释曲线。样本的数量越多，得到的结果就越准确。

在一些实例中，所述非自然环境样本为来自动物体的微生物群，包括但不限于动物体的口腔、胃、肠、阴道、皮肤等的微生物群样本。

在一些实例中，所述自然环境样本为来自土壤、水体、空气的样本。

在一些实例中，所述的分析方法还包括功能注释：基于基因数据分析结果，将Unigenes比对到各类数据库进行基因功能注释和功能丰度统计。

在一些实例中，功能注释具体为将Unigenes序列通过DIAMOND软件，阈值evalue<＝1e-5，比对到数据库，同时结合基因丰度表格计算不同数据库比对结果的功能丰度信息。

本发明的第二个方面，提供：

微生物群宏基因组数据的binning分析方法，包括：

按本发明第一个方面的操作得到组装结果contigs；

去嵌合体校正：将样本测序的reads比对contig，计算一条contig上两个窗口的reads覆盖度，基于覆盖度差异程度判断是否为嵌合体，并进行切分校正；

binning：选择长度大于1.5kb的contig，将样本测序的高质量reads再次比对contig，计算每条contig的覆盖度和GC含量，使用metaBAT2进行binning，即进行contig聚类，软件鉴定判断后，将具有相似丰度和核酸组成的contig归为一个bin，即一个潜在的基因组。

在一些实例中，所述的binning分析方法进一步包括高质量bin筛选，具体包括：基于单拷贝基因集，评估每个bin的完整度、污染度。

在一些实例中，使用CheckM软件计算bin的完整度、污染度；完整度高于80％且污染度低于5％归类为高质量bin。

在一些实例中，基于contig的reads比对结果，统计每个bin的contig数量、测序深度、长度、GC特征并进行可视化分析。

在一些实例中，所述的binning分析方法还包括：

高质量bin物种注释，包括使用CheckM软件预测获得每个bin的基因序列信息，并使用DIAMOND软件将基因序列比对Nr库，提取基因在种水平的物种注释，并基于物种分类层级关系，统计获得基因对应的各分类水平的物种注释；或

高质量bin基因注释，包括基于基因序列，使用DIAMOND软件将基因序列比对Nr库，提取基因的功能描述信息。

在一些实例中，所述的binning分析方法还包括：

使用气泡图展示每个bin的基因特征，即以横轴展示基因长度、以纵轴展示基因序列比对的相似度、以点大小展示基因序列比对的得分、以点颜色展示基因在种水平的物种注释，点大小表示数据库比对的score得分，点越大，表示注释越可靠；和/或

选择目标bin构建bin基因组圈图，构建操作包括：

策略B1：使用bin的contig，即宏基因组拼接获得的contig，在完整度允许条件下，比对查找近缘物种，基于近缘物种特征，绘制基因组圈图；或

策略B2：将样本测序的reads再次比对回该bin的contig，提取比对上该bin的reads，使用soapdenovo软件进行重组装，基于新的组装结果绘制基因组圈图。

在一些实例中，样本测序的reads为其高质量reads。

本发明的第三个方面，提供：

一种微生物群宏基因组的分析系统，包括：

数据存储装置，用于存储待分析的高通量测序数据；

数据处理装置，用于对高通量测序数据进行组装和分析，其中，数据组装包括：

根据样本来源和测序数据量的大小，选择不同的数据组装策略，得到组装结果contigs，优选的，组装完成后过滤掉长度在500bp以下的短序列，其中：

策略1：若样本为非自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略2：若样本为非自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略3：若样本为自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

策略4：若样本为自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

数据分析包括：

基因数据分析：将S2步骤获得的数据结果进行基因数据分析，分析内容选自基因预测、基因丰度统计、核心基因和泛基因分析；

结果输出装置，用于输出处理后的数据。

在一些实例中，所述基因预测的操作包括：利用MetaGeneMark软件对contigs进行基因预测，获得基因序列，然后采用CD-HIT软件，将基因序列相似度大于95％，序列比对区域大于90％的基因聚类为一个cluster，95％identity、90％coverage，选取最长的基因作为每个cluster的代表序列，称为Unigene，所有代表序列即为获得的非冗余基因集合Unigene catalogue；

在一些实例中，所述基因丰度统计的操作包括：利用bowtie2将高质量数据的reads重新比对Unigene，计算各样本中，每个基因的reads比对数目；过滤掉在各个样品中reads支持数目都≤2的基因，获得用于后续步骤分析的基因集合；基于比对上的reads数目及基因长度出发，按公式

计算得到各基因在各样品中的丰度信息，式中，r为比对上某基因k的reads数目，L为基因k的长度。

在一些实例中，所述核心基因和泛基因分析的操作包括：基于各样品的基因丰度表，从所有样本中，随机抽取1个样本，统计基因数目，多次随机抽取后，获得1个样本时的基因数目分布盒型图；再从所有样本中随机抽取2个样本，统计两样本基因交集和并集的基因数量，多次随机抽取后，获得2个样本时的基因数目分布盒型图；重复抽取足够，如10个样本，最终获得不同数目样品组合时的基因数目分布，得到核心基因和泛基因稀释曲线。

在一些实例中，所述数据分析还包括功能注释：基于基因数据分析结果，将Unigenes比对到各类数据库进行基因功能注释和功能丰度统计。

在一些实例中，功能注释具体为将Unigenes序列通过DIAMOND软件，阈值evalue<＝1e-5，比对到数据库，同时结合基因丰度表格计算不同数据库比对结果的功能丰度信息。

本发明的第四个方面，提供：

一种微生物群宏基因组数据的binning分析系统，包括：

数据存储装置，用于存储待分析的高通量测序数据；

数据处理装置，用于对高通量测序数据进行组装和分析，其中，数据组装包括：

根据样本来源和测序数据量的大小，选择不同的数据组装策略，得到组装结果contigs，优选的，组装完成后过滤掉长度在500bp以下的短序列，其中：

策略1：若样本为非自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略2：若样本为非自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略3：若样本为自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

策略4：若样本为自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

数据分析包括：

去嵌合体校正：将样本测序的reads比对contig，计算一条contig上两个窗口的reads覆盖度，基于覆盖度差异程度判断是否为嵌合体，并进行切分校正；

binning：选择长度大于1.5kb的contig，将样本测序的高质量reads再次比对contig，计算每条contig的覆盖度和GC含量，使用metaBAT2进行binning，即进行contig聚类，软件鉴定判断后，将具有相似丰度和核酸组成的contig归为一个bin，即一个潜在的基因组。

结果输出装置，用于输出处理后的数据。

在一些实例中，所述进一步包括高质量bin筛选，具体包括：基于单拷贝基因集，评估每个bin的完整度、污染度；优选的，使用CheckM软件计算bin的完整度、污染度；完整度高于80％且污染度低于5％归类为高质量bin。

在一些实例中，还包括高质量bin物种注释，包括使用CheckM软件预测获得每个bin的基因序列信息，并使用DIAMOND软件将基因序列比对Nr库，提取基因在种水平的物种注释，并基于物种分类层级关系，统计获得基因对应的各分类水平的物种注释；或

高质量bin基因注释，包括基于基因序列，使用DIAMOND软件将基因序列比对Nr库，提取基因的功能描述信息。

在一些实例中，还包括所述数据分析还包括基于contig的reads比对结果，统计每个bin的contig数量、测序深度、长度、GC特征并进行可视化分析。

在一些实例中，使用气泡图展示每个bin的基因特征，即以横轴展示基因长度、以纵轴展示基因序列比对的相似度、以点大小展示基因序列比对的得分、以点颜色展示基因在种水平的物种注释，点大小表示数据库比对的score得分，点越大，表示注释越可靠；和/或

选择目标bin构建bin基因组圈图，构建操作包括：

策略B1：使用bin的contig，即宏基因组拼接获得的contig，在完整度允许条件下，比对查找近缘物种，基于近缘物种特征，绘制基因组圈图；或

策略B2：将样本测序的reads再次比对回该bin的contig，提取比对上该bin的reads，使用soapdenovo软件进行重组装，基于新的组装结果绘制基因组圈图。

在一些实例中，样本测序的reads为其高质量reads。

本发明的有益效果是：

本发明一些实例的方法，基于样本的具体情况选择不同的数据组装策略，数据的组装更有针对性，大幅提高了运算效率，对运行资源的占用率显著下降。

本发明一些实例的方法，得到的bin结果更为准确可靠。

本发明一些实例的方法，可以不同层级的注释信息，不仅可以提供基于bin的功能锁定目标bin的分析策略，还为用户建立了样本特有的所有高质量bin物种与基因的对应关系，可作为宏基因组、宏转录组等分析的个性化参考基因组(现有数据库中不存在的基因组，具有极强的样本针对性)，为客户批量开展bin的多组学分析提供了个性化的策略，可大大提高研究的针对性和研究深度。

本发明一些实例的方法，可以将结果可视化展示，便于直观显示结果。

与现有技术相比，本发明是一种针对微生物群落，同时开展群落整体的生物信息分析和不依赖分离培养的“单菌”基因组分析，一方面，在宏基因组层面，提供了更贴合样本特征、测序数据量的高效优质组装算法，并包含丰富全面的信息分析内容，个性新颖的可视化；另一方面，实现了传统宏基因组分析由群落到单菌的质变，方案包含可提高准确性的数据校正，包含全面完善的bin信息汇总，有利于更方便高效的筛选到有价值的目标bin，还包含系统完善的目标bin后续分析的挖掘思路。

附图说明

图1是本发明实施例中组装的contig长度分布图；

图2是本发明实施例中binning示意图；

图3是本发明实施例中基因注释气泡图。

具体实施方式

本发明一些实例使用的软件简介如下：

高质量测序数据的获取：

对原始高通量测序数据进行过滤，去除含adapter的reads；去除N的比例大于10％的reads；

数据经过过滤后通过分析碱基的组成及质量分布，进行碱基质量分析，去除低质量reads(质量值Q≤20的碱基数占整个read的50％以上)。一种具体方法是将过滤前后数据信息统计报表中，并将绘制各样品过滤前后碱基组成分布图，直观展示数据质量情况。

最后保留的数据称为高质量数据，后续分析均基于高质量数据进行。

Kmer长度的设置：

如果K-mer太长，可以拼接的更长，但是因为需要匹配的序列也变长了，所以拼接出的数量变少了。如果K-mer太短，因为需要匹配的序列短，所以会导致可以往下一步拼接的可能非常多，无法判断，而导致拼接中断，随意拼接的很短。具体如何设置Kmer的长度，是一项具有挑战性的任务。

自然环境样本包括但不限于土壤、水体、空气、污泥等非生物体来源的样本。非自然环境样本，特别指来源于生物体的样本，包括但不限于动物和人的胃、肠道，口腔，阴道，皮肤等样本。

宏基因组分析：

步骤S1：数据质控。对宏基因组测序所得原始数据进行过滤，并统计过滤前后碱基组成、reads数量、碱基质量。后续生物信息分析均基于过滤后的高质量数据进行；

步骤S2：多策略组装。提供4类组装策略，将reads组装为连续长序列(contig)，可根据样本、数据特征选择。组装完成后过滤掉长度在500bp以下的短序列(图1)，使用N50、N90、总长、平均长度等指标评估组装效果，其中以N50评估最为主要。后续分析基于500bp以上的contig进行。

策略1，若样本为肠道等非自然环境样本，当单样本测序数据量小于20G时，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127。

策略2，若样本为肠道等非自然环境样本，当单样本测序数据量超过20G时，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127。

策略3，若样本为土壤、水体等自然环境样本，当单样本测序数据量小于20G时，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141。

策略4，若样本为土壤、水体等自然环境样本，当单样本测序数据量超过20G时，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141。

步骤S3：基因预测。将S2步骤获得的数据结果进行基因预测分析，包括基因预测、基因丰度统计、核心基因(Core gene)和泛基因(Pan gene)分析，具体内容为：

S3.1：基因预测：利用MetaGeneMark软件对contigs进行基因预测，获得基因序列，然后采用CD-HIT软件，将基因序列相似度大于95％，序列比对区域大于90％的基因聚类为一个cluster，95％identity、90％coverage，选取最长的基因作为每个cluster的代表序列，称为unigene。所有代表序列即为获得的非冗余基因集合(Unigene catalogue)。

S3.2：基因丰度统计。利用bowtie2将高质量数据的reads重新比对unigene，计算各样本中，每个基因的reads比对数目。过滤掉在各个样品中reads支持数目都≤2的基因，获得用于后续步骤分析的基因集合。基于比对上的reads数目及基因长度出发，按公式计算得到各基因在各样品中的丰度信息，r为比对上某基因k的reads数目，L为基因k的长度。

S3.3：构建和绘制Core基因和Pan基因稀释曲线。

基于各样品的基因丰度表，从所有样本中，随机抽取1个样本，统计基因数目，多次随机抽取后，获得1个样本时的基因数目分布盒型图；再从所有样本中随机抽取2个样本，统计两样本基因交集和并集的基因数量，多次随机抽取后，获得2个样本时的基因数目分布盒型图；重复抽取到10个样本时，最终获得不同数目样品组合时的基因数目分布，绘制了Core基因和Pan基因稀释曲线。

步骤S4：功能注释。将基因丰度表中的Unigenes比对到各类数据库进行基因功能注释和功能丰度统计。

具体操作为将Unigenes序列通过DIAMOND软件(阈值evalue<＝1e-5)比对到数据库，同时结合基因丰度表格计算不同数据库比对结果的功能丰度信息，并使用柱形图、circos图、热图进行可视化，以进行系统丰富的组间功能差异分析和比较。也可以使用其他的可视化图。

优选的，数据库包括但不限于KEGG、eggNOG、CAZy、CARD、VFDB、PHI中的一个或多个，基于各种数据库比对注释结果分析预测样本中微生物群落的功能特征。

步骤S5：多策略物种注释。提供两种物种分析策略，可根据分析需求选择。

策略S5-1，基于reads，使用Kaiju软件进行物种注释。首先将高质量reads翻译为氨基酸序列，从终止子处断开，然后挑选高得分氨基酸序列比对NCBI的Refseq数据库，获得序列的物种分类注释，然后基于每个物种的reads支持数目，获得不同分类水平下物种在各样本中的丰度表，并图示化展示。

策略S5-2，基于基因序列，使用DIAMOND软件比对至Nr数据库，按照LCA算法，获得基因的物种注释信息，并结合基因丰度表，统计不同分类水平下物种在各样本中的丰度表，并图示化展示。

优选的，图示化展示包括但不限于物种分布堆叠图、物种分布热图、物种分布Circos图。

步骤S6：进行组间多元统计比较，包括但不限于：样本相关性分析、PCA主成分分析、PCoA主坐标分析、NMDS分析、UPGMA分类树分析、Anosim差异分析、Adonis分析。

步骤S7：进行差异分析，包括但不限于：welch’s T检验、方差分析(ANOVA，Analysis of Variance)、三元图分析、MetaStats差异物种分析、LefSe差异物种分析。

binning分析

该分析流程包括如下步骤：

步骤S8：去嵌合体校正。对步骤S2所获得的组装结果(contig)，将样本测序的高质量reads比对contig，计算一条contig上两个窗口的reads覆盖度，基于覆盖度差异程度判断是否为嵌合体，并进行切分校正。

步骤S9：binning。选择长度大于1.5kb的contig，将样本测序的高质量reads再次比对contig，计算每条contig的覆盖度和GC含量，使用metaBAT2进行binning，即进行contig聚类，软件鉴定判断后，将具有相似丰度和核酸组成的contig归为一个bin，即一个潜在的基因组。

步骤S10：bin质量评估。使用CheckM软件，基于单拷贝基因集，评估每个bin的完整度、污染度。基于contig的reads比对结果，统计每个bin的contig数量、测序深度、长度、GC特征等。默认按完整度高于80％且污染度低于5％的阈值，筛选高质量的bin。

使用GC-depth气泡图，直观呈现每个样本中高质量bin整体的聚类效果，以点大小展示contig长度，以点颜色展示contig对应的bin分类。

使用GC-depth分布图，详细展示每个bin的contig分布特征，以反映bin的污染程度、bin大小等信息。

步骤S11：多策略高质量bin筛选

汇总文献中有参考价值的多种筛选阈值的组合，便于客户根据自己的研究方向，灵活选择筛选策略，提高bin数据的使用效率和研究深度。

1)研究方向为针对目标bin开展基因组分析，即需要精选目标bin范围，尽可能获得高质量的有价值的基因组，可选完整度高于80％～90％且污染度低于5％～10％的bin。

2)研究方向为获得宏基因组、宏转录组分析的参考基因组，则对基因组的质量要求较低，仅需要获得高维度基因(宏基因组中上百万的基因)的分类，即需要尽可能多的参考物种数量，可选完整度高于50％～80％且污染度低于10％～40％的bin。

统计各样本所有高质量bin的contig组成信息，并从每个样本中挑选完整度排名top20的bin，绘制binning示意图，直观展示一个样本中高质量bin的特征(图2)。图中，横轴表示GC含量，纵轴表示coverage(测序深度)。图中点表示contig，点颜色表示contig所属的bin，大小表示contig长度。图形可呈现样本中高质量bin的聚类特征。

步骤S12：所有高质量bin物种注释：

步骤S12.1：使用CheckM软件预测获得每个bin的基因序列信息，并使用DIAMOND软件将基因序列比对Nr库，提取基因在种水平的物种注释，并基于物种分类层级关系，统计获得基因对应的门、纲、目、科、属等各分类水平的物种注释。

步骤S12.2：统计一个bin中注释为相同种的基因总长度占bin基因总长度的比例，输出占比大于0.1的物种名称，作为该bin门到种水平的物种注释。

步骤S12.3：统计一个bin中注释为相同属的基因总长度占bin基因总长度的比例，输出占比大于0.1的物种名称，作为该bin门到属水平的物种注释。

步骤S13：高质量bin基因注释。基于基因序列，使用DIAMOND软件将基因序列比对Nr库，提取基因的功能描述信息。

使用气泡图展示每个高质量bin的基因特征，即以横轴展示基因长度、以纵轴展示基因序列比对的相似度、以点大小展示基因序列比对的得分、以点颜色展示基因在种水平的物种注释(本发明特有)，点大小表示数据库比对的score得分，点越大，表示注释越可靠。从图中不同颜色点的分布和大小，可初步了解每个bin注释到的物种数量、基因特征等(图3)，可综合呈现一个bin的物种注释情况。

bin后续分析

A.目标bin基因组分析

综合考虑上述分析获得的bin大小、测序深度、质量、物种注释、基因注释等信息，由用户挑选感兴趣的一个或多个bin，提供将bin作为“单菌基因组”继续深入分析的功能，并且结合bin的特征丰富调整分析细节。

为填补针对目标bin后续数据分析的空白，本发明经过调查汇总，确定了多种目标bin的分析策略。该分析流程包括如下步骤：

步骤S14：目标bin基因组圈图。提供两种策略，可根据bin数据情况选择。

策略S14-1，原始contig。考虑到bin的特殊性(纯分析手段预测，无实验分离培养)故提供新的细菌基因组获得方法。使用bin的contig，即宏基因组拼接获得的contig，在完整度允许条件下，比对查找近缘物种，基于近缘物种特征，绘制基因组圈图(bin基因组分析特有)。

策略S14-2，序列重组装。将样本测序的高质量reads再次比对回该bin的contig，提取比对上该bin的reads，使用soapdenovo软件进行重组装，基于新的组装结果绘制基因组圈图(bin基因组特有，与传统二代测序细菌基因组分析直接使用测序reads组装不同，本流程需要将样本所有reads比回bin的每个contig，数据量巨大，并且基于物种新颖程度，设定了宽、松两种比对参数，获得属于该bin的reads，然后进行组装)。

步骤S15：基因预测。

策略S15-1，直接提取S12.1中目标bin基因序列开展分析。

策略S15-2，使用MetaGeneMark软件对bin基因组序列重新进行基因预测。

步骤S16：比较基因组分析。根据目标bin基因组间或者和数据库中挑选的近缘物种间，基于基因组序列，进行共线性分析、core-pan基因分析、基因家族分析、系统发育树分析等。

步骤S17：bin基因结构分析。根据目标bin中的contig序列或新组装的基因组序列，分析目标基因上下游序列特征、基因信息等。

步骤S18：bin基因功能特征分析。基因序列可比对Swiss-Prot、GO、KEGG、COG、Pfam、PHI、CAZy、CARD、VFDB数据库进行功能注释等分析。

B.高质量bin分析。以宏转录组结合为例，bin基因表达量鉴定：

将每个bin作为一个物种，即一个参考基因组，基于bin的基因注释，统计汇总每个bin的基因、KEGG等不同功能分类层级的功能表达量，并统计bin在样本所有基因表达量的均值，使用该均值校正bin基因表达量，即每个基因的表达量比上均值，作为最终bin中每个基因的表达量，以排除样本中物种数量高低的影响，便于科学进行物种间基因表达活性的对比。可基于表达量特征，开展群落中特殊单个物种的功能表达规律探究。

下面结合实施例，进一步说明本发明的技术方案。

在行业内，N50可作为基因组拼接效果好坏的一个判断标准，且N50越大表明拼接效果越好。

实施例1(组装策略2)

非自然环境样本，不同人肠道4例样本(总数据量约为6Gb)的组装情况，对肠道样本目前现有技术较好的组装效果为文献Han M,Yang P,Zhong C,et al.The Human GutVirome in Hypertension[J].Frontiers in Microbiology,2018,9:中Assembly of theHuman Gut Metagenomic Data分析部分，文献中N50为4152bp，本实施例和参考文献采用kmer组装参数为27，37，47，57，67，77，87，97，107，117，127，采用本实施例组装的N50均高于现有技术，组装效果更好。

	数量	总长度	平均长度	N50	N90	最大长度	最小长度	GC
									样本1	27273	87611145	3212.38	4785	1263	371303	1000	48.54％
样本2	26592	88292244	3320.26	4967	1329	124147	1000	57.51％
									样本3	17137	58929791	3438.75	5418	1315	277887	1000	49.90％
样本4	22432	74417541	3317.47	5009	1320	190304	1000	51.53％

注：因文献中未公开具体的质控细节，且是所有样本的组装结果，与本实施例的数据处理不完全相同。

实施例二(组装策略4)：

与现有技术文献(Zhang M,Pan L,Huang F,et al.Metagenomic analysis ofcomposition,function and cycling processes of microbial community in water,sediment and effluent of Litopenaeus vannamei farming environments underdifferent culture modes[J].Aquaculture,2019,506:280-293.)中3例自然环境样本即水体样本(GW，HE，HW)数据，以本发明技术方法流程，现有参数，和文献数据比对，结果如下：

(1)GW水体，在文献中记载总长度为288496，N50为954，将样本数据以现有技术参数和本发明技术进行比对，结果数据如下，本发明技术N50为1027，组装效果强于现有技术参数。

(2)HE水体，在文献中记载总长度为356232，N50为914，将样本数据以现有技术参数和本发明技术进行比对，结果数据如下，本发明技术N50为919，组装效果强于现有技术参数。

(3)HW水体，在文献中记载总长度为274626，N50为1122，将样本数据以现有技术参数和本发明技术进行比对，结果数据如下，本发明技术N50为1355，组装效果强于现有技术参数。

实施例三(组装策略3)

自然环境样本-沉积物样本3例，单例数据量为6G，组内样本reads混合组装，混合组装效果较好。

	数量	总长度	平均长度	N50	N90	最大长度	最小长度	GC
									样本1	177935	289272894	1625.72	2560	618	755937	500	62.61％
样本2	177055	287062202	1621.32	2559	617	755937	500	62.63％
									样本3	174799	283115279	1619.66	2579	616	755937	500	62.64％
混合组装	325995	578768707	1775.39	3070	646	1275998	500	62.52％

Claims (10)

1.一种微生物群宏基因组的分析方法，包括：

测序数据过滤：对宏基因组测序所得原始数据进行过滤，并统计过滤前后碱基组成、reads数量、碱基质量，得到高质量测序数据；

数据组装：根据样本来源和测序数据量的大小，选择不同的数据组装策略，得到组装结果contigs，优选的，组装完成后过滤掉长度在500bp以下的短序列，其中：

策略1：若样本为非自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略2：若样本为非自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略3：若样本为自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

策略4：若样本为自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

基因数据分析：将S2步骤获得的数据结果进行基因数据分析，分析内容选自基因预测、基因丰度统计、核心基因和泛基因分析；优选的，基因预测的操作包括：利用MetaGeneMark软件对contigs进行基因预测，获得基因序列，然后采用CD-HIT软件，将基因序列相似度大于95％，序列比对区域大于90％的基因聚类为一个cluster，95％identity、90％coverage，选取最长的基因作为每个cluster的代表序列，称为Unigene，所有代表序列即为获得的非冗余基因集合Unigene catalogue；

基因丰度统计的操作包括：利用bowtie2将高质量数据的reads重新比对Unigene，计算各样本中，每个基因的reads比对数目；过滤掉在各个样品中reads支持数目都≤2的基因，获得用于后续步骤分析的基因集合；基于比对上的reads数目及基因长度出发，按公式

计算得到各基因在各样品中的丰度信息，式中，r为比对上某基因k的reads数目，L为基因k的长度；

核心基因和泛基因分析的操作包括：基于各样品的基因丰度表，从所有样本中，随机抽取1个样本，统计基因数目，多次随机抽取后，获得1个样本时的基因数目分布盒型图；

再从所有样本中随机抽取2个样本，统计两样本基因交集和并集的基因数量，多次随机抽取后，获得2个样本时的基因数目分布盒型图；重复抽取足够样本，最终获得不同数目样品组合时的基因数目分布，得到核心基因和泛基因稀释曲线；

优选的，非自然环境样本为来自动物体的微生物群；所述自然环境样本为来自土壤、水体、空气的样本。

2.根据权利要求1所述的分析方法，其特征在于：还包括功能注释：基于基因数据分析结果，将Unigenes比对到各类数据库进行基因功能注释和功能丰度统计；

优选的，功能注释具体为将Unigenes序列通过DIAMOND软件，阈值evalue<＝1e-5，比对到数据库，同时结合基因丰度表格计算不同数据库比对结果的功能丰度信息。

3.微生物群宏基因组数据的binning分析方法，包括：

按权利要求1的操作得到组装结果contigs；

去嵌合体校正：将样本测序的reads比对contig，计算一条contig上两个窗口的reads覆盖度，基于覆盖度差异程度判断是否为嵌合体，并进行切分校正；

binning：选择长度大于1.5kb的contig，将样本测序的高质量reads再次比对contig，计算每条contig的覆盖度和GC含量，使用metaBAT2进行binning，即进行contig聚类，软件鉴定判断后，将具有相似丰度和核酸组成的contig归为一个bin，即一个潜在的基因组；或

进一步包括高质量bin筛选，具体包括：基于单拷贝基因集，评估每个bin的完整度、污染度；优选的，使用CheckM软件计算bin的完整度、污染度；完整度高于80％且污染度低于5％归类为高质量bin。

4.根据权利要求3所述的binning分析方法，其特征在于：基于contig的reads比对结果，统计每个bin的contig数量、测序深度、长度、GC特征并进行可视化分析。

5.根据权利要求3所述的binning分析方法，其特征在于：还包括：

高质量bin物种注释，包括使用CheckM软件预测获得每个bin的基因序列信息，并使用DIAMOND软件将基因序列比对Nr库，提取基因在种水平的物种注释，并基于物种分类层级关系，统计获得基因对应的各分类水平的物种注释；或

高质量bin基因注释，包括基于基因序列，使用DIAMOND软件将基因序列比对Nr库，提取基因的功能描述信息。

6.根据权利要求3所述的binning分析方法，其特征在于：还包括：

使用气泡图展示每个bin的基因特征，即以横轴展示基因长度、以纵轴展示基因序列比对的相似度、以点大小展示基因序列比对的得分、以点颜色展示基因在种水平的物种注释，点大小表示数据库比对的score得分，点越大，表示注释越可靠；和/或

选择目标bin构建bin基因组圈图，构建操作包括：

策略B1：使用bin的contig，即宏基因组拼接获得的contig，在完整度允许条件下，比对查找近缘物种，基于近缘物种特征，绘制基因组圈图；或

策略B2：将样本测序的reads再次比对回该bin的contig，提取比对上该bin的reads，使用soapdenovo软件进行重组装，基于新的组装结果绘制基因组圈图；优选的，样本测序的reads为其高质量reads。

7.一种微生物群宏基因组的分析系统，包括：

数据存储装置，用于存储待分析的高通量测序数据；

数据处理装置，用于对高通量测序数据进行组装和分析，其中，数据组装包括：

根据样本来源和测序数据量的大小，选择不同的数据组装策略，得到组装结果contigs，优选的，组装完成后过滤掉长度在500bp以下的短序列，其中：

策略1：若样本为非自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略2：若样本为非自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略3：若样本为自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

策略4：若样本为自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

数据分析包括：

基因数据分析：将S2步骤获得的数据结果进行基因数据分析，分析内容选自基因预测、基因丰度统计、核心基因和泛基因分析；优选的，基因预测的操作包括：利用MetaGeneMark软件对contigs进行基因预测，获得基因序列，然后采用CD-HIT软件，将基因序列相似度大于95％，序列比对区域大于90％的基因聚类为一个cluster，95％identity、90％coverage，选取最长的基因作为每个cluster的代表序列，称为Unigene，所有代表序列即为获得的非冗余基因集合Unigene catalogue；

基因丰度统计的操作包括：利用bowtie2将高质量数据的reads重新比对Unigene，计算各样本中，每个基因的reads比对数目；过滤掉在各个样品中reads支持数目都≤2的基因，获得用于后续步骤分析的基因集合；基于比对上的reads数目及基因长度出发，按公式

计算得到各基因在各样品中的丰度信息，式中，r为比对上某基因k的reads数目，L为基因k的长度；

核心基因和泛基因分析的操作包括：基于各样品的基因丰度表，从所有样本中，随机抽取1个样本，统计基因数目，多次随机抽取后，获得1个样本时的基因数目分布盒型图；

再从所有样本中随机抽取2个样本，统计两样本基因交集和并集的基因数量，多次随机抽取后，获得2个样本时的基因数目分布盒型图；重复抽取足够样本，最终获得不同数目样品组合时的基因数目分布，得到核心基因和泛基因稀释曲线；

结果输出装置，用于输出处理后的数据。

8.根据权利要求7所述的分析系统，其特征在于：所述数据分析还包括功能注释：基于基因数据分析结果，将Unigenes比对到各类数据库进行基因功能注释和功能丰度统计；

优选的，功能注释具体为将Unigenes序列通过DIAMOND软件，阈值evalue<＝1e-5，比对到数据库，同时结合基因丰度表格计算不同数据库比对结果的功能丰度信息。

9.一种微生物群宏基因组数据的binning分析系统，包括：

数据存储装置，用于存储待分析的高通量测序数据；

数据处理装置，用于对高通量测序数据进行组装和分析，其中，数据组装包括：

根据样本来源和测序数据量的大小，选择不同的数据组装策略，得到组装结果contigs，优选的，组装完成后过滤掉长度在500bp以下的短序列，其中：

策略1：若样本为非自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略2：若样本为非自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为27、37、47、57、67、77、87、97、107、117、127；

策略3：若样本为自然环境样本，单样本测序数据量小于20G，将组内样本reads混合组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

策略4：若样本为自然环境样本，单样本测序数据量超过20G，单样本单独组装；使用MEGAHIT软件，组装参数设置为：kmer长度为21、41、61、81、101、121、141；

数据分析包括：

去嵌合体校正：将样本测序的reads比对contig，计算一条contig上两个窗口的reads覆盖度，基于覆盖度差异程度判断是否为嵌合体，并进行切分校正；

binning：选择长度大于1.5kb的contig，将样本测序的高质量reads再次比对contig，计算每条contig的覆盖度和GC含量，使用metaBAT2进行binning，即进行contig聚类，软件鉴定判断后，将具有相似丰度和核酸组成的contig归为一个bin，即一个潜在的基因组；或进一步包括

高质量bin筛选，具体包括：基于单拷贝基因集，评估每个bin的完整度、污染度；优选的，使用CheckM软件计算bin的完整度、污染度；完整度高于80％且污染度低于5％归类为高质量bin；或包括

高质量bin物种注释，包括使用CheckM软件预测获得每个bin的基因序列信息，并使用DIAMOND软件将基因序列比对Nr库，提取基因在种水平的物种注释，并基于物种分类层级关系，统计获得基因对应的各分类水平的物种注释；或

高质量bin基因注释，包括基于基因序列，使用DIAMOND软件将基因序列比对Nr库，提取基因的功能描述信息；

结果输出装置，用于输出处理后的数据。

10.根据权利要求9所述的binning分析系统，其特征在于：所述数据分析还包括基于contig的reads比对结果，统计每个bin的contig数量、测序深度、长度、GC特征并进行可视化分析；

优选的使用气泡图展示每个bin的基因特征，即以横轴展示基因长度、以纵轴展示基因序列比对的相似度、以点大小展示基因序列比对的得分、以点颜色展示基因在种水平的物种注释，点大小表示数据库比对的score得分，点越大，表示注释越可靠；和/或

选择目标bin构建bin基因组圈图，构建操作包括：

策略B1：使用bin的contig，即宏基因组拼接获得的contig，在完整度允许条件下，比对查找近缘物种，基于近缘物种特征，绘制基因组圈图；或

策略B2：将样本测序的reads再次比对回该bin的contig，提取比对上该bin的reads，使用soapdenovo软件进行重组装，基于新的组装结果绘制基因组圈图；优选的，样本测序的reads为其高质量reads。

Images