CN111192630A

CN111192630A - 一种宏基因组数据挖掘方法

Info

Publication number: CN111192630A
Application number: CN201911343764.XA
Authority: CN
Inventors: 张俊亚; 魏源送
Original assignee: Research Center for Eco Environmental Sciences of CAS
Current assignee: Research Center for Eco Environmental Sciences of CAS
Priority date: 2019-12-24
Filing date: 2019-12-24
Publication date: 2020-05-22
Anticipated expiration: 2039-12-24
Also published as: CN111192630B

Abstract

本发明涉及一种宏基因组数据挖掘方法，包括以下步骤：1)从KEGG数据库中获取目标代谢途径的全部基因信息，并建立特异性数据库；2)建立目标代谢途径的特异性数据库的映像文件；3)基于获得的目标代谢途径的特异型数据库，对宏基因组测序获得的clean reads进行数据库快速比对，得到各样品的比对结果；4)对各样品的比对结果进行排序、统计和整合；5)对各样品的注释结果进行均一化处理，根据均一化结果进行不同样品间的定量分析。本发明能够快速建立指定代谢途径的特异性数据库用于后续分析，并能够对数据进行均一化和后处理，进而定量化比较不同样品中代谢途径相关基因差异，因而，可以广泛应用于宏基因组数据挖掘领域。

Description

一种宏基因组数据挖掘方法

技术领域

本发明属于生物信息学分析领域，特别是涉及一种宏基因组数据挖掘方法。

背景技术

宏基因组测序得到越来越广泛的应用，其数据的挖掘技术不断更新，在宏基因组的生物信息学分析过程中，数据库的使用是后续功能分析的根本。现阶段国内外针对宏基因组数据的分析缺乏特异性，针对特定领域的数据库建设并不完善，分析结果在不同样本之间无法进行定量或半定量分析。传统的分析方法多为：双末端测序→拼接成重叠群(contigs)→开放阅读框(Open reading frame,ORF)注释→数据分析。在这个过程中会损失大量的测序序列。如一般宏基因组双末端测序(5G数据)会得到5000万左右的读条(reads)数量，拼接后一般会获得25万左右的contigs(＞500bp)，而其中可用于ORF注释的contigs约15万左右。以抗生素抗性基因(antibiotic resistance genes,ARGs)的研究为例，最终注释成ARGs的contigs一般为600个左右，无法对不同样品间ARGs的丰度进行定量比较，偏重于定性分析。

现阶段以质控后读条(clean reads)进行直接比对，能够将测序结果充分利用，并且得到的数据量大，能够对样品之间差异进行定量比较。这种研究方法在ARGs的相关研究中已经得到了广泛使用、验证和证实。然而，限制这种方法使用的瓶颈是特异性数据库的建立、分析和使用。现阶段存在的生物信息学数据库的特点是大而冗余，如著名的nr数据库，涵盖了已知的所有功能序列信息；eggnog数据库涵盖了已知的蛋白质序列信息；kegg数据库涵盖了已知的代谢途径、酶功能及序列信息；Cazy是糖代谢中涉及到的功能序列；而特异性功能数据库尚缺乏，如甲烷代谢的数据库、丙酸代谢数据库等；这种特异性小型数据库特别适合于小领域的研究，追求的是精准，如CARD的ARGs数据库、Ncyc的氮循环数据库、VFDB的毒力因子数据库等；这些小领域研究适合的数据库往往存在于大数据库之中，但这种特异性小型数据库的建立，往往如大海捞针，收集起来尤其繁琐。

发明内容

针对上述问题，本发明的目的是提供一种宏基因组数据挖掘方法，该方法能够实现快速高效的特异性数据库的构建，并通过基于reads的结果注释、整合、均一化处理和统计分析，实现不同样本间数据的可量化比较。

为实现上述目的，本发明采取以下技术方案：一种宏基因组数据挖掘方法，其包括以下步骤：

1)从KEGG数据库中获取目标代谢途径的全部基因信息，并建立特异性数据库DB.fasta；

2)建立目标代谢途径的特异性数据库DB.fasta的映像文件；

3)基于获得的特异型数据库DB.fasta，对宏基因组测序获得的clean reads进行比对，得到各样品的比对结果；

4)根据得到的各样品的比对结果以及特异型数据库DB.fasta的映像文件，进行排序、统计和整合；

5)对各样品的比对结果进行均一化处理，根据均一化处理结果进行不同样品间的定量分析。

进一步的，所述步骤1)中，从KEGG数据库中获取目标代谢途径的全部基因信息，并建立特异性数据库DB.fasta的方法，包括以下步骤：

1.1)从KEGG数据库中选定目标代谢途径，获得目标代谢途径的map序列号，并将获得的map序列号保存到ko_ID.txt文件；

1.2)获取KEGG数据库上的物种分类信息；

1.3)从步骤1.1)得到的ko_ID.txt文件中对map序列号进行识别，获得目标代谢途径的所有核酸和氨基酸序列，并将得到的所有核酸和氨基酸序列信息保存在ko_pathway_information.txt文件中；

1.4)根据步骤1.2)得到的物种分类信息，从步骤1.3)得到的目标代谢途径的所有核酸和氨基酸序列中去掉真核生物的基因序列，并根据最终产生的基因编号，通过TBtools工具中的序列提取命令，获得细菌和古菌的序列，作为最终的目标代谢途径的特异性数据库DB.fasta。

进一步的，所述步骤2)中，建立目标代谢途径的特异性数据库DB.fasta的映像文件的方法，包括以下步骤：

2.1)从建立的目标代谢途径的特异性数据库DB.fasta中，获得数据库的索引文件DB.fasta.fai；所述数据库索引文件DB.fasta.fai中，第一列是特异性数据库中的基因名，第二列为该基因对应的氨基酸序列长度；

2.2)通过共有的基因名排序，将步骤2.1)得到的索引文件DB.fasta.fai以及步骤1.3)得到的序列信息文件ko_pathway_information.txt进行合并，形成目标代谢途径的特异性数据库的mapping文件DB.txt。

进一步的，所述步骤3)中，基于获得的目标代谢途径的特异型数据库DB.fasta，对宏基因组测序获得的clean reads进行数据库比对，得到各样品的比对结果的方法，包括以下步骤：

3.1)基于获得的目标代谢途径的特异型数据库DB.fasta，构建适用于diamond软件的数据库；

3.2)基于构建的适用于diamond软件的数据库，对宏基因组测序获得的cleanreads序列文件进行比对，得到比对结果。

进一步的，所述步骤4)中，根据得到的各样品的比对结果以及特异型数据库DB.fasta的映像文件，进行排序、统计和整合的方法，包括以下步骤：

4.1)根据步骤2)中获得的特异型数据库的mapping文件DB.txt第一列中的基因名进行排序，并逐一统计计算各样品中每个基因比对到的reads个数，获得单个样品的数据库注释信息；

4.2)将多个样品的比对结果按样品名进行分列排序；

4.3)根据基因名排序同mapping文件进行合并，获得包含所有样品的完整基因定量化注释信息。

进一步的，所述步骤5)中，对各样品的比对结果进行均一化处理，根据均一化处理结果进行不同样品间的定量分析的方法，包括以下步骤：

5.1)对不同样品的比对结果进行均一化处理；

5.2)按照目标代谢途径基因对应的酶编号，对均一化后的注释信息做进一步统计加和，用于酶水平上基因的定量分析。

进一步的，所述步骤5.1)中，对不同样品的比对结果进行均一化处理时，所采用的均一化公式为：

式中，N_{targetgene-likesequence}为目标同源基因的数量；L_{referencesequence}为数据库中参考基因长度；L_reads宏基因组测序获得reads长度；N_{16S sequence}为宏基因组中16s rRNA的reads数量；L_{16S sequence}为16s rRNA比对数据库中的平均长度。

本发明由于采取以上技术方案，其具有以下优点：1、本发明通过自主研发的计算机语言命令从已有KEGG大数据库中获得目标代谢途径的特异性数据库，用于后续代谢途径的特异性分析。2、本发明根据所获得的特异性数据库，基于宏基因组测序结果所获得cleanreads进行比对，无需组装，大幅降低计算机运算成本。3、本发明可以用于定量比较不同样本间的基因差异。4、本发明通过数据的均一化处理，适应于不同测序深度下目标代谢途径基因的结果比较。5、本发明通过数据库构建→reads水平注释→均一化等步骤，实现了基于宏基因组测序结果目标代谢途径的个性化、可量化、可重复的分析。因此，本发明可以广泛应用于宏基因组数据挖掘领域。

附图说明

图1a和图1b为本发明实施例1甲烷代谢分析结果，其中，图1a为甲烷代谢途径中检测到的基因酶编号，图1b为对应酶水平上基因的丰度变化；

图2a和图2b为本发明实施例2氮代谢分析结果，其中，图2a为氮代谢途径中检测到的基因酶编号，图2b为对应酶水平上基因的丰度变化；

图3a和图3b为本发明实施例3硫代谢分析结果，其中，图3a为硫代谢途径中检测到的基因酶编号，图3b为对应酶水平上基因的丰度变化。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明提供的一种宏基因组数据挖掘方法，可用于个性化定制数据库，并针对特异代谢途径进行定量化解析。具体的，包括以下步骤：

1)通过自主开发的计算机语言命令个性化获得目标代谢途径的全部基因信息，并建立特异性数据库DB.fasta；

2)建立目标代谢途径的特异性数据库DB.fasta的映像(mapping)文件；

3)基于建立的目标代谢途径的特异型数据库DB.fasta，对通过宏基因组测序获得的clean reads进行数据库快速比对，得到各样品的比对结果；

上述步骤1)中，根据目标代谢途径建立特异型数据库的方法，包括以下步骤：

1.1)从已有KEGG数据库中选定目标代谢途径，获得目标代谢途径的map序列号(图片编号)，并将获得的map序列号保存到ko_ID.txt文件。

1.2)通过运行自主开发的计算机语言命令(spec_extract.pl)，获得KEGG数据库上的物种分类信息。

1.3)从步骤1.1)得到的ko_ID.txt文件中对map序列号进行识别，获得目标代谢途径的所有核酸和氨基酸序列，并将得到的所有核酸和氨基酸序列信息保存在ko_pathway_information.txt文件中。运行产生的tmp和tmp_seq建议不要删除，以后再次运行就不用再进行重新下载，大幅削减时间。

1.4)根据步骤1.2)得到的物种分类信息，从步骤1.3)得到的目标代谢途径的所有核酸和氨基酸序列中去掉真核生物的基因序列，并根据最终产生的基因编号(gene ID)，通过TBtools工具中的序列提取命令(Amazing fasta extractor)功能，特异性获得细菌和古菌的序列，作为最终使用的目标代谢途径的特异性数据库DB.fasta。

上述步骤2)中，建立目标代谢途径的特异性数据库的映象(mapping)文件的方法，包括以下步骤：

2.1)从建立的目标代谢途径的特异性数据库DB.fasta中，获得数据库的索引文件DB.fasta.fai。其中，获取索引文件的命令为samtools faidx DB.fasta；获得的数据库索引文件DB.fasta.fai中，第一列是特异性数据库中的基因名，第二列为该基因对应的氨基酸序列长度。

2.2)通过共有的基因名排序，将步骤2.1)中得到的索引文件DB.fasta.fai以及步骤1.3)得到的序列信息文件ko_pathway_information.txt进行合并，形成目标代谢途径的特异性数据库的mapping文件DB.txt。

上述步骤3)中，基于获得的目标代谢途径的特异型数据库DB.fasta，对宏基因组测序获得的clean reads进行数据库快速比对，得到各样品的比对结果的方法，包括以下步骤：

3.1)运行diamond makedb--in DB.fasta-d DB_nr命令，基于获得的目标代谢途径的特异型数据库DB.fasta，构建适用于diamond的数据库。其中，构建方法是本领域人员公知技术，在此不再赘述。

3.2)基于构建的适用于diamond的数据库，对宏基因组测序获得的clean reads序列文件进行快速比对，得到各样品的比对结果。

具体方法为：运行diamond blastx-d DB_nr–q命令输入序列文件.gz-o，得到输出结果.txt--evalue 1e-5--query-cover 75--id 90-k 1。其中，输入序列文件可以是压缩文件，通过合理设置比对条件，包括临界值，覆盖度，一致性等相关参数的设定，获得数据库的比对结果，也即序列的注释信息。

其中，本发明中，将基于蛋白质组数据库的比对条件设置为E-value≤10^-5，Identity≥90％，Coverage≥25AA，认定clean reads为目标基因。

上述步骤4)中，根据得到的各样品的比对结果以及特异型数据库DB.fasta的映像文件，进行排序、统计和整合的方法，包括以下步骤：

4.1)根据步骤2)中获得的特异型数据库的mapping文件DB.txt第一列中的基因名，进行排序，然后逐一统计计算各样品中每个基因比对到的reads个数，然后提取获得单个样品的数据库注释信息。

对注释到目标序列的reads进行统计的命令为：

for i in`cut-f 1<DB.txt>|sort-u`

do echo"echo'"$i"'>>raw.txt"

echo"grep-c'"$i"'*_*.txt>>raw.txt"

done|sh

4.2)将多个样品比对结果，按照样品名进行分列排序，运行：

cat raw.txt|awk-F':”{if(NF＝＝1){print}else{print$2}}'|xargs-n<样品数+1>|sed's//\t/g'-|sed"1i$(head-<样品数+1>raw.txt|awk-F':''{if(NF＝＝2){print$1}else{print'\t'}}'|xargs|sed's//\t/g'-)"-|awk'{if(NR＝＝1){print"ID\t",$0}else{print$0}}'>DB.xls

上述步骤5)中，对各样品的比对结果进行均一化处理，根据均一化处理结果进行不同样品间的定量分析的方法，包括以下步骤：

5.1)对不同样品的注释结果进行均一化处理，以消除不同样品测序深度不同带来的影响。

由于不同样品的测序深度可能有所不同，因此需要通过内参基因如16s RNA进行均一化，以便数据结果在不同样品、不同测序深度下进行比较。通过以下公式进行数据均一化：

式中，N_{targetgene-likesequence}为目标同源基因的数量；L_{referencesequence}为数据库中参考基因长度；L_reads宏基因组测序获得reads长度；N_{16S sequence}为宏基因组中16s rRNA的reads数量；L_{16S sequence}为16s rRNA比对数据库中的平均长度。以上各数值均可从步骤4)中得到的比对结果中得到。

具体的，通过ARGs-OAP程序获得样本所需均一化的结果，运行

./argoap_pipeline_stageone_version2-i<测序结果文件目录>-o<目标保存目录>-m meta-data.txt-s-n 8-f fa

结果产生的meta_data_online中有每个样品均一化所需要的数据，包括按照reads数/16s/cell number均一化所需要的数据。

5.2)按照目标代谢途径基因的对应的酶编号，对均一化后的注释信息做进一步统计加和，用于酶水平上基因的定量分析。

下面通过具体实施例对本发明方法做进一步介绍。

实施例1：

宏基因组测序中甲烷代谢(Methane mechanism)功能基因的数据挖掘。

宏基因组测序结果：12个双末端测序，测序深度5G；

目标：研究不同氨氮抑制条件对甲烷代谢的影响。

1.构建Methane mechanism的特异性数据库

1)运行perl kegg_pathway_extract.pl--ko_ID_file ko_ID.txt#ko_ID.txt文件(map00680)#

2)根据物种分类信息，去掉真核生物的基因序列，并根据最终产生的gene ID，通过TBtools工具Amazing fasta extractor功能，特异性获得细菌和古菌的序列，作为最终使用的特影响数据库Methane_mechanism.fasta。

2.建立Methane mechanism特异性数据库的mapping文件

1)从Methane mechanism特异性数据库中，获取Methane mechanism特异性数据库的索引文件，运行命令为samtools faidx Methane_mechanism.fasta

2)通过共有的基因名排序，将索引文件Methane_mechanism.fasta.fai和序列信息文件ko_pathway_information.txt进行合并,形成Methane mechanism特异性数据库的mapping文件Methane_mechanism.txt。

3.Methane mechanism数据库比对

diamond makedb--in Methane_mechanism.fasta-d Methane_mechanism_nr

diamond blastx-d Methane_mechanism_nr-q输入序列文件.gz-o输出结果.txt--evalue 1e-5--query-cover 75--id 90-k 1

4.比对结果整合

for i in`cut-f 1<Methane_mechanism.txt>|sort-u`；do echo"echo'"$i"'>>raw.txt"；echo"grep-c'"$i"'*_*.txt>>raw.txt"；done|sh

cat raw.txt|awk-F':”{if(NF＝＝1){print}else{print$2}}'|xargs-n 13|sed's//\t/g'-|sed"1i$(head-13raw.txt|awk-F':”{if(NF＝＝2){print$1}else{print'\t'}}'|xargs|sed's//\t/g'-)"-|awk'{if(NR＝＝1){print"ID\t",$0}else{print$0}}'>Methane_mechanism.xls

5.数据的均一化处理

./argoap_pipeline_stageone_version2-i<测序结果文件目录>-o<输出问价保存目录>-m meta-data.txt-s-f fa

最后按照公式进行均一化。

如图1a和图1b所示，为目标样本中甲烷代谢途径的结果。图1a中，加粗框线显示的为目标代谢途径中检测到的基因酶编号，图1b为酶水平上基因的丰度变化。

实施例2：

宏基因组测序中氮代谢(Nitrogen mechanism)功能基因的数据挖掘。

宏基因组测序结果：12个双末端测序，测序深度5G；

目标：研究不同氨氮抑制条件对厌氧消化中氮代谢的影响。

1.构建Nitrogen mechanism数据库

1)perl kegg_pathway_extract.pl--ko_ID_file ko_ID.txt#ko_ID.txt文件(map00910)#

2)根据物种分类信息，去掉真核生物的基因序列，并根据最终产生的gene ID，通过TBtools工具Amazing fasta extractor功能，特异性获得细菌和古菌的序列，作为最终使用的数据库Nitrogen_mechanism.fasta。

2.建立Nitrogen mechanism数据库mapping文件

1)samtools faidx Nitrogen_mechanism.fasta

2)合并Nitrogen_mechanism.fasta.fai和ko_pathway_information.txt,形成mapping文件Nitrogen_mechanism.txt。

3.Nitrogen mechanism数据库比对

diamond makedb--in Nitrogen_mechanism.fasta-d Nitrogen_mechanism_nr

diamond blastx-d Nitrogen_mechanism_nr-q输入序列文件.gz-o输出结果.txt--evalue 1e-5--query-cover 75--id 90-k 1

4.比对结果整合

for i in`cut-f 1<Nitrogen_mechanism.txt>|sort-u`；do echo"echo'"$i"'>>raw.txt"；echo"grep-c'"$i"'*_*.txt>>raw.txt"；done|sh

cat raw.txt|awk-F':”{if(NF＝＝1){print}else{print$2}}'|xargs-n 13|sed's//\t/g'-|sed"1i$(head-13raw.txt|awk-F':”{if(NF＝＝2){print$1}else{print'\t'}}'|xargs|sed's//\t/g'-)"-|awk'{if(NR＝＝1){print"ID\t",$0}else{print$0}}'>Nitrogen_mechanism.xls

5.数据的均一化处理

最后按照公式进行均一化。

如图2a和图2b所示，为目标样本中氮代谢途径的结果。图2a中，加粗框线显示的为目标代谢途径中检测到的基因酶编号，图2b为对应酶水平上基因的丰度变化。

实施例3：

宏基因组测序中硫代谢(Sulfur mechanism)功能基因的数据挖掘。

宏基因组测序结果：12个双末端测序，测序深度5G；

目标：研究不同氨氮抑制条件对厌氧消化中硫代谢的影响。

1.构建Sulfur mechanism数据库

1)perl kegg_pathway_extract.pl--ko_ID_file ko_ID.txt#ko_ID.txt文件(map00920)#

2)根据物种分类信息，去掉真核生物的基因序列，并根据最终产生的gene ID，通过TBtools工具Amazing fasta extractor功能，特异性获得细菌和古菌的序列，作为最终使用的数据库Sulfur_mechanism.fasta。

2.建立Sulfur mechanism数据库mapping文件

1)samtools faidx Sulfur_mechanism.fasta

2)合并Sulfur_mechanism.fasta.fai和ko_pathway_information.txt,形成mapping文件Sulfur_mechanism.txt。

3.Sulfur mechanism数据库比对

diamond makedb--in Sulfur_mechanism.fasta-d Sulfur_mechanism_nr

diamond blastx-d Sulfur_mechanism_nr-q输入序列文件.gz-o输出结果.txt--evalue 1e-5--query-cover 75--id 90-k 1

4.比对结果整合

for i in`cut-f 1<Sulfur_mechanism.txt>|sort-u`；do echo"echo'"$i"'>>raw.txt"；echo"grep-c'"$i"'*_*.txt>>raw.txt"；done|sh

cat raw.txt|awk-F':”{if(NF＝＝1){print}else{print$2}}'|xargs-n 13|sed's//\t/g'-|sed"1i$(head-13raw.txt|awk-F':”{if(NF＝＝2){print$1}else{print'\t'}}'|xargs|sed's//\t/g'-)"-|awk'{if(NR＝＝1){print"ID\t",$0}else{print$0}}'>Sulfur_mechanism.xls

5.数据的均一化处理

最后按照公式进行均一化。

如图3a和图3b所示，为目标样本中硫代谢途径的结果。图3a中，加粗框线显示的为目标代谢途径中检测到的基因酶编号，图3b为对应酶水平上基因的丰度变化。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims

1.一种宏基因组数据挖掘方法，其特征在于包括以下步骤：

2)建立特异性数据库DB.fasta的映像文件；

2.如权利要求1所述的一种宏基因组数据挖掘方法，其特征在于：所述步骤1)中，从KEGG数据库中获取目标代谢途径的全部基因信息，并建立特异性数据库DB.fasta的方法，包括以下步骤：

1.2)获取KEGG数据库上的物种分类信息；

3.如权利要求2所述的一种宏基因组数据挖掘方法，其特征在于：所述步骤2)中，建立特异性数据库DB.fasta的映像文件的方法，包括以下步骤：

4.如权利要求1～3任一项所述的一种宏基因组数据挖掘方法，其特征在于：所述步骤3)中，基于获得的目标代谢途径的特异型数据库DB.fasta，对宏基因组测序获得的cleanreads进行数据库比对，得到各样品的比对结果的方法，包括以下步骤：

3.2)基于构建的适用于diamond软件的数据库，对宏基因组测序获得的clean reads序列文件进行比对，得到比对结果。

5.如权利要求3所述的一种宏基因组数据挖掘方法，其特征在于：所述步骤4)中，根据得到的各样品的比对结果以及特异型数据库DB.fasta的映像文件，进行排序、统计和整合的方法，包括以下步骤：

4.2)将多个样品的比对结果按样品名进行分列排序；

6.如权利要求1～3任一项所述的一种宏基因组数据挖掘方法，其特征在于：所述步骤5)中，对各样品的比对结果进行均一化处理，根据均一化处理结果进行不同样品间的定量分析的方法，包括以下步骤：

5.1)对不同样品的比对结果进行均一化处理；

7.如权利要求6所述的一种宏基因组数据挖掘方法，其特征在于：所述步骤5.1)中，对不同样品的比对结果进行均一化处理时，所采用的均一化公式为：

式中，N_{targetgene-likesequence}为目标同源基因的数量；L_{referencesequence}为数据库中参考基因长度；L_reads宏基因组测序获得reads长度；N_16Ssequence为宏基因组中16s rRNA的reads数量；L_16Ssequence为16s rRNA比对数据库中的平均长度。