一种用于宏基因组测序数据的微生物物种与功能组成分析方法

技术领域

本发明涉及微生物基因分析领域，尤其涉及一种用于宏基因组测序数据的微生物物种与功能组成分析方法。

背景技术

随着下一代高通量测序技术(Next Generation Sequencing，NGS)的不断发展，人们对于微生物群落方面的研究也越来越全面和深入。区别于常见的靶向微生物核糖体RNA基因的扩增子测序技术，宏基因组学是将整个群落系统中全体微生物的基因组作为研究对象，基于鸟枪法测序技术，全面展现整个群落的物种组成和功能潜能组成，进而阐明微生物群落的作用机制。然而，由于测序样本类型的多样，样本量规模、测序深度的多变，以及宿主、污染基因组的多少等因素，以及分析本身的复杂性，研究者们发开出了数量庞大的各类软件，以及配套的更为复杂的各种分析模式、参数和数据库。目前，只有少数几个流程软件或分析网站提供了自动化的宏基因组分析方法或服务。其中，以MG-RAST(https://www.mg-rast.org/)和IMG/M(https://img.jgi.doe.gov/cgi-bin/m/main.cgi)等(宏)基因组数据库网站为代表，它们在提供数据上传和存储服务的同时，也提供了有限种类的分析项目和计算资源。而另一类软件，如HUMAnN2(http://www.huttenhower.org/humann2)以及发明专利“一种宏基因组数据分析方法及系统”(CN 108334750 B)等，则是基于非拼接序列数据的数据库比对注释，也就是说，其物种或基因的鉴定，主要依赖于数据库对目标样本的微生物基因组的覆盖度。然而，对于环境样本，如土壤，底泥，极端环境，我们对于微生物群落的认识是极其有限的，因而此类分析流程，适用的样本类型也是有限，主要以人体相关样本为宜。

这些现有的宏基因组分析流程具有如下缺陷：

(1)未包含拼接步骤，无法满足基于拼接产生的叠连群序列的进一步分析。

(2)物种分类学注释依靠非拼接序列数据的比对分析，未包含序列拼接并基于拼接后的数据进行相应鉴定，因而结果假阳性高，依赖于专用数据库的构建，适用面有限。

(3)对于样本中基因的鉴定，主要采用比对方式，依赖于专用数据库的构建，假阳性高，灵敏度低，适用范围有限。

(4)未包含宿主和/或污染基因组序列去除步骤，或去除步骤依赖于宿主和/或污染基因组的存在。

发明内容

本发明的提供一种用于宏基因组测序数据的微生物物种与功能组成分析方法。

本发明的方案是：

一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：

1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；

2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；

3)对剔除非目物种后的序列进行拼接，获得叠连群序列；

4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；

5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；

6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；

7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；

8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

作为优选的技术方案，步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或trimmomatic软件其中一种将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2或bmtagger其中一种软件，去除比对到宿主基因组上的序列。

序列质控：检查原始宏基因组数据的测序质量；使用移动划窗法将原始数据中低质量片段切除，并弃去过短序列；若已知宿主基因组，则将比对到宿主基因组上的序列剔除。

作为优选的技术方案，使用fastp软件将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；

作为优选的技术方案，使用bmtagger软件去除比对到宿主基因组上的序列。

作为优选的技术方案，步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k-mer检索或使用kaiju及蛋白序列数据库进行局部相似性比对，对序列数据进行物种注释；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列或者注释到自定义的物种的序列。

非拼接序列物种注释：进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一。

作为优选的技术方案，使用kraken2及核酸序列数据库，对质控后的序列数据进行物种注释，并依据注释结果剔除非目的序列，使用bracken计算物种组成丰度表。

作为优选的技术方案，步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或先对各样本按分组分别合并拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或对所有样本合并拼接；所述拼接用的拼接软件为megahit或metaspades；使用minimap2或bowtie2软件找出各样本无法比对上叠连群的序列。

序列拼接：对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余。

作为优选的技术方案，所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2或cd-hit，所述MMseqs2采用easy-linclust模式；计算各样本非冗余的叠连群序列丰度采用基于比对的htseq-count或不基于比对的salmon。

叠连群序列丰度计算：计算叠连群丰度，并去除丰度为0的序列，对余下叠连群序列进行相似性聚类，获得非冗余的叠连群序列的丰度表。

作为优选的技术方案，叠连群核酸序列聚类软件为MMseqs2。

作为优选的技术方案，计算叠连群丰度采用salmon。

作为优选的技术方案，所述步骤5)中比对用软件为minimap2或bowtie2。

叠连群序列物种注释：使用blastn算法，将非冗余的叠连群序列与NCBI-nt数据库进行比对，采用最近祖先算法获取叠连群序列的物种信息；再结合非冗余的叠连群序列的丰度表，获得物种组成丰度表二；同时基于物种信息，将非拼接序列物种注释模块中注释得到的物种分为已验证存在的物种和疑似存在物种。

作为优选的技术方案，所述步骤6)中基因预测软件为metagenemark、prodigal或FragGeneScan其中一种；蛋白序列聚类软件为cd-hit或MMseqs2，所述MMseqs2的聚类模式为easy-linclust或easy-cluster模式。

基因预测：预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集。

作为优选的技术方案，蛋白序列聚类软件为MMseqs2。

作为优选的技术方案，MMseqs2的聚类模式为easy-cluster模式。

作为优选的技术方案，步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI-nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2或diamond软件。

功能注释：将非冗余蛋白序列集与蛋白或功能数据库进行比对，获得蛋白序列的功能注释信息，也即是基因序列的功能注释信息。

作为优选的技术方案，使用MMseqs2将非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对。

作为优选的技术方案，所述步骤8)中计算基因序列丰度采用基于比对的htseq-count或不基于比对的salmon，所述比对用软件为minimap2或bowtie2。

基因丰度计算：计算基因序列丰度，也即是冗余蛋白序列的丰度，进而通过蛋白序列聚类的对应表，获得非冗余蛋白序列的丰度表，再结合这些蛋白的功能注释信息，最终获得功能单元的丰度表。

作为优选的技术方案，使用salmon计算基因序列丰度。

作为优选的技术方案，基于KEGG数据注释获得的KO丰度表，使用MinPath推断存在的KEGG代谢通路并计算其丰度表。

本发明提供了一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；3)对剔除非目物种后的序列进行拼接，获得叠连群序列；4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

本发明的优点：

满足了拼接产生的叠连群序列的进一步分析，避免了结果假阳性高，依赖于专用数据库的构建，适用面窄的问题；灵敏度高，提高了准确性；

将复杂的宏基因组分析流程拆解总结为八大模块，再基于这些模块而非软件搭建流程，允许在特定模块中，选择调用不同的分析策略、软件或参数以应对各类情形。

对于物种的分类学鉴定和丰度定量，同时包含了：基于非拼接序列的比对或k-mer检索，用于快速鉴定物种并准确定量；基于拼接叠连群序列的blastn比对，用于严格判别丰度较高的物种或关注的物种的有无，降低物种鉴定假阳性。

对于功能(潜能)的鉴定和丰度定量，同时包含了：基于非拼接序列的比对或k-mer检索，用于快速鉴定物种并准确定量；基于拼接叠连群序列的blastn比对，用于严格判别丰度较高的物种或关注的物种的有无，降低物种鉴定假阳性。

输出的分析结果文件全面，使用者可依据不同需求进行相应选择，为后续个性化的统计分析提供可能。

附图说明

图1为本发明实施例2的流程示意图。

具体实施方式

为了弥补以上不足，本发明提供了一种用于宏基因组测序数据的微生物物种与功能组成分析方法以解决上述背景技术中的问题。

一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：

1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；

2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；

3)对剔除非目物种后的序列进行拼接，获得叠连群序列；

4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；

5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；

6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；

7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；

8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或trimmomatic软件其中一种将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2或bmtagger其中一种软件，去除比对到宿主基因组上的序列。

步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k-mer检索或使用kaiju及蛋白序列数据库进行局部相似性比对；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列或者注释到自定义的物种的序列。

步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或先对各样本按分组分别合并拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；或对所有样本合并拼接；所述拼接用的拼接软件为megahit或metaspades；使用minimap2或bowtie2软件找出各样本无法比对上叠连群的序列。

所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2或cd-hit，所述MMseqs2采用easy-linclust模式；计算各样本非冗余的叠连群序列丰度采用基于比对的htseq-count或不基于比对的salmon。

所述步骤5)中比对用软件为minimap2或bowtie2。

所述步骤6)中基因预测软件为metagenemark、prodigal或FragGeneScan其中一种；蛋白序列聚类软件为cd-hit或MMseqs2，所述MMseqs2的聚类模式为easy-linclust或easy-cluster模式。

步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI-nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2或diamond软件。

所述步骤8)中计算基因序列丰度采用基于比对的htseq-count或不基于比对的salmon，所述比对用软件为minimap2或bowtie2。

为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。

实施例1：

一种用于宏基因组测序数据的微生物物种与功能组成分析方法，包括下列步骤：

1)切除原始数据中的接头序列片段、低质量片段，滤除过短序列、含模糊碱基序列；若已知宿主基因组，则将宿主序列剔除；

2)使用上述获得的数据进行物种注释，并统计物种序列数即为丰度，再基于注释结果剔除注释到非目的物种的序列；

3)对剔除非目物种后的序列进行拼接，获得叠连群序列；

4)对叠连群序列进行相似性聚类，并计算各样本非冗余的叠连群序列丰度，并去掉总丰度为零的序列；

5)使用blastn算法，对非冗余的叠连群序列进行核酸数据库比对，并采用共同祖先算法获取拼接序列的物种注释信息，再基于注释信息，将步骤2)中注释得到的物种分为已验证存在的物种与疑似存在物种；

6)预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集；

7)将非冗余蛋白序列集与各类蛋白注释数据库进行比对，获得基因序列和蛋白序列的功能注释信息；

8)计算基因序列丰度，再通过基因/蛋白对应功能信息，获得功能丰度表。

步骤1)中使用FastQC检查原始数据的测序质量情况；使用fastp或将原始数据中低质量片段切割，并滤除接头序列和过短序列，获得高质量序列；使用bowtie2软件，去除比对到宿主基因组上的序列。

步骤2)中进行物种注释为进行非拼接序列的物种注释，将质控后的序列基于物种注释数据库进行k-mer检索或局部相似性比对，获取序列的物种注释信息以及物种丰度表一，其中使用kraken2及核酸序列数据库进行k-mer检索；使用bracken计算物种组成丰度表；非目的序列默认为注释到后生动物及绿色植物的序列。

步骤3)中序列拼接为对剔除非目物种后的序列进行拼接，获得叠连群序列，并去冗余；所述拼接方式包括各样本单独拼接，再分别找出各样本无法比对上的序列，对无法比对上的序列进行混合拼接；所述拼接用的拼接软件为megahit；使用minimap软件找出各样本无法比对上叠连群的序列。

所述叠连群核酸序列行相似性聚类采用的聚类软件为MMseqs2，所述MMseqs2采用easy-linclust模式；计算各样本非冗余的叠连群序列丰度采用不基于比对的salmon。

所述步骤5)中比对用软件为minimap2。

所述步骤6)中基因预测软件为metagenemark；蛋白序列聚类软件为MMseqs2，所述MMseqs2的聚类模式为easy-linclust模式。

步骤7)中功能注释的蛋白或功能数据库包括eggNOG、KEGG、CAZy、GO、NCBI-nr、swissprot；所述非冗余蛋白序列集与功能数据库的参考蛋白序列集进行比对使用使用MMseqs2软件。

所述步骤8)中计算基因序列丰度采用不基于比对的salmon，所述比对用软件为minimap2。

实施例2：

在本发明实施例中，使用模拟宏基因组数据进行分析，模拟数据物种组成见附表1，其中已知宿主基因组为人基因组。

在步骤S101中，首先使用FastQC检查原始数据的测序质量情况；使用Cutadapt识别3’端潜在的接头序列，并在识别的接头处截断。这里。要求与接头序列(R1:AGATCGGAAGAGCACACGTCTGAACTCCAGTCA；R2:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT)的匹配长度至少达到3bp，且允许至多20％的碱基错配率。再使用fastp软件将其中的低质量片段切割，具体为采用滑动窗口法对序列进行质量筛查：窗口大小为5bp，从5'端第一个碱基位置开始移动，要求窗口中碱基平均质量大于等于Q20(即碱基平均测序准确率大于99％)，从第一个平均质量值低于Q20的窗口的3'端碱基处截断序列；经上述质量筛查后，去除序列长度小于50bp的序列，以及含有模糊碱基的序列。使用bmtagger软件去除比对到宿主——人基因组上的序列。

在步骤S102中，进行非拼接序列的物种注释，将步骤S101获得的有效序列基于物种注释数据库进行搜索或比对，获取序列的物种注释信息。这里使用kraken2及核酸序列数据库，对质控后的序列数据进行物种注释，并依据注释结果剔除非目的序列，使用bracken计算物种组成丰度表。非目的序列设置为注释到后生动物及绿色植物的序列。使用Kraken2进行物种注释。Kraken2通过使用精准k-mer匹配来实现高精度且快速的物种注释。数据库使用NCBI的RefSeq基因组数据库中的细菌、古菌、真菌、原生动物、病毒基因组，以及后生动物和绿色植物基因组数据作为参考构建数据库，调整置信度为0.5，进行注释分析。基于注释结果，再剔除注释到非目的物种的序列，这里的非目的物种为Metazoa(后生动物)和Viridiplantae(绿色植物)。最后，使用Bracken计算物种组成丰度表

在步骤S103中，先采用megahit的meta-sensitive模式，先对各样本单独拼接，再采用minimap2软件分别找出各样本无法比对上的序列，对这些序列再采用megahit的meta-sensitive模式进行混合拼接。

在步骤S104中，先使用软件MMseqs2的easy-linclust模式将(合并的)叠连群序列集按相似度95％，对齐区域覆盖度90％(占短序列的比例)进行聚类去冗余，获得非冗余的叠连群集合。丰度计算采用salmon软件，采用k-mer检索算法，计算冗余叠连群序列丰度，并依据聚类的对应关系，获得非冗余的叠连群序列的丰度表。

在步骤S105中，使用blastn，将非冗余的叠连群序列与ncbi-nt数据库比对，E值设置为0.00001，我们保留比对结果中E值最小的5个结果，采用Blast2lca软件中的“最近共同祖先(Lowest Common Ancestor,LCA)”算法，将参考序列分化为不同物种分枝前的最后一级共同分类，作为目标序列的物种分类注释信息。并进一步结合非冗余的叠连群序列的丰度表，获得物种组成丰度表二。基于该信息，再将步骤二中注释得到的物种分为已验证存在的物种和疑似存在物种。

在步骤S106中，进行基因预测。预测非冗余的叠连群序列中的基因区域，获得基因序列及其翻译的蛋白序列，再对蛋白序列进行相似性聚类，获得非冗余蛋白序列集。对非冗余的叠连群序列采用MetaGeneMark软件识别其中的开放阅读框，并预测其中的编码区域，从而获得对应的基因序列文件，蛋白序列文件，GTF(gene transfer format)文件以及GFF(general feature format)文件。再使用软件MMseqs2的easy-cluster模式将(合并的)蛋白序列集按相似度95％，对齐区域覆盖度90％(占短序列的比例)进行去冗余，获得非冗余的蛋白序列集。

在步骤S107中，采用MMseqs2的search模式，将非冗余蛋白序列集与eggnog数据库蛋白参考序列进行比对，获得蛋白序列的eggNOG Orthologous Groups功能注释信息。基于eggNOG注释信息，使用eggnog-mapper，推定蛋白序列的GO的词条，获得其GO功能注释信息。采用MMseqs2的search模式，将非冗余蛋白序列集与KOBAS软件的蛋白数据库进行比对，获得蛋白序列对应的KO信息。采用MMseqs2的search模式，将非冗余蛋白序列集与KOBAS软件的蛋白数据库进行比对，获得蛋白序列对应的KO信息。采用MMseqs2的search模式，将非冗余蛋白序列集与CAZy数据库的蛋白序列进行比对，获得蛋白序列的碳水化合物活性酶注释信息。采用MMseqs2的search模式，将非冗余蛋白序列集与NCBI-nr和swissprot数据库进行比对，获得蛋白序列的多种功能注释信息。

在步骤S108中，使用salmon计算基因序列丰度，也即是冗余的蛋白序列丰度；再通过聚类的对应表，获得非冗余蛋白序列的丰度表，结合S107得到的蛋白的功能信息，获得功能丰度表。基于各样本KO的丰度，使用MinPath推断各样本存在的KEGG代谢通路并计算其丰度。

自此，获得了基于分类水平注释的物种丰度表，物种分为已验证存在的物种和疑似存在物种。同时，也获得了eggNOG、KEGG、CAZy、GO等功能数据库注释的功能丰度表。

附表1如下：

以上显示和描述了本发明的基本原理、主要特征及本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。