肠道菌群：宏基因组测序分析流程解读（上）
上回给大家讲述了16S测序分析流程，应各位小伙伴们的要求，本期的宏基因组测序分析来啦~
- 概述
- 宏基因组测序实验流程
- 宏基因组分析流程
- 文章分析思路
- 常见问题
- 相关名词解释
- 参考文献
16S rDNA扩增子测序（16S rDNA amplicons sequencing）和鸟枪法宏基因组测序（whole metageome shotgun sequencing）都是研究微生物组的重要方法，那么问题来了：这两者到底有什么区别呢？什么情况下做16S测序？什么情况下做宏基因组测序？
概述
先要来说说什么是宏基因组测序：
宏基因组 ( Metagenome) (也称元基因组) 。是由 Handelsman 等 1998 年提出的新名词, 其定义为“the genomes of the total microbiota found in nature” , 即生境中全部微小生物遗传物质的总和。它包含了可培养的和未可培养的微生物的基因, 目前主要指环境样品中的细菌和真菌的基因组总和。
宏基因组学 (或元基因组学, metagenomics) 就是一种以环境样品中的微生物群体基因组为研究对象, 以功能基因筛选和/或测序分析为研究手段, 以微生物多样性、 种群结构、 进化关系、 功能活性、 相互协作关系及与环境之间的关系为研究目的的新的微生物研究方法。
宏基因组测序（Metagenomics Sequencing）是对环境样品中全部微生物的总DNA进行高通量测序，主要研究微生物种群结构、基因功能、微生物之间的相互协作关系以及微生物与环境之间的关系。宏基因组测序研究摆脱了微生物分离纯培养的限制，扩展了微生物资源的利用空间，为环境微生物群落的研究提供了有效工具。
高通量测序平台
16S扩增子测序和宏基因组测序的主要区别如下：
测序原理不同
16S rDNA基因存在于所有细菌的基因组中，具有高度的保守性。该序列包含9个高变区和10个保守区（如下图），通过对某一段高变区序列（V4区或V3-V4区）进行PCR扩增后进行测序，得到300-500bp左右的序列。
宏基因组测序则是将微生物基因组DNA随机打断成500bp的小片段，然后在片段两端加入通用引物进行PCR扩增测序，再通过组装的方式，将小片段拼接成较长的序列。
研究目的不同
16S测序主要研究群落的物种组成、物种间的进化关系以及群落的多样性。
宏基因组测序在16S测序分析的基础上还可以进行基因和功能层面的深入研究，宏基因组测序可以回答这样的问题"who is there?"和 "what are they doing?"。
物种鉴定分辨率不同
16S测序得到的序列很多注释不到种水平，而宏基因组测序则能鉴定微生物到种水平甚至菌株水平。
对于16S测序而言，任何一个高变区或几个高变区，尽管具有很高的特异性，但是某些物种（尤其是分类水平较低的种水平）在这些高变区可能非常相近，能够区分它们的特异性片段可能不在扩增区域内。
宏基因组测序通过对微生物基因组随机打断，并通过组装将小片段拼接成较长的序列。因此，在物种鉴定过程中，宏基因组测序具有较高的优势。
Tips：通常情况下，建议同时结合宏基因组测序和16S测序两种技术手段，可以更高效、更准确地研究微生物群落组成结构、多样性以及功能情况。
如果样本污染宿主DNA比较严重，例如肠道粘膜样本，直接宏基因组测序会产生大量的宿主污染，为了降低实验成本，可以使用16S测序。
如果想快速鉴定未知病原感染，直接通过metagenome测序可以鉴定是细菌、真菌或者是病毒感染。
宏基因组测序实验流程
从环境（如土壤、海洋、淡水、肠道等）中采集实验样本，将原始采样样本或已提取的 DNA 样本低温运输（0℃以下），对样品进行样品检测。检测合格的 DNA 样品，进行文库构建以及文库检测，检测合格的文库将采用 Illumina 高通量测序平台进行测序，测序得到的下机数据(Raw Data)将用于后期信息分析。
为保证测序数据的准确性、可靠性，对样品检测、建库、测序每一个生产步骤都严格把控，从根本上确保高质量数据的产出，具体的实验流程图如下：
DNA样品检测
对 DNA 样品的检测主要包括 3 种方法：
(1) 琼脂糖凝胶电泳（AGE）分析 DNA 的纯度和完整性；
(2) Nanodrop 检测 DNA 的纯度（OD260/280 比值）；
(3) Qubit 对 DNA 浓度进行精确定量；
文库构建及质检
检测合格的 DNA 样品用超声波破碎仪随机打断成长度约为350bp的片段，经末端修复、3’端加A、加测序接头、纯化、片段选择、PCR 扩增等步骤完成整个文库制备。
文库构建完成后，先用电泳及Nanodrop进行初步定量，对浓度>=15ng/ul的文库进行Qubit定量，用毛细管电泳对文库的插入片段大小进行检测，插入片段大小符合预期后，使用qPCR方法对文库的有效浓度进行准确定量（文库有效浓度＞3nM），以保证文库上机质量。
上机测序
建库质检合格后，把不同文库按照有效浓度及目标下机数据量的需求，混合后进行Illumina测序。
宏基因组分析流程
下面先放一张宏基因组分析流程图，供小伙伴们快速了解一下。
测序数据预处理
•常规数据质控工具
ØFASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/commandline.html)
ØCutadapt (https://github.com/marcelm/cutadapt)
ØTrimmomatic (http://www.usadellab.org/cms/?page=trimmomatic)
ØSickle (https://github.com/najoshi/sickle)
•去除宿主污染
ØSoapAligner (http://soap.genomics.org.cn/soapaligner.html)
A member of the SOAP (Short Oligonucleotide Analysis Package). It is an updated version of SOAP software for short oligonucleotide alignment.
ØBowtie (http://sourceforge.net/projects/bowtie-bio/files/)
An ultrafast, memory-efficient short read aligner, which aligns short DNA sequences (reads) to the human genome at a rate of over 25 million 35-bp reads per hour.
采用 Illumina 测序平台测序获得的原始数据(Raw Data)存在一定比例低质量数据，里面含有带接头的、重复的，以及测序质量很低的reads，这些 reads 会影响组装和后续分析，为了保证后续分析的结果准确可靠，需要对原始的测序数据进行预处理，获取用于后续分析的有效数据(Clean Data)。
可以使用过滤软件 Trimmomatic，可以从任意一段切除低质量的碱基，同时还支持滑窗过滤，根据情况设定滑窗的大小和阈值，当滑窗内的碱基质量与设定的阈值进行比较，如果数值低于阈值则切除整个滑窗的碱基。高通量测序一般会包括接头序列以及引物片段，可以使用 Trimmomatic来去除这些序列。
例如具体处理步骤如下：
1) 去除所含低质量碱基（质量值≤38）超过一定比例（默认设为 40bp）的 reads；
2) 去除 N 碱基达到一定比例的 reads（默认设为10bp）；
3) 如果样品存在宿主污染，需与宿主数据库进行比对，过滤掉可能来源于宿主（可以采用 SoapAligner 软件，参数设置: identity ≥ 90%, -l 30, -v 7, -M 4,-m 200,-x 400）的 reads；
数据处理信息统计表
注：RawData 表示下机原始数据；CleanData 表示过滤得到的有效数据；CleanQ20，CleanQ30 表示 CleanData 中测序错误率小于0.01(质量值大于 20)和 0.001(质量值大于 30)的碱基数目的百分比；Clean_GC(%) 表示 CleanData 中碱基的 GC 含量；Effective(%) 表示有效数据( CleanData )与原始数据( RawData )的百分比。
物种注释metaphlan2
主页：http://segatalab.cibio.unitn.it/tools/metaphlan2/
从Clean read出发，使用metaphlan2软件分析，获得不同分类层级的物种丰度表。
MetaPhlAn2是分析微生物群落(细菌、古菌、真核生物和病毒)组成的工具，它在宏基因组研究中非常有用，只需一条命,即可获得微生物的物种丰度信息。同时配合自带的脚本可进一步统计和可视化。
MetaPhlAn2整理了超过17000个参考基因组，包括13500个细菌和古菌，3500个病毒和110种真核生物，汇编整理了>1百万类群特异的标记基因，可以实现：
精确的分类群分配
准确估计物种的相对丰度
种水平精度
株鉴定与追踪
超快的分析速度
Duy Tin Truong, Eric Franzosa, Timothy L Tickle, Matthias Scholz, George Weingart, Edoardo Pasolli, Adrian Tett, Curtis Huttenhower, and Nicola SegataMetaPhlAn2 for enhanced metagenomic taxonomic profilingNature Methods 12, 902–903 (2015)
看看软件主要能干什么
分析人类微生物组(HMP)数据，在样本间进行层级聚类
使用实例
定量物种谱
一条命令就可以从原始数据得到物种的相对丰度！！
最常用的是使用双端压缩fastq文件，参数--nproc调用多个线程，并输出比较结果
metaphlan2.py --bowtie2out metagenome.bowtie2.bz2 --nproc 20 --input_type fastq <(zcat metagenome_1.fq.gz metagenome_2.fq.gz) > profiled_metagenome1.txt
例如人类肠道的数据，一般几个小时内就出结果了。
输出结果为各层级物种相对丰度值，可以直接作为lefse的输入文件。
#SampleID Metaphlan2_Analysisk__Bacteria 100.0k__Bacteria|p__Bacteroidetes 58.71768k__Bacteria|p__Firmicutes 27.73866k__Bacteria|p__Verrucomicrobia 9.51729k__Bacteria|p__Proteobacteria 3.20978k__Bacteria|p__Actinobacteria 0.81659k__Bacteria|p__Bacteroidetes|c__Bacteroidia 58.71768k__Bacteria|p__Firmicutes|c__Bacilli 17.31386k__Bacteria|p__Firmicutes|c__Clostridia 10.26799k__Bacteria|p__Verrucomicrobia|c__Verrucomicrobiae 9.51729k__Bacteria|p__Proteobacteria|c__Deltaproteobacteria 2.43682k__Bacteria|p__Actinobacteria|c__Actinobacteria 0.81659k__Bacteria|p__Proteobacteria|c__Betaproteobacteria 0.62445k__Bacteria|p__Firmicutes|c__Erysipelotrichia 0.15682k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria 0.14851k__Bacteria|p__Bacteroidetes|c__Bacteroidia|o__Bacteroidales 58.71768k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales 17.31386
批量处理多个样品
假如我们有1-20一共20个样本，调用for循环并转后台并行处理
for i in `seq 1 20`;dometaphlan2.py --nproc 5 --input_type fastq <(zcat Sample_${i}.1.fq.gz Sample_${i}.2.fq.gz) > metaphlan2_${i}.txt &done
样品物种谱合并
mergemetaphlantables.py命令可以将每个样品结果表合并，程序位于程序的utils目录中，请自行添加环境变量或使用绝对路径。合并时支持输入文件多个文件空格分隔，或使用通配符(如下)。可结合sed删除样本名中共有部分，精简样品名方便可视化简洁美观
merge_metaphlan_tables.py metaphlan2*.txt | sed 's/metaphlan2_//g' > merged.txt
获得了矩阵表，下面可以进行各种统计分析与可视化啦！
种相对丰度概况
从不同分类层级的相对丰度表出发，选取出在各样品中的最大相对丰度排名前 15 的物种，绘制出各样品对应的物种注释结果在不同分类层级上的相对丰度柱形图。
基于GraPhlAn的分类学组成信息可视化
功能注释humann2
Species-level functional profiling of metagenomes and metatranscriptomes
Nature Methods, Article, 2018-10-30
DOI: http://dx.doi.org/10.1038/s41592-018-0176-y
第一作者： Eric A. Franzosa, Lauren J. McIver
通讯作者：Curtis Huttenhower
主要单位：哈佛医学院统计系；哈佛和麻省理工博德(Broad)研究所
humman2是一套分析流程，它包括调用metaphlan流程来分析物种组成，和自身分析功能基因和代谢通路组成。HUMAnN2 Workflow
humann2 --threads 20 --input test.fastq --output humann2_out/humann2_renorm_table --input abundance.tsv --output abundance_relab.tsv --units relabhumann2_join_tables --input humann2_out --output humann2_pathabundance.tsv --file_name pathabundance_relabmerge_metaphlan_tables.py *tsv > metaphlan2_merged.tsv
Metagenome 组装和注释
Metagenome 组装
两种常见策略：
1 基于序列overlap关系进行拼接
2 基于de Bruijn图进行组装
由于现阶段的主流测序方法是二代短片段测序，序列短而且数目庞大，如果利用overlap关系直接进行组装，这要求每对reads之间都进行一次序列比较，这会很耗费时间，而且结果并不可靠。为迎合二代测序的特点，一种基于k-mer的de Bruijn组装策略则成为更有效的解决方法。
能用于宏基因组组装的软件有很多，如metavelvet，SOAPdenovo，megahit，ibda-ud，metaspades等。
1、SOAPdenovo
SOAPdenovo is a novel short-read assembly method that can build a de novo draft assembly for the human-sized genomes， which specially designed to assemble Illumina GA short reads.
2、MEGAHIT
MEGAHIT is a single node assembler for large and complex metagenomics NGS reads, such as soil. Compare to SOAPdenovo, it generates longer contigs and consumes less memory.
3、IDBA-UD
IDBA-UD is a iterative De Bruijn Graph De Novo Assembler for Short Reads Sequencing data with Highly Uneven Sequencing Depth.
4、SPAdes和metaSPAdes
An assembly toolkit containing various assembly pipelines, which works with Illumina or IonTorrent reads and is capable of providing hybrid assemblies using PacBio, Oxford Nanopore and Sanger reads.
如下是三种主流软件分析，运行所消耗时间、内存比较:
目前MEGAHIT在现有组装软件中，资源消耗基本上是最低的，因此很适合宏基因组中的复杂环境样品。
还有上面介绍过的软件SPAdes，无论单菌、宏基因组还是宏病毒组都表现不错 metaSPAdes: a new versatile metagenomics assembler）中显示即使在复杂环境（土壤），组装效果也大大优于megahit、IDBA-UD等，但是没有megahit资源消耗低。
以下使用 SOAPdenovo[1]进行Metagenome 组装为例。
从各样品质控后的Clean reads出发，组装主要分为3步：单样本组装，多样本组装结果合并和丰度过滤。
1) 经过预处理后得到Clean Data，使用SOAPdenovo[1]组装软件进行组装分析(Assembly Analysis)；
2) 对于单个样品，首先选取一个K-mer(默认选取55)进行组装,得到该样品的组装结果；组装参数：-d1,-M3,-R,-u,-F
3) 将组装得到的Scaffolds从N连接处打断，得到不含N的序列片段，称为Scaftigs(i.e.,continuous sequences within scaffolds)；
4) 将各样品质控后的CleanData采用SoapAligner软件比对至各样品组装后的Scaftigs上，获取未被利用上的PE reads；比对参数：-u,-2,-m200
5) 将各样品未被利用上的reads放在一起，进行混合组装，考虑到计算消耗和时间消耗，只选取一个kmer进行组装(默认-K55)，其他组装参数与单样品组装参数相同；
6) 将混合组装的Scaffolds从N连接处打断，得到不含N的Scaftigs序列；
7) 对于单样品和混合组装生成的Scaftigs，过滤掉500bp以下的片段，并进行统计分析和后续基因预测；
各样品组装结果基本信息统计
说明：SampleID为样品名称；Total Length表示组装得 Scaftigs 的总长； Scaftigs 表示组装得到的Scaftigs 总条数； N50，N90 表示将 Scaftigs 按照长度进行排序，然后由长到短加和，当加和值达到 Scaftigs 总长的 50%，90%时的 Scaftigs 的长度值；Max len 表示组装得到的最长 Scaftigs 的长度值；Min len表示组装得到的最长 Scaftigs 的长度值。
样品的scaftigs长度分布统计如下：
各样品scaftigs数目统计：
组装结果质量评估工具 MetaQUAST
基因预测及丰度分析
基因预测及丰度分析基本步骤
1) 从各样品及混合组装的 Scaftigs（>=500bp）出发，采用MetaGeneMark进行ORF (OpenReading Frame) 预测，并从预测结果出发，过滤掉长度小于 100nt的信息；预测参数：采用默认参数进行
2) 对各样品及混合组装的 ORF 预测结果，采用CD-HIT软件进行去冗余，以获得非冗余的初始 genecatalogue（此处，操作上，将非冗余的连续基因编码的核酸序列称之为 genes），默认以 identity 95%,coverage 90% 进行聚类，并选取最长的序列为代表性序列；采用参数：-c 0.95, -G 0, -aS 0.9, -g 1, -d 0
3) 采用SoapAligner，将各样品的 Clean Data 比对至初始 gene catalogue，计算得到基因在各样品中比对上的reads 数目；比对参数：-m 200, -x 400, identity ≥ 95%
4）过滤掉在各个样品中支持 reads 数目<=2 的基因，获得最终用于后续分析的 genecatalogue（Unigenes）
5）从比对上的 reads 数目及基因长度出发，计算得到各基因在各样品中的丰度信息
6）基于 gene catalogue 中各基因在各样品中的丰度信息，进行基本信息统计，core-pan 基因分析，样品间相关性分析，及基因数目韦恩图分析。
Gene catalogue 基本信息
说明：ORFs NO. 表示 gene catalogue 中基因的数目；integrity:start 表示只含有起始密码子的基因数目及百分比；integrity:end 表示只含有终止密码子的基因数目及百分比；integrity:none 表示没有起始密码子也没有终止密码子的基因数目及百分比；integrity:all 表示完整基因（既有起始密码子也有终止密码子）数目的百分比；Total Len.(Mbp) 表示 gene catalogue 中基因的总长，单位是百万；Average Len. 表示 gene catalogue 中基因的平均长度；GC Percent 表示预测的 gene catalogue 中基因的整体 GC 含量值。
core-pan 基因分析
从基因在各样品中的丰度表出发，可以获得各样品的基因数目信息，通过随机抽取不同数目的样品，可以获得不同数目样品组合间的基因数目，由此我们构建和绘制了 Core 和 Pan 基因的稀释曲线，图片展示如下:
组间基因数目差异分析
为了考察组与组间的基因数目差异情况，绘制了组间基因数目差异箱图，展示结果如下：
说明：图中，横坐标为各个分组信息；纵坐标为基因数目。
基因数目韦恩图分析
为了考察指定样品（组）间的基因数目分布情况，分析不同样品（组）之间的基因共有、特有信息，绘制了韦恩图(Venn Graph)或花瓣图，展示结果如下：
说明：图中，每个圈代表一个样品；圈和圈重叠部分的数字代表样品之间共有的基因个数；没有重叠部分的数字代表样品的特有基因个数。
基于基因丰度的样品间相关性分析
样品间基因丰度相关性是检验实验可靠性和样本选择是否合理性的重要指标。相关系数越接近1，表明样品之间基因丰度模式的相似度越高。
说明：图中不同颜色代表相关系数的高低，相关系数与颜色间的关系见有图图例说明。
基于物种丰度的降维分析
1) PCA分析
主成分分析PCA（Principal component analysis）是一种研究数据相似性或差异性的可视化方法，通过一系列的特征值和特征向量进行排序后，选择主要的前几位特征值，采取降维的思想，PCA 可以找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样品点之间的相互位置关系，只是改变了坐标系统。PCA 可以观察个体或群体间的差异。下图每一个点代表一个样本，相同颜色的点来自同一个分组，两点之间距离越近表明两者的群落构成差异越小。
说明：种水平PCA分析，横坐标表示第一主成分，百分比则表示第一主成分对样品差异的贡献值；纵坐标表示第二主成分，百分比表示第二主成分对样品差异的贡献值；图中的每个点表示一个样品，同一个组的样品使用同一种颜色表示。
2) NMDS分析
NMDS（Non-metric Multidimensional scaling）非度量多维尺度分析是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。适用于无法获得研究对象间精确的相似性或相异性数据，仅能得到他们之间等级关系数据的情形。其基本特征是将对象间的相似性或相异性数据看成点间距离的单调函数，在保持原始数据次序关系的基础上，用新的相同次序的数据列替换原始数据进行度量型多维尺度分析。其特点是根据样品中包含的物种信息，以点的形式反映在多维空间上，而对不同样品间的差异程度，则是通过点与点间的距离体现的，最终获得样品的空间定位点图。
组间差异物种的LEfSe分析
为了筛选组间具有显著差异的物种Biomarker，首先通过秩和检验的方法检测不同分组间的差异物种并通过LDA（线性判别分析）实现降维并评估差异物种的影响大小，即得到LDA score；组间差异物种的LEfSe分析结果包括三部分，分别是LDA值分布柱状图，进化分支图（系统发育分布）和组间具有统计学差异的Biomarker在不同组中丰度比较图。差异物种的LDA值分布图和进化分支图如下：
常用功能数据库注释
功能注释主要是基于功能同源的序列往往具有序列相似性的原理，将去冗余后的基因蛋白序列，与不同的蛋白功能数据库进行序列比对， 然后用比对到的序列的功能作为目标序列的功能。
从 gene catalogue 出发，进行代谢通路（KEGG）[6,7]，同源基因簇（eggNOG）[8]，碳水化合物酶（CAZy）的功能注释和丰度分析。基于物种丰度表和功能丰度表，可以进行丰度聚类分析，PCA和NMDS 降维分析，Anosim分析，样品聚类分析；当有分组信息时，可以进行Metastat和LEfSe多元统计分析以及代谢通路比较分析，挖掘样品之间的物种组成和功能组成差异。
目前常用的功能数据库主要有：
KEGG，全称Kyoto Encyclopedia of Genes and Genomes，是一个关于基因功能注释方面的综合性数据库，包括基因的功能，分类，代谢通路(KEGG Pathway数据库， 是KEGG最核心的数据库)等诸多方面的信息。
KEGG 数据库于 1995 年由 Kanehisa Laboratories 推出 0.1 版，目前发展为一个综合性数据库，其中最核心的为KEGG PATHWAY 和 KEGG ORTHOLOGY 数据库。在 KEGG ORTHOLOGY 数据库中，将行使相同功能的基因聚在一起，称为Ortholog Groups (KO entries)，每个 KO 包含多个基因信息，并在一至多个 pathway 中发挥作用。
KEGG Pathway数据库将生物代谢通路划分为6大类(A级分类)，分别为：新陈代谢（Metabolism）、遗传信息处理（Genetic Information Processing）、 环境信息处理（Environmental Information Processing）、细胞过程（Cellular Processes）、生物体系统（Organismal Systems）、 人类疾病（Human Diseases），其中每大类又被系统分类为B、C、D3个级别。 其中B级分类目前包括有 43 种子功能；C级分类即为代谢通路图；D级分类为每个代谢通路图的具体注释信息。
eggNOG 数据库是利用 Smith-Waterman 比对算法对构建的基因直系同源簇 (Orthologous Groups) 进行功能注释， eggNOG（v4.5）目前涵盖了2031个物种，构建了包含25个大类，约19万个直系同源簇。
CAZy 数据库是研究碳水化合物酶的专业级数据库，主要涵盖 6 大功能类：糖苷水解酶（Glycoside Hydrolases ，GHs）， 糖基转移酶（Glycosyl Transferases，GTs），多糖裂合酶（Polysaccharide Lyases，PLs），碳水化合物酯酶（Carbohydrate Esterases，CEs）， 辅助氧化还原酶(Auxiliary Activities , AAs)和碳水化合物结合模块（Carbohydrate-Binding Modules， CBMs）。 其中每一个大类有可以分类很多小的家族，比如CE1，CE2等等，注释结果中的CE0表示没有小家族分类的结果。
功能注释基本步骤
1）使用DIAMOND软件将 Unigenes 与各功能数据库进行比对（blastp，evalue ≤ 1e-5）；
2）比对结果过滤：对于每一条序列的 比对结果，选取 score 最高的比对结果（one HSP > 60 bits）进行后续分析；
3）从比对结果出发，统计不同功能层级的相对丰度（各功能层级的相对丰度等于注释为该功能层级的基因的相对丰度之和），其中，KEGG 数据库划分为 5 个层级，eggNOG 数据库划分为 3 个层级，CAZy 数据库划分为 3 个层级，各数据库的详细划分层级如下所示：
| 数据库名称 | 划分层级 | 该层级的描述 |
|---|---|---|
| KEGG | Level 1 | KEGG 代谢通路第一层级 6 大代谢通路； |
| KEGG | Level 2 | KEGG 代谢通路第二层级 43 种子 pathway； |
| KEGG | Level 3 | KEGG pathway id（例：ko00010）； |
| KEGG | ko | KEGG ortholog group (例：K00010)； |
| eggNOG | Level 1 | 24大功能类 |
| eggNOG | Level 2 | ortholog group description |
| eggNOG | og | ortholog group ID |
| CAZy | Level 1 | 6大功能类 |
| CAZy | Level 2 | CAZy family |
| CAZy | Level 3 | EC number |
4）从功能注释结果及基因丰度表出发，获得各个样品在各个分类层级上的基因数目表，对于某个功能在某个样品中的基因数目，等于在注释为该功能的基因中，丰度不为 0 的基因数目；
5）从各个分类层级上的丰度表出发，进行注释基因数目统计，相对丰度概况展示，丰度聚类热图展示，PCA和NMDS降维分析，基于功能丰度的Anosim组间（内）差异分析，代谢通路比较分析，组间功能差异的Metastat和LEfSe分析。
注释基因数目统计
从 Unigenes 注释结果出发，绘制各个数据库的注释基因数目统计图，展示结果如下图所示：
各数据库注释基因数目统计图
说明：从上至下依次为CAZy，eggNOG 的 Unigenes 注释数目统计图。横坐标轴是各数据库中 level1 各功能类的代码，代码的解释见对应的图例说明。
功能相对丰度概况
从各数据库 level1 的相对丰度表出发，绘制出各个数据库中，各样品对应的在 level1 层级上的丰度统计图。
KEGG:
Eggnog:
功能注释在 level1 上的相对丰度柱形图
说明：从上至下依次为 Kegg,eggNOG 的结果展示。纵轴表示注释到某功能类的相对比例；横轴表示样品名称；各颜色区块对应的功能类别见右侧图例。
功能相对丰度聚类分析
根据所有样品在各个数据库中的功能注释及丰度信息，选取丰度排名前 35 的功能及它们在每个样品中的丰度信息绘制热图，并从功能差异层面进行聚类。
Kegg:
Eggnog:
功能丰度聚类热图
说明：从上至下依次为 KEGG, eggNOG 的结果展示。横向为样品信息；纵向为功能注释信息；
基于功能丰度的降维分析
基于不同数据库在各个分类层级的功能丰度进行 PCA 和 NMDS 降维 分析，如果样品的功能组成越相似，则它们在降维图中的距离越接近。基于KEGG的KO功能丰度进行PCA和NMDS 分析的结果展示如下：
抗生素抗性基因注释
抗生素的滥用导致人体和环境中微生物群落发生不可逆的变化，对人体健康和生态环境造成风险，因此抗性基因的相关研究受到了研究者的广泛关注。
目前，抗生素抗性基因的研究主要利用ARDB（Antibiotic Resistance Genes Database）数据库和CARD（Comprehensive Antibiotic Resistance Database）数据库来注释抗生素抗性基因。通过该数据库的注释，可以找到抗生素抗性基因（Antibiotic Resistance Genes ，ARGs）及其抗性类型（Antibiotic Resistance Class）以及这些基因所耐受的抗生素种类（ Antibiotic）等信息。
抗性基因丰度概况
依据耐受的抗生素种类丰度结果，提取丰度前15抗生素种类绘制柱状图如下：
不同抗生素抗性基因在各样品中的丰度柱形图
注：横轴为样品，纵轴为抗生素抗性基因
抗性基因类型丰度聚类分析
根据抗性基因类型的丰度表，绘制聚类热图：
注：横轴为样品，纵轴为抗性基因类型，不同颜色达标图例中对应的丰度范围。
组间抗性基因数目差异分析
注：横轴为分组，纵轴为抗性基因类型数目
物种间相关性网络分析
机器学习算法
分类器是通过特征工程中的特征选择的方法挑选出若干不同级别的bio-markers，并根据bio-markers来构建的一种分类模型，它采用分类评价指标来评估模型效果的好坏，目的是用于识别疾病与健康人群。
分类器是来自机器学习中的分类算法，包括Logistictic回归、Support Vector Machine (SVM)、Random Forest、人工神经网络、朴素贝叶斯等。除了分类算法，特征选择对分类性能也会有直接的影响，特征选择就是在我们获得的给定特征集中选择出与分类相关的特征子集的过程。
分类评价指标是用来评价分类器分类性能好坏的一种指标，常见的评价指标有正确率、错误率、灵敏度、特异度、精度、Matthews correlation coefficient (MCC)、Receiver operating characteristic (ROC)。其中MCC和ROC不受分组样本数量影响；ROC是根据一系列不同的二分类方式，以真阳性率(灵敏度)为纵坐标，假阳性率(1-特异度)为横坐标绘制的曲线，直观描绘了分类器在TP和FP间的trade-off。
AUC (area under curve)是ROC曲线下方面积之和，它以数值的形式来评估分类器的好坏，AUC的取值范围在(0.5, 1.0)范围内，AUC的值越大，分类的性能就越好。若分类评价指标最终显示分类器性能不佳，则需重新构建分类器直至分类性能达到稳定为止。
分类器的引入对于构建用于临床诊断的基于肠道微生物的检测方法具有重要意义。
文章分析思路
文章分析思路：整体概述——物种组成——功能组成——关联分析——因果关系验证
文章主要结果
样本背景、数据评估、物种、功能描述
物种组成整体差异和分组具体差异
功能组成整体差异和分组具体差异
环境因子与微生物组的关系
菌群或单菌移植到无菌小鼠体内，验证因果关系
常见问题
16S扩增子和宏基因组分析结果为什么存在差别？
1、两者的实验方法存在较大差异：
16S是先扩增后测序，而且不同物种DNA的扩增有偏好；在宏基因组测序中，得到的相对物种丰度的差异与测序深度、DNA提取以及测序的方法都密切相关；
2、两者采用的物种注释方法及数据库都存在着一定差别：
16S采用的是将16S rDNA与Greengenes/Silva数据库进行比对注释，只能注释到细菌或古菌；而宏基因组则是将预测得到的基因与NR数据库或特定的marker基因数据库比对从而进行注释，宏基因组注释得到的物种信息更为全面，不仅包括细菌，还包括真菌、古菌以及病毒等。
很多文章在物种丰度水平更多使用16s的结果，而功能注释则使用metagenome的结果。
宏基因组分析策略
宏基因组研究主要有两种分析策略：
第一种有参考基因组，直接将reads比对到参考基因组上进行研究。
第二种为无参考基因组，需要通过序列组装、基因预测、再进行物种功能注释，是一种基于de novo组装的研究。
我们平时所做的宏基因组研究大多数是把两种研究方法结合起来。
为什么基因组的组装不完整？
宏基因组组装的效果主要跟以下几个因素有关：样本的测序数据量，物种多样且丰度分布不均匀等，这些因素都会造成宏基因组组装比细菌等单物种的组装更加困难，这也是目前宏基因组研究中有待突破的重点。
随着测序读长不断增加，测序质量的不断提高，将三代测序数据应用到宏基因组中将是未来研究的一大方向。
为什么不把所有样本合并在一起进行组装？
不同样本中物种丰度差异很大，如果把所有样本都混合在一起，对服务器的要求很高（例如需要大内存服务器），将会大大增加数据的复杂度，组装效果可能会更差。
宏基因组测序的测序技术选择，二代测序还是三代？
二代测序的短读长限制了contig的长度，基因之间的物理位置很难准确确定。
随着三代测序技术成熟，基于三代PacBio的SMRT单分子测序和Nanopore长读长特性可以跨越大的重复区域，使得宏基因组拼接组装效果有了很大程度的提高，能组装出更长的contig，增加了完整基因的数量，可以对复杂的重复区域进行分析。
宏基因组测序深度的选择？
1）由于受到测序深度及测序成本的影响，在现在的宏基因组文章中，测序数据量一般选择大于5G，可以测出样品中绝大多数的微生物，但是对于一些低丰度的物种，因为测序深度的原因，确实很有可能会组装不出来；
2）在宏基因组分析中，也一般多关注的是较高丰度物种的组成情况，如果要对低丰度物种进行特殊分析，一般需要加大测序数据量，或者在前期提取过程中经过一些特殊的处理，尽可能的富集出多的低丰度物种，再进行测序分析。
关于宏基因组的测序量，这个要根据样品的复杂度来看：
如果diversity比较低，例如像人肠道之类的， 每个样品测个5-10Gb的数据就可以了，对于复杂样品，例如土壤，致使测序几十到几百G，也可能也会深度不足。
如果不知道样品的diversity，一般建议先对样品做一个 16S rRNA测序的survey，来看一下你的样品里面大致有多少种微生物。
抗性基因和毒力基因分析所用数据库是什么？
随着人们对抗性基因相关研究的广泛关注，我们宏基因组的标准分析中推出了抗性基因的相关分析。并且，由于自2009年ARDB数据库再无更新，因此我们目前所用的抗性基因数据库为CARD数据库。
VFDB（Virulence Factors Database）数据库用来注释毒力基因。
环境因子与微生物组的关系
RDA 或者CCA定义是基于对应分析发展而来的一种排序方法，将对应分析与多元回归分析相结合。此分析是主要用来反映菌群与环境因子之间关系。RDA是基于线性模型，CCA 是基于单峰模型。分析可以检测环境因子、样品、菌群三者之间的关系或者两两之间的关系。
注：图中数字表示样品名，不同颜色或形状表示不同环境或条件下的样本组；箭头表示环境因子；图中蓝色上三角表示不同属的细菌；物种与环境因子之间的夹角代表物种与环境因子间的正、负相关关系（锐角：正相关；钝角：负相关；直角：无相关性）；由不同的样本向各环境因子做垂线，投影点越相近说明样本间该环境因子属性值越相似，即环境因子对样品的影响程度相当。
如何进行肠型分析？
2011年肠型（Enterotypes）概念第一次被提出，肠型被分为3种Bacteroides (enterotype 1), Prevotella (enterotype 2) and Ruminococcus (enterotype 3)，文中指出：
1）在不同的肠道样本中，均有明显的类别存在；
2）3类肠型中的优势菌，在不同样本中是不同的。
影响肠型结果的各种因素：
1）不同聚类方法、距离计算方式最终得到的肠型会有不同
2）不同OTU的层级得到的肠型不同
3）不同测序区域得到的肠型不同
4）OTU的计算策略得到的肠型不同
5）测序策略得到的肠型不同
肠型分类器的在线使用
肠型在线分类器主页: http://enterotypes.org/基于MetaHIT数据集，建立了一个基于属水平的在线分类器(基于HMP1和中国二型糖尿病 )
宏基因组测序能否鉴定到菌株水平？
CAG (co-abundance gene group)
CAG方法是2014年在《自然∙生物技术》杂志上报道出来的人体肠道宏基因组数据分析方法。这种方法的思想和MGS的思想有一定相同之处，都利用了这样一种思想或理论：在同一个株或种中的基因，这些基因的丰度在群体中具有丰度一致性。为了解决使用MGS方法时仅关注基因标记的局限性，canopy算法被用来鉴定共丰富基因簇(CAG)，CAG分析基于基因丰度鉴定物种，并且每个CAG可以被看作是部分微生物或完整微生物。
当基因丰度表包含足够数量的样本时，可以应用canopy算法来识别CAG。canopy算法使用Pearson相关系数和Spearman相关系数作为阈值，在对只有1个基因的canopies进行第一轮聚类和过滤之后，根据簇的平均丰度，使用Pearson相关系数作为阈值再次被应用到canopy算法中，第二轮簇可以包含重叠基因。
因此，对于出现在多于1个簇中的基因，基因及其相关簇之间的距离能够被确定，对于每个重叠基因，选择最近的簇。最后，选择含有超过700个基因的簇最可能是含有细菌基因组部分的CAG。
基于CAG分析的结果，可以在人肠道微生物群中鉴定出可能显著影响宿主的未分类或难以理解的微生物(也称为生物暗物质)。在过去的研究中，生物暗物质一般被忽略，因此，许多疾病相关的肠道微生物可能未被检测到。CAG中的重叠群和基因可以产生无法从当前的细菌基因组数据库获得的大量新信息。
两款鉴定到菌株的软件
PanPhlAn：PanPhlAn - Strain-level resolution in Metagenomics（一种基于菌株特异性泛基因组学的分析方法）
ConStrains: ConStrains identifies microbial strains in metagenomic datasets
相关名词解释
De novo测序 从头测序，无需任何参考序列，直接对一个物种进行测序，然后进行拼接、组装成该物种的基因组序列图谱。
建库 在基因组随机打断的片段上机前，为保证在测序时有足够的数据强度支持，所进行的基于PCR反应的片段扩增，区别于通常说的基于Fosmid/BAC等载体的建立文库。
Paired-End测序 双末端测序，对插入片段两端进行测序，产生具有Paired-End关系的reads。
测序策略sequences strategyPE100：即Paired-End (100,100)，采用双末端测序法对插入片段两端进行读长为100 bp的高通量测序。
读长 测序仪所能获取的实际长度，即reads的长度，例如：90bp、125bp、300bp。
Raw data（reads）即下机数据，原始的reads。
Clean data（reads）即下机数据经过去接头污染、过滤低质量reads、去duplication（针对大片段数据）等之后，实际用于组装及分析的数据。
Contig序列 由来源于同一基因组，具有overlapping关系的reads拼接而成的片段。
Scaffold序列 又称super Contig，通过reads的Paired-End关系将Contig连接，并用N补充其中的内洞而形成。
N50/N90 将组装所得片段（Scaffold/Contig）按照从长到短排序并累加求和，累加值达到基因组总长度一半时的片段长度即是该组装结果的N50值，通常用来衡量组装情况；N90与之类似，即累加长度为基因组总长90%时，该片段的长度。
基因集（Gene Catalog）由所有样本中检测到的基因构成的集合。
非冗余基因集（Non-redundant Gene Catalog）将所有样本的基因以一定的阈值进行聚类，每个类别取最长的基因作为代表序列。
物种分类学注释：将基因集序列与参考数据库比对，并根据分类学谱系系统得到该序列的物种分类地位，从而分析得到整个基因集的物种信息；
COG功能注释：将基因集序列与COG参考数据库比对，获得基因集的同源聚类蛋白群的相应功能注释信息；
KEGG功能注释：将基因集序列与KEGG的基因组数据库比对，获得基因集在代谢、酶、反应和功能模块方面的注释信息。
CAZy数据库注释：将基因集序列与特殊数据库比对，如 CAZy数据库（Carbohydrate-Active enZYmes Database，碳水化合物活性酶）
CARD数据库（Comprehensive Antibiotic Resistance Database，抗生素抗性基因）， 获得基因集特殊的功能信息
参考文献
[1] Luo et al.: SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. GigaScience 2012 1:18.
[2] Li J, Jia H, Cai X, et al. An integrated catalog of reference genes in the human gut microbiome[J]. Nature biotechnology, 2014, 32(8): 834-841.
[3] Oh J, Byrd A L, Deming C, et al. Biogeography and individuality shape function in the human skin metagenome[J]. Nature, 2014, 514(7520): 59-64.
[4] Zhu, Wenhan, Alexandre Lomsadze, and Mark Borodovsky. "Ab initio gene identification in metagenomic sequences." Nucleic acids research 38.12 (2010): e132-e132
[5] Duy Tin Truong, Eric A Franzosa,et al. MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nature Methods 12, 902-903 (2015)
[6 Kanehisa M, Goto S, Hattori M, Aoki-Kinoshita KF, Itoh M, Kawashima S, et al. (2006). From genomics to chemical genomics: new developments in KEGG. Nucleic Acids Res 34(Database issue): D354–7.
[7 Kanehisa, M., Goto, S., Sato, Y., Kawashima, M., Furumichi, M., and Tanabe, M.; Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 42, D199–D205 (2014).
[8] Powell S, Forslund K, Szklarczyk D, et al. eggNOG v4. 0: nested orthology inference across 3686 organisms[J]. Nucleic acids research, 2013: gkt1253.
[9] Liu B, Pop M. ARDB—Antibiotic Resistance Genes Database[J]. Nucleic Acids Research, 2009, 37(Database issue):D443-7.
[10] Li W, Godzik A: Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 2006, 22(13):1658-1659.
[11] Fu L, Niu B, Zhu Z, Wu S, Li W: CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 2012, 28(23):3150-3152.
相关阅读：