一种基于宏基因组学分析精准识别水体中未知微生物群落的方法
阅读说明：本技术 一种基于宏基因组学分析精准识别水体中未知微生物群落的方法 (Method for accurately identifying unknown microbial community in water body based on metagenomic analysis ) 是由 陈倩 钟浩辉 刘树枫 倪晋仁 于 2021-01-25 设计创作，主要内容包括：本发明公开了一种基于宏基因组学分析精准识别水体中未知微生物群落的方法,所述方法包括步骤(1)自水样提取宏基因组DNA；(2)DNA测序；(3)根据目标群落选择性构建参考数据库；(4)对测序数据进行组装获得组装数据；(5)对组装数据进行分箱；(6)对分箱数据进行质量测试,标记MAGs的质量,并计算测序深度；(7)根据构建的参考数据库对组装数据进行注释；(8)对MAGs做进化关系分析,进而做出宏基因组群落结构分析。本发明不依赖更新速度较慢的大型软件配套数据库,适合处理未知物种较多的样本,且可检测极低丰度物种,检测全面、快速。(The invention discloses a method for accurately identifying an unknown microbial community in a water body based on metagenomic analysis, which comprises the following steps of (1) extracting metagenomic DNA from a water sample; (2) DNA sequencing; (3) selectively constructing a reference database according to the target community; (4) assembling the sequencing data to obtain assembly data; (5) binning the assembly data; (6) performing quality test on the boxed data, marking the quality of the MAGs, and calculating the sequencing depth; (7) annotating the assembly data according to the constructed reference database; (8) and carrying out evolutionary relationship analysis on the MAGs so as to further carry out metagenome community structure analysis. The method does not depend on a large software matched database with a slow updating speed, is suitable for processing samples with more unknown species, can detect extremely low abundance species, and is comprehensive and rapid in detection.)
技术领域
本发明涉及生物信息学方法在环境监测领域的应用，具体涉及一种基于宏基因组学分析精准识别水体中未知微生物群落的方法。
背景技术
近年来随着测序技术的发展，对环境微生物的研究方式日新月异。早期环境微生物的研究基本基于纯培养或形态生理生化等表型数据，但因为培养方式只能获取环境微生物中极小的部分，仅观察表型特征也难以研究其系统发生，这个时期对环境微生物的群落研究往往难以得到完整图景。
高通量测序技术的发展为微生物群落研究提供了新的方向。单基因水平上，传统的方法包括基于引物的标记基因如16S扩增子研究等，因为受引物限制较大，一部分微生物基因序列和常用引物无法匹配，导致群落结构研究中此类微生物被忽视，不能全面反映群落结构；同时所得片段较小，种系区分率不佳，难以精细分类；且物种缺乏完整基因组无法预测其代谢。而近年来发展起的对环境中所有DNA片段测序的宏基因组技术能够在基因组水平上研究群落结构，同时能将群落结构和每个分类所拥有的基因功能结合研究，比起基于引物测序标记基因的研究方式能提供更多相关信息。
中国专利申请202010628901.0公开了一种优化的宏基因组binning分析微生物群落的方法，包括对测序数据进行过滤，得到高质量测序数据，然后根据样本的来源及测序数据量的大小，选择不同的组装策略得到contigs，接着进行基因数据分析。在宏基因组层面，提供了更贴合样本特征、测序数据量的高效优质组装算法，并包含丰富全面的信息分析内容，个性新颖的可视化。有利于更方便高效地筛选到有价值的目标bin。
然而，常用的流程化分析软件和方法在宏基因组研究上依赖已有数据库，运行时需要检测全部数据速度较慢；而数据库维护上往往具有多年的滞后性，现行稳定的软件版本和数据库通常不包括最近两三年鉴定分类的微生物；且不同软件不同数据库结构往往不一致或非常见结构，自己手动为特定软件补充数据库较为复杂。
因此，确立一种基于宏基因组，能简便地跟随最新论文研究分类成果精准分类，操作可自定义程度高的群落结构监测与识别方法，可以有效地提升我们对环境微生物的群落结构了解。
发明内容
为解决现有技术中存在的不足，本发明的目的在于提供一种基于目标群落选择性构建参考数据库，并在此基础上优化数据处理流程，进而实现对水体中微生物群落的精准识别。该方法不依赖更新速度较慢的大型软件配套数据库，适合处理未知物种较多的样本，且可检测极低丰度物种，检测全面、快速。
本发明采用如下的技术方案：
一种基于宏基因组学分析精准识别水体中未知微生物群落的方法，所述方法包括步骤：
(1)自水体取样，提取待测样本宏基因组DNA；
(2)对DNA测序，并进行测序数据质量控制；
(3)根据目标群落选择性构建参考数据库；
(4)对质控后测序数据进行组装获得组装数据，并获取其测序深度数据；
(5)对组装数据进行分箱：将数据划分在多个单独文件中，每个划分结果视为同一种生物的基因组，称为“由宏基因组组装的基因组”(Metagenome-Assembled Genome，简称MAG)；
(6)对分箱数据进行质量测试，标记MAGs的质量，并计算测序深度；
(7)根据构建的参考数据库对组装数据进行注释，获取单拷贝标记基因的注释与测序深度，计算出宏基因组测序的平均深度(测序细胞数)和MAGs相对丰度；
(8)对MAGs做进化关系分析，进而做出宏基因组群落结构分析；
所述步骤(4)，采用MEGAHIT或metaSPAdes软件对质控后数据进行组
装获取重叠群contigs，采用BBMap软件将质控后读段以最佳相似度比对mapping到组装好的重叠群上，获取每条重叠群的测序深度，
所述测序采用Illumina测序，读段长度为150bp，比对时相似度阈值设为97％以上。
所述步骤(3)，通过论文或NCBI taxonomybrowser最新基因组数据选取参考基因组，手动构建参考数据库。
优选地，根据分类需要，从论文资料或者NCBI的分类浏览器(https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi)选取所需参考分类级别，获得txid号，从NCBI的Assembly(https://www.ncbi.nlm.nih.gov/assembly/)根据该txid打包下载所需分类下所有基因组genomic sequence格式的fasta文件。
可选地，一部分较新或质量较低的基因组无法从assembly网页打包下载，可根据提示的GCA或GCF号，选用ftp服务器(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/)手动或爬取下载。
可选地，为进一步加快分类速度，可进一步手动挑选精简下载的参考基因组，同一分类级别下分类尽量分散，例如门级别分类的话挑选的参考基因组应位于不同纲下。
所述步骤(1)，采样方式为水样原位过滤，所用滤膜为0.22微米；使用成熟DNA提取试剂盒提取DNA。
所述步骤(2)，DNA测序采用Illumina HiSeq双端150bp×2，下机读段reads
使用fastp软件进行质量控制获取质控读段。
优选地，下机数据去除嵌合，质量分数低于35的以及含有未知碱基(N)超过10％的读段，后续使用fastp软件(v0.21.0)可视化读段统计数据，根据双端测序重合部分自动纠错，修正重叠部分不一致碱基中质量分数较低的，去除长度低于140bp的读段，并切除平均质量分数最低的最后一位碱基。
所述步骤(5),使用分箱软件metaBAT2、CONCOCT或MaxBin2中一种，或者3者整合的metaWRAP流程，输入步骤(4)中的重叠群测序深度以及重叠群数据自身，对重叠群数据进行分箱，获得MAGs。
优选地，选择metaSPAdes进行组装(v3.14.1，--meta选项)获取重叠群。
可选地，30Gbp读段(约60GB fastq文件)运行metaSPAdes约需要200-500GB内存。若样本数据量较大，系统内存较小，选择MEGAHIT进行组装(v1.2.9)获取重叠群。
优选地，使用BBMap软件(v38.86)对组装好后的重叠群做比对获取sam文件；使用samtools软件(v1.10)将sam文件转换为二进制bam文件并排序。使用metabat2软件(v2.12.1)中自带的jgi_summarize_bam_contig_depths计算排序bam中各重叠群的测序深度。
优选地，使用metabat2软件(v2.12.1)，输入测序深度和重叠群数据进行分箱，得到MAGs。
可选地，若时间和计算资源充裕，使用CONCOCT(v1.1.0)或MaxBin2(v2.2.7)或者整合的metaWRAP(v1.3)流程，输入测序深度和重叠群数据进行分箱，得到MAGs。
所述步骤(6)，采用CheckM软件对分箱结果进行质量检测，标记高质量与低质量分箱bin，并记录由CheckM测算的粗略分类位置信息；根据步骤(4)中测序深度计算出MAGs的测序深度。
可选地，根据CheckM粗略分类结果挑选感兴趣的MAGs。
优选地，使用CheckM软件(v1.1.2)软件对分箱产生的MAGs进行质量检测，获取其完整度和污染度数据，记录每个MAG的质量。同时该软件也会给出每个MAGs的粗略分类信息。
可选地，根据分类和MAG质量信息，剔除完整度<50％，污染度>10％且MAG测序深度低的MAG。
所述步骤(6)，根据公式：
计算每个MAG的测序深度。
所述步骤(7)，使用prodigal软件预测组装好的重叠群和参考基因组中的氨基
酸序列，并使用HMMER软件建立参考用的多种单拷贝标记基因的hmm文件，再使用hmmsearch预测单拷贝标记基因，然后根据参考基因计算宏基因组总测序深度(可看作宏基因组中测序的细胞数目)和MAGs的相对丰度。
优选地，使用prodigal软件(v2.6.3)预测重叠群中的基因，获取氨基酸序列文件，并根据MAG中含有的重叠群名称获取每个MAG的氨基酸序列文件。
优选地，选取10种核糖体蛋白作为单拷贝标记基因代表，经研究发现它们是原核生物基因组中分布最为广泛，单拷贝性佳，序列保守，且都属于蛋白质合成功能大类中，适合用于系统发生分析，如下：
rpS10 rpL4 rpL2 rpL3 rpL22 rpL14 rpS5 rpS2 rpL1 rpS9
可选地，其他常用单拷贝标记基因如上述10种以外的核糖体蛋白、fusA、recA、gyrB等也可以使用。
优选地，从Pfam数据库下载各个基因的hmm文件。
可选地，从其他数据库下载各个基因的氨基酸序列fasta文件，使用HMMER软件(v3.3.1)中的hmmbuild从fasta文件计算出hmm文件。
优选地，使用HMMER软件(v3.3.1)中的hmmsearch从宏基因组氨基酸序列文件中预测出所需的单拷贝标记基因。
所述步骤(7)，依据公式：
计算每个宏基因组样本测序深度以及每个MAG相对丰度。
所述步骤(8)，根据步骤(6)中粗略分类信息和参考基因组，再选用步骤(7)中得
到的多个单拷贝标记基因，用MAFFT对齐后，使用进化树软件如iqtree，建立MAGs与参考基因组的进化树，得到MAGs的进化位置，总结出分类信息。
优选地，用上述同样的方法获取参考基因组中的单拷贝标记基因氨基酸序列，将它们和MAG的同种基因复制到一个fasta文件中。若MAG中存在同一种标记基因多拷贝则可能是污染，多条序列均舍弃不用。使用MAFFT(v7.471)多序列比对，对齐各个基因的序列，再将每个基因组所包含的单拷贝标记基因氨基酸序列按同样的基因前后顺序串接(concatenating)，若有基因缺失则用空缺符号代替，删除掉缺失基因种数超过总数目一半的基因组，保证串接后每个基因组仅有一条串接序列，长度相等且序列位置比对合适。
优选地，使用iqtree软件(v2.1.2)，输入串接后的fasta文件，对参考基因组及MAGs做系统发生分析建立极大似然树(ML树)。根据分类位置可为每个基因组详细分类。
可选地，根据实际需要使用更复杂的替代模型增加可靠性。
可选地，若数据量大建树时间会极长，可以选用软件运行初期产生的生物邻接树(bioNJ树)做近似分类筛选高相似度MAGs。bioNJ树在门之间分类可能有较大误差，但是对和参考基因组相似度较高的MAGs分类较好。剔除分类较好的MAGs后再完整计算ML树得到其他MAGs的分类位置。
优选地，统计MAGs分类信息，相对丰度，得到样本的群落结构信息。
与现有技术相比，本发明的有益效果在于：
1.本发明可从微生物丰度较低的水体中进行富集采样，利用二代测序技术、宏基因组组装、分箱(Binning)等技术手段，基于多种单拷贝标记基因(绝大多数微生物中有且仅有单拷贝，且适合用于系统发生分析的基因)的分类，在基因组水平上识别出微生物的群落结构。
2.本发明基于目标群落选择性构建参考数据库，并在此基础上优化数据处理流程，进而实现对水体中微生物群落的精准识别。
3.本发明不依赖更新速度较慢的大型软件配套数据库，适合处理未知物种较多的样本，或处理大规模样本时仅挑选出自身希望研究的特定小类群，具有检测全面，快速，可检测极低丰度物种，可检测最新发现物种的特点，适用于大规模采样调查研究；同时获取的宏基因组信息有利于后续流程分析，适合于分析研究特殊环境水体。
4.本发明基于宏基因组检测群落结构，检测精度较高，比起传统常用方法：16SrRNA基因扩增子的群落检测法，其不依赖引物、能更全面地反映环境中的微生物；其不依赖培养，且比培养的方法检测更为灵敏。在检测群落结构方面更加全面，细致。
5.本发明得到的宏基因组分类信息可以在后续分析中与宏基因组功能分析结合，比传统通过群落结构推测功能的研究方法准确度更高。更适合精准研究特定功能与群落结构的关系。
6.本发明所采用的宏基因组流程经过优化，选用软件均为分析中较为稳定软件，在对MAGs分类进行鉴定的步骤时，和同类的宏基因组分析流程相比，因其选用的标记基因经过精选，分析准确率在未明显降低的情况下，速度比采用更多种类标记基因的研究方法更快(相较约30-50种标记基因的常见方法时间缩短约2/3-4/5，若选用标记基因更少则准确率有所下降)。且系统发生分析时不依赖软件的配套数据库，具有更加灵活性，适合最新发布的分类单元下宏基因组物种结构鉴定；且比对使用的基于结构域的HMMER比常用的基于相似度的DIAMOND或BLAST等比对方式更适合保守序列，比对出的基因更适合用于进化分析。综合来说该宏基因组分析流程比传统流程，在鉴定未知物种较多、采样站点较为特殊的宏基因组群落结构时更具优势。
附图说明
图1为本发明整个分析处理流程图；
图2为本发明利用包含参考基因组的数据库和所得中高质量的MAGs建立TACK超门ML进化树。
具体实施方式
下面结合附图对本申请作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本申请的保护范围。
本实施例是基于宏基因组学分析精准识别水体中未知微生物群落的方法，具体涉及在TACK古菌的门级别上对地下水宏基因组进行检测。TACK古菌是近年来提出的一个古菌超门，内含10个左右门，其中多个为近年来新鉴定的门，在很多数据库中缺乏对应分类级别。如图1所示，具体步骤如下：
(1)自水体取样，提取待测样本宏基因组DNA。
(2)对DNA测序，并进行测序数据质量控制。
(3)根据目标群落选择性构建参考数据库。
(4)对质控后测序数据进行组装获得组装数据，并获取其测序深度数据。
(5)对组装数据进行分箱：将数据划分在多个单独文件中，每个划分结果视为同一种生物的基因组，称为“由宏基因组组装的基因组”MAG。
(6)对分箱数据进行质量测试，标记MAGs的质量，并计算测序深度。
(7)根据构建的参考数据库对组装数据进行注释，获取单拷贝标记基因的注释与测序深度，计算出宏基因组测序的平均深度和MAGs相对丰度。
(8)对MAGs做进化关系分析，进而做出宏基因组群落结构分析。
步骤(1)，采样方式为水样原位过滤，所用滤膜为0.22微米；使用成熟DNA提取试剂盒提取DNA。
步骤(2)，DNA测序采用Illumina HiSeq双端150bp×2，下机读段reads使用fastp软件进行质量控制获取质控读段。
优选地，下机数据去除嵌合，质量分数低于35的以及含有未知碱基(N)超过10％的读段，后续使用fastp软件(v0.21.0)可视化读段统计数据，根据双端测序重合部分自动纠错，修正重叠部分不一致碱基中质量分数较低的，去除长度低于140bp的读段，并切除平均质量分数最低的最后一位碱基。
步骤(3)，通过论文或NCBI taxonomy browser最新基因组数据选取参考基因组，手动构建参考数据库。
优选地，根据分类需要，从论文资料或者NCBI的分类浏览器(https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi)选取所需参考分类级别，获得txid号，从NCBI的Assembly(https://www.ncbi.nlm.nih.gov/assembly/)根据该txid打包下载所需分类下所有基因组genomic sequence格式的fasta文件。
可选地，一部分较新或质量较低的基因组无法从assembly网页打包下载，可根据提示的GCA或GCF号，选用ftp服务器(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/)手动或爬取下载。
可选地，为进一步加快分类速度，可进一步手动挑选精简下载的参考基因组，同一分类级别下分类尽量分散，例如门级别分类的话挑选的参考基因组应位于不同纲下。
步骤(4)，采用MEGAHIT或metaSPAdes软件对质控后数据进行组装获取重叠群contigs，采用BBMap软件将质控后读段以最佳相似度比对mapping到组装好的重叠群上，获取每条重叠群的测序深度。
测序采用Illumina测序，读段长度为150bp，比对时相似度阈值设为97％以上。
步骤(5),使用分箱软件metaBAT2、CONCOCT或MaxBin2中一种，或者3者整合的metaWRAP流程，输入步骤(4)中的重叠群测序深度以及重叠群数据自身，对重叠群数据进行分箱，获得MAGs。
优选地，选择metaSPAdes进行组装(v3.14.1，--meta选项)获取重叠群。
可选地，30Gbp读段(约60GB fastq文件)运行metaSPAdes约需要200-500GB内存。若样本数据量较大，系统内存较小，选择MEGAHIT进行组装(v1.2.9)获取重叠群。
优选地，使用BBMap软件(v38.86)对组装好后的重叠群做比对获取sam文件；使用samtools软件(v1.10)将sam文件转换为二进制bam文件并排序。使用metabat2软件(v2.12.1)中自带的jgi_summarize_bam_contig_depths计算排序bam中各重叠群的测序深度。
优选地，使用metabat2软件(v2.12.1)，输入测序深度和重叠群数据进行分箱，得到MAGs。
可选地，若时间和计算资源充裕，使用CONCOCT(v1.1.0)或MaxBin2(v2.2.7)或者整合的metaWRAP(v1.3)流程，输入测序深度和重叠群数据进行分箱，得到MAGs。
步骤(6)，采用CheckM软件对分箱结果进行质量检测，标记高质量与低质量分箱bin，并记录由CheckM测算的粗略分类位置信息；根据步骤(4)中测序深度计算出MAGs的测序深度。
可选地，根据CheckM粗略分类结果挑选感兴趣的MAGs。
优选地，使用CheckM软件(v1.1.2)软件对分箱产生的MAGs进行质量检测，获取其完整度和污染度数据，记录每个MAG的质量。同时该软件也会给出每个MAGs的粗略分类信息。
可选地，根据分类和MAG质量信息，剔除完整度<50％，污染度>10％且MAG测序深度低的MAG。
步骤(6)，根据公式：
计算每个MAG的测序深度。
步骤(7)，使用prodigal软件预测组装好的重叠群和参考基因组中的氨基酸序列，并使用HMMER软件建立参考用的多种单拷贝标记基因的hmm文件，再使用hmmsearch预测单拷贝标记基因，然后根据参考基因计算宏基因组总测序深度和MAGs的相对丰度。
优选地，使用prodigal软件(v2.6.3)预测重叠群中的基因，获取氨基酸序列文件，并根据MAG中含有的重叠群名称获取每个MAG的氨基酸序列文件。
优选地，选取10种核糖体蛋白作为单拷贝标记基因代表，经研究发现它们是原核生物基因组中分布最为广泛，单拷贝性佳，序列保守，且都属于蛋白质合成功能大类中，适合用于系统发生分析，如下：
rpS10 rpL4 rpL2 rpL3 rpL22 rpL14 rpS5 rpS2 rpL1 rpS9
可选地，其他常用单拷贝标记基因如上述10种以外的核糖体蛋白、fusA、recA、gyrB等也可以使用。
优选地，从Pfam数据库下载各个基因的hmm文件。
可选地，从其他数据库下载各个基因的氨基酸序列fasta文件(尽可能全面包括古菌细菌)，使用HMMER软件(v3.3.1)中的hmmbuild从fasta文件计算出hmm文件。
优选地，使用HMMER软件(v3.3.1)中的hmmsearch从宏基因组氨基酸序列文件中预测出所需的单拷贝标记基因。
步骤(7)，依据公式：
计算每个宏基因组样本测序深度以及每个MAG相对丰度。
步骤(8)，根据步骤(6)中粗略分类信息和参考基因组，再选用步骤(7)中得到的多个单拷贝标记基因，用MAFFT对齐后，使用进化树软件如iqtree，建立MAGs与参考基因组的进化树，得到MAGs的进化位置，总结出分类信息。
优选地，用上述同样的方法获取参考基因组中的单拷贝标记基因氨基酸序列，将它们和MAG的同种基因复制到一个fasta文件中。若MAG中存在同一种标记基因多拷贝则可能是污染，多条序列均舍弃不用。使用MAFFT(v7.471)多序列比对，对齐各个基因的序列，再将每个基因组所包含的单拷贝标记基因氨基酸序列按同样的基因前后顺序串接(concatenating)，若有基因缺失则用空缺符号代替，删除掉缺失基因种数超过总数目一半的基因组，保证串接后每个基因组仅有一条串接序列，长度相等且序列位置比对合适。
优选地，使用iqtree软件(v2.1.2)，输入串接后的fasta文件，对参考基因组及MAGs做系统发生分析建立极大似然树(ML树)。根据分类位置可为每个基因组详细分类。
可选地，根据实际需要使用更复杂的替代模型增加可靠性。
可选地，若数据量大建树时间会极长，可以选用软件运行初期产生的生物邻接树(bioNJ树)做近似分类筛选高相似度MAGs。bioNJ树在门之间分类可能有较大误差，但是对和参考基因组相似度较高的MAGs分类较好。剔除分类较好的MAGs后再完整计算ML树得到其他MAGs的分类位置。
优选地，统计MAGs分类信息，相对丰度，得到样本的群落结构信息。
更具体的步骤如下：
1.1样品采集
1.2样品处理
1.3DNA提取
2.1宏基因组测序
使用Illumina HiSeq型号测序仪，片段长度为，双端测序150bp×2，下机数据去除嵌合，质量分数于35的序列。
2.2质量控制
使用fastp软件进一步控制质量，去除140bp以下读段，删除最末尾一个碱基因为其质量太低，根据双端测序重叠部分纠错(选取质量较高的一侧的碱基)，具体命令为(方括号内文件名根据需要替换)：
fastp-w 16-c-l 140-t 1-h[report.html]-j[jsonfile.json]-i[reads_R1.fastq.gz]-I[reads_R2.fastq.gz]-o[clean_reads_R1.fastq.gz]-O[clean_reads_R2.fastq.gz]
3.1构建数据库
从NBCI的taxonomy browser选取参考基因组，获取TACK古菌类群的txid:1783275。从右侧assembly进入网页，可打包下载所有所需基因组。一部分基因组会提示可用(valid)但不范围内(not in scope)，这些可以根据提示GCA或GCF号到ftp服务器下载，例如GCA_aaabbbccc.1的基因组可以进入ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/aaa/bbb/ccc/下的文件夹里下载后缀名为genomic.fna.gz的文件。
3.2挑选精简数据库
在taxonomybrowser中勾选assembly选项，该页面中从每个门中挑选出所属类群不同的20多个基因组，进一步根据所列组装信息、组装质量等精选出5-10个。
或根据已有论文，对特定门类研究较仔细的情况下可以挑选常见的进化树用参考基因组。
4.1宏基因组组装
因为宏基因组数据量较大，选用MEGAHIT作为组装软件，命令为：
megahit--k-min 19--k-max 139--k-step 10-m 200000000000-t 14-1[clean_reads_R1.fastq.gz]-2[clean_reads_R2.fastq.gz]-o[output_dir]--out-prefix[output_filename_prefix]
组装完成后可以在[output_dir]中找到[output_filename_prefix]开头的.contigs.fa文件即为组装结果。
4.2测序深度
使用BBMap做读段比对，输入组装结果：
BBMap in＝[clean_reads_R1.fastq.gz]in2＝[clean_reads_R2.fastq.gz]threads＝14outm＝[mapped_reads.sam]ref＝[prefix.contig.fa]nodiskminid＝0.95
接着用samtools做格式转换并排序方便后续处理：
samtools [email protected]–bS[mapped_reads.sam]>[mapped_reads.bam]
samtools [email protected][sorted_reads.sbam][mapped_reads.bam]
然后使用metabat2中自带的脚本统计每条重叠群的测序深度，该脚本自动设置相似度阈值为97％：
jgi_summarize_bam_contig_depths--outputDepth[depth_file.tab][sorted_reads.sbam]
5.1分箱
选用metabat2进行分箱：
metabat2-i[prefix.contig.fa]-o[bins_output_dir]/[prefix]-a[depth_file.tab]-t14
运行完成后会在[bins_output_dir]文件夹下面生成前缀名为prefix，后缀名为.fa的MAG文件。
6.1质量测试
使用checkM软件中的lineage_wf流程：
checkm lineage_wf-t 14-x fa[bins_output_dir][checkm_output_dir]>[checkm_result.tab]
完成后可以在[checkm_result.tab]文件里看到各个MAG文件的质量进行筛选。本实施例中选择筛去质量低于50％，污染大于10％的MAGs，
6.2计算测序深度
获取每个MAGs中重叠群的名称，在linux下进入[bins_output_dir]，输入：
ls*.fa|xargs grep'>'>[all_bin_contigs.tab]
用表格软件等打开[all_bin_contigs.tab][depth_file.tab]统计每条重叠群信息，根据公式：
计算出每个MAG的测序深度。
6.3挑选TACK古菌
checkM流程做出了粗略分类，但因为其数据库近几年未更新，TACK古菌在checkM中仅能注释到古菌一个界上。从上述筛选后的MAGs中挑选出所有古菌界的MAGs后作为潜在TACK古菌进入后续分析。
7.1标记基因注释
使用prodigal软件，对[.contig.fa]文件做预测：
prodigal-i[.contig.fa]-a[predicted_amino_acids.faa]
使用下载的10种核糖体蛋白的[ribo.hmm]文件分别预测标记基因，比对阈值选择1e-5：
hmmsearch--tblout[hmm_out.tab]-E 1e-5[ribo.hmm][predicted_amino_acids.faa]
7.2计算测序深度
根据此前得到的每个重叠群测序深度，[depth_file.tab]文件以及10种[hmm_out.tab]中的信息，依据公式：
计算每个宏基因组样本测序深度以及每个MAG相对丰度。
8.1获取参考基因组及MAG的标记基因
根据[hmm_out.tab]的信息，若有重叠群属于古菌MAG，则可以挑选出来，将属于同一个MAG的所有10种核糖体蛋白先挑选出来，剔除核糖体蛋白种类数不满5种的MAG，将剩下的核糖体蛋白氨基酸fasta序列做好各自所属MAG标记后，各自放在一个文件[small/large_ribo_no.x.faa]里，共10种。
采用7中步骤获取标记基因氨基酸序列，储存到对应的标记基因文件中。
8.2标记基因对齐串接
使用MAFFT软件多序列比对：
mafft--maxiterate 1000--thread 14--localpair[small/large_ribo_no.x.faa]>[small/large_ribo_no.x.alignment.faa]
用文本编辑器打开[small/large_ribo_no.x.alignment.faa]去除所有序列间换行，然后粘贴到表格软件中排序，将10种标记基因放到同一个表格中然后串接，空缺部分全用-符号按同样长度补全，保证得到的每个基因组所含的每条串接后氨基酸序列长度一致。将结果以fasta格式保存至[all_MAG_ribo.faa]。
8.3建立进化树
使用iqtree软件预测模型：
iqtree-s[all_MAG_ribo.faa]-m MF
根据得到的模型[model]，建树第一步：
iqtree-s[all_MAG_ribo.faa]-m[model]
短暂运行后能很快得到一棵bioNJ树，立刻停止运行，根据bioNJ树结果，删去TACK超门以外的MAGs，保存到[TACK_ribo.faa]。树为newick格式，可视化软件可选用itol、mega等等。继续建树：
iqtree-s[TACK_ribo.faa]-m[model]-nt 10
8.4信息分析解读
得到ML树后，使用进化树可视化软件itol、mega等观察树形态，根据树内每个MAG的位置确定其分类。如图2所示，结尾为binxxx的均是MAGs。本示例研究的地下水中TACK古菌分布于两个门：深古菌门(Bathyarchaeota)以及奇古菌门(Thaumarchaeota)中，奇古菌又分为两个主要类群，具有氨氧化功能的氨氧化古菌(AOA)以及不能氨氧化的奇古菌(non-AOA)。其它用作参考但并未发现有MAGs所属的门均予以折叠表示。
根据MAG的相对丰度数据，即可计算出每个样本内TACK古菌各自的丰度、多样性等。
本发明基于目标群落选择性构建参考数据库，并在此基础上优化数据处理流程，进而实现了对水体中微生物群落的精准识别。
以上通过具体实施例详细描述了本发明的一种具体实践方式，本领域技术人员应当理解，本发明并不仅限于上述实施例，例如研究范围可以扩大到整个细菌域，整个古菌域等，分类水平也可根据参考基因组选取密度和对该类微生物的分类详略个别细化至纲目科属种水平。在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。