肠道菌群：16S测序分析流程解读

-   概述

-   16S测序实验流程

-   16S分析流程

-   结果怎么看

-   常见问题

-   参考文献

概述

你知道吗？居住在我们肠道内的细菌数量，是人体细胞总数的10倍之多！我们每天排出的粪便中，50%以上的干重量是由这些细菌及其“尸体”构成的。因此有人打趣的说，从数量上来看，我们人类并不应该被称为人类，而应被称作细菌。如此庞大的细菌群体驻扎在肠道内，构成了一个极为复杂的群体，被称作肠道菌群。肠道菌群被认为与人类健康息息相关，据估计，每个人的肠道菌群可包括～500种细菌。

那么，如何鉴别肠道菌群的种类和丰度呢？实验室常用的研究方法主要有哪些呢？

• 常规的基于分离培养的方法，该方法较为成熟，依然被许多肠道菌群的研究者采用，但是肠道菌群大多不能培养且比较耗时耗力，无法全面的分析对于种类和数量巨大的肠道微生态系统；

• 质谱技术，基于对菌体蛋白质分析的质谱技术具有在复杂体系中同时区分不同蛋白组分的特点，为肠道菌群的快速、准确、高通量的鉴别提供了可靠方法；

• 对于不容易分离培养的微生物，基因检测更显优势，常用的有基于分子杂交技术的分析方法、基于 DNA 指纹图谱的分析方法、基于 DNA 测序的检测方法等。

近年来，随着高通量测序技术的发展，16S rRNA基因测序技术在细菌的鉴定与分类研究中发挥着越来越重要的作用。

16S rRNA基因普遍存在于细菌细胞，在细菌基因组中位于核糖体小亚基（约1540 bp），该区域兼顾保守性和高变性，含有10 个保守区域(Conserved Regions)和 9 个高变区域(variable Regions)，保守区可用于设计引物进行目的片段的扩增，而通过对高变区的分析可以辨别细菌种类。因此，16S rRNA基因被认为是最适于细菌系统发育学研究和物种分类鉴定。目前用于16S rRNA基因深度测序的区域主要有V4区，V3-V4区、和V4-V5区等。

16S测序实验流程

将检测合格的环境微生物DNA样本对其指定区域进行PCR扩增、文库制备、文库质检、定量，使用设定的TAG序列进行样本区分。采用Illumina Hiseq 2500高通量测序平台对检测合格的文库进行测序。实验流程如下图所示：

16S分析流程

16S分析流程主要包括：Hiseq/Miseq测序获得的Paired-end (PE) reads拼接成一条序列，对目标序列进行质控过滤，过滤后的序列与参考数据库作比对，去除嵌合体序列得到最终得优化序列。基于优化序列进行OTU聚类分析和物种分类注释，基于OTU聚类结果进行多样性指数分析，基于分类学信息进行物种结构分析和物种差异分析。

数据质控与优化

数据质量评估

根据Fastq文件对测序样品进行数据质量评估（可以参考生信入门：Fasta与Fastq格式文件详解详细了解fastq文件）。单个样品的碱基质量分布如下图所示：

注： 横坐标为reads的碱基位置，纵坐标为单碱基平均Phred值。前250bp 为双端测序序列的Read1的碱基质量值分布情况，后250bp为read2 的碱基质量值分布情况。

也可以使用一文搞定细菌基因组De Novo测序分析中提到的fastqc和fastp进行质控。

如果想确定拿到手里的序列fastq序列就是16S的一部分，而不是其他的神马鬼，可以在NCBI的blast上比一比（可参考生信入门：序列比对之blast在线和本地使用）。

序列拼接

根据PE reads之间的overlap采用Flash软件对数据进行拼接。

1. #注意参数-M, --max-overlap=NUM 的设置，根据测序策略和PCR产物的长短决定。

2. flash -o extendedFrags test_1.fas tq test_2.fastq -M 180 -t 8

数据过滤

由于测序过程中会引入错误或者不可靠碱基，严重影响后续分析结果准确性。因此，采用fastx-toolkit工具过滤数据，只保留高质量（Q值>=25）的碱基比例大于等于90%的reads。

1. fastq_quality_filter -h

2. usage: fastq_quality_filter [-h] [-v] [-q N] [-p N] [-z] [-i INFILE] [-o OUTFILE]

3. Part of FASTX Toolkit 0.0.13.2 by A. Gordon ([e-mail])

4. [-h] = This helpful help screen.

5. [-q N] = Minimum quality score to keep.

6. [-p N] = Minimum percent of bases that must have [-q] quality.

7. [-z] = Compress output with GZIP.

8. [-i INFILE] = FASTA/Q input file. default is STDIN.

9. [-o OUTFILE] = FASTA/Q output file. default is STDOUT.

10. [-v] = Verbose - report number of sequences.

11. If [-o] is specified, report will be printed to STDOUT.

12. If [-o] is not specified (and output goes to STDOUT),

13. report will be printed to STDERR.

1. fastq_quality_filter -i extendedFrags.fastq -p 90 -q 25 -Q33

嵌合体序列检测

嵌合体在遗传学上用以指不同遗传性状嵌合或混杂表现的个体。嵌合体序列的出现会导致测序结果中产生一些实际并不存在的核酸序列，影响结果的可靠性。

因此，可以用usearch 64-bit软件进行嵌合体序列的检测及过滤。

usearch 64-bit是收费的，主要优势是支持大内存处理海量数据，它是超快的序列分析软件，在序列比对、聚类、操作等多领域广泛应用。

1. usearch8.0 -uchime_ref input.fasta -db gold.fasta -strand plus -nonchimeras good.fasta

2. # -db 推荐的数据为"Gold" database（http://drive5.com/uchime/uchime_download.html），ITS的数据，推荐使用UNITE数据库。

OTU聚类及物种注释

OTU（Operational Taxonomic Units，操作分类单元）是在系统发生分析或群体遗传研究中的一个假定的分类单元，通过一定的距离度量方法计算两两不同序列之间的距离度量或相似性，继而设置特定的分类阈值，获得同一阈值下的距离矩阵，进行聚类操作，形成不同的分类单元。

OTU聚类的目的：

1. 每个OTU（97%）在级别上对应一个种/属；

2. 每个OTU挑选出一个代表序列参与后续分析，节约计算资源；

3. 减少PCR或测序过程中引入的错误（错误序列与其来源序列较为相似，会聚成一个OTU）。

根据97%的序列相似性水平，利用QIIME软件包中的Uclust 方法进行OTU聚类分析。然后基于Silva（例如：Release132 https://www.arb-silva.de/documentation/release-132）参考数据库，对每个样品的OTUs进行物种分类学（Taxonomy）注释。

OTU聚类分析

OTU聚类部分结果如下表所示：

物种结构分析

统计所有样本在门、纲、目、科、属各层次上的分类结果。基于丰度前十五的物种采用累积柱状图比较样本间的物种组成差异，并进行列表展示。全样本在门层次的群落结构分析结果如下图所示：

物种丰度聚类热图

基于属水平上每个样品的物种丰度信息，选取丰度前50的属，根据各样品的丰度信息对样品和物种进行聚类，绘制热图。便于观察样品对应物种含量得高低。结果如下图所示：

单样本多级物种组成图

利用KRONA软件对物种注释结果进行可视化展示， 圆圈图得各层依次代表物种得分类级别，扇形的大小代表注释物种的比例。

Alpha多样性分析

Alpha多样性是指一个特定区域或生态系统内的多样性。多样性指数是反映样本中微生物丰富度和均匀度的综合指标。

Alpha多样性分析主要包括Shannon多样性指数、Simpson指数、Chao1丰富度指数等多样性指数的计算.

稀释曲线分析

采用随机抽样的方法，从样本中随机抽取一定数量的个体，统计出这些个体所代表物种数目，并以个体数与物种数来构建曲线，即稀释性曲线（rarefaction curve）。它可以用来比较测序数量不同的样本物种的丰富度，也可以用来说明样本的取样大小是否合理。

稀释曲线是利用已测得16S rRNA序列中已知的各种OTU的相对比例，来计算抽取n个（n小于测得reads序列总数）reads时出现OTU数量的期望值，然后根据一组n值（一般为一组小于总序列数的等差数列）与其相对应的OTU数量的期望值（此处采用chao1算法估计），当曲线趋向平坦时，说明测序数据量合理，更多的数据量只会产生少量新的OTU，反之则表明继续测序还可能产生较多新的OTU。

Chao1：用chao1算法估计群落中含OTU数目的指数，在生态学中常用来估计物种总数，由Chao在1984年最早提出。计算公式如下：

其中chao1是估计的OTU数；是实际观察的OTU数；F1是只有一条序列的OTU数目（如“singletons”）；F2是只有两条序列的OTU数目（如“doubletons”）。绘制各样本稀释曲线如下图所示：

由图可知，当测序量（横坐标）增加再多，对OTU数（纵坐标）的增长贡献不大，说明测序量已经足够，增加测序量不能提供更多的新的OTU。

Shannon-Wiener曲线

Shannon-Wiener指数是由Shannon, C. E. and Weaver W在1948年提出。反映样品中微生物多样性的指数，利用各样品的测序量在不同测序深度时的微生物多样性指数构建曲线，以此反映各样本在不同测序数量时的微生物多样性。

绘制Shannon曲线如下图所示：

Specaccum物种累积曲线

物种累积曲线( species accumulation curves) 是用于描述随着样品数量的加大而物种增加的状况，是调查样品的物种组成和预测样品中物种丰度的有效工具，被广泛用于样品量是否充分的判断以及物种丰富度( species richness) 的估计。因此，通过物种累积曲线不仅可以判断样品量是否充分，在样品量充分的前提下，运用物种累积曲线还可以对物种丰富度进行预测（默认在样品量大于10个分析）。

Beta多样性分析

Beta Diversity是对不同样本/不同组间样本的微生物群落构成进行比较分析。

可以用OTUs的丰度信息进行样本间距离计算， 也可以用OTUs之间的系统发生关系进行计算。

UniFrac是用于比较生物群落的一种距离度量，它在计算过程中包含了遗传距离（phylogenetic distances）的信息, 根据构建的进化树枝的长度计量两个不同环境样品之间的差异。

UniFrac 分析分为weighted unifrac 和unweighted unifrac 两种度量方法，两者之间差异在于是否计入不同环境样品的序列相对丰度。

Unweighted UniFrac只考虑了物种有无的变化，因此结果中，距离为0表示两个微生物群落间OTU的种类一致。而Weighted UniFrac则同时考虑物种有无和物种丰度的变化，距离为0则表示群落间OTU的种类和数量都一致。

实际应用的合理性：在环境样本的检测中，由于影响因素复杂，群落间物种的组成差异更为剧烈，因此往往采用非加权方法进行分析。但如果要研究对照与实验处理组之间的关系，例如研究短期青霉素处理后，人肠道的菌落变化情况，由于处理后群落的组成一般不会发生大改变，但群落的丰度可能会发生大变化，因此更适合用加权方法去计算。

PCoA分析

主坐标分析（Principal Coordinate Analysis，PCoA）将多维数据进行降维，是一种研究数据相似性或差异性的可视化方法，基于距离矩阵来寻找主坐标，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，有效地找出数据中最“主要”的元素和结构，对样本之间的关系进行描述。首先随机取样计算各样本之间的unifrac距离，再根据距离矩阵绘制二维PCoA图，如果样本的距离越接近（即物种的丰度和构成越相似），则它们在PCoA图中的距离越近。结果可见下图 。

NMDS分析

NMDS（Non-metric Multidimensional scaling）非度量多维尺度分析是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。适用于无法获得研究对象间精确的相似性或相异性数据，仅能得到他们之间等级关系数据的情形。其基本特征是将对象间的相似性或相异性数据看成点间距离的单调函数，在保持原始数据次序关系的基础上，用新的相同次序的数据列替换原始数据进行度量型多维尺度分析。其特点是根据样品中包含的物种信息，以点的形式反映在多维空间上，而对不同样品间的差异程度，则是通过点与点间的距离体现的，最终获得样品的空间定位点图。

差异分析

LEfSe分析

LEfSe分析即LDA Effect Size分析，是一种用于发现高维生物标识和揭示基因组特征的分析，使用non-parametric factorial Kruskal-Wallis (KW) sum-rank test（非参数因子克鲁斯卡尔—沃利斯和秩验检）检测具有显著丰度差异特征，并找到与丰度有显著性差异的类群。最后，LEfSe采用线性判别分析（LDA）来估算每个组分（物种）丰度对差异效果影响的大小。能够在组与组之间寻找具有统计学差异的生物标识（Biomarker），即组间差异显著的物种。

LDA值分布柱状图：

通过LDA分析（线性回归分析）统计不同组别中有显著作用微生物类群的LDA分值，展示了LDA 分值大于设定阈值的物种，即具有统计学差异的biomaker。柱状图的长度代表显著差异物种的影响大小；

1. #第一步：需要把丰度信息的表格转化成LEfSe识别的格式。这一步会生成.in结尾的文件

2. format_input.py inputfile inputfile.in -c 1 -s -1 -u 2 -o 1000000

3. #第二步：统计显著差异的biomarker、统计子组组间差异、统计effect sizes（LDA score），生成.res格式的文件

4. run_lefse.py -l 2 inputfile.in inputfile.res

5. #第三步：基于第二步的数据，生成多种输出文件，例如png和pdf。

6. plot_res.py --dpi 400 inputfile.res inputfile.png

7. #进化分支图

8. plot_cladogram.py --dpi 400 inputfile.res inputfile.cladogram.png --format png

进化分支图

由内至外辐射的圆圈代表了由门至属（或种）的分类级别。在不同分类级别上的每一个小圆圈代表该水平下的一个分类，小圆圈直径大小与相对丰度大小呈正比。

着色原则：无显著差异的物种统一着色为黄色，差异物种 Biomarker跟随组进行着色，红色节点表示在红色组别中起到重要作用的微生物类群，绿色节点表示在绿色组别中起到重要作用的微生物类群，其它圈颜色意义类同。图中英文字母表示的物种名称在右侧图例中进行展示。

结果怎么看

1. 先看报告中数据量是多少，有多少条序列，一般至少需要1万条序列进行后续分析。看看得到多少个OTU，以及丰度前10的微生物是什么以及在不同分组中丰度的变化。

2. 在OTU 分析文件中看各微生物丰度文件，一般里面有门、纲、目、科、属水平的丰度，目前主流看门水平和属水平。

3. 看Alpha多样性中稀释曲线、Shannon或Rank abundance曲线是否饱和，goods_coverage是否98%以上（数据量是否够）。看看每个样本的多样性 Alpha多样性，四大指数（Shannon、simpson、ace、chao1），看看组间是否有多样性差异，一般指数间变化是一致的，如果不一致可以选择其中的一两种作为结果。

4. Beta多样性中PCA、PCOA、NMDS降维分析的图，看看不同分组是否能分开，如果肉眼能看到明显的分离，说明组间的差异很大。很多情况下，不一定有明显的分离，这时需要看看adonis或者anosim的P值，看组间是否有统计学的差异。

5. 看看LEfSe和Metastats分析出来的显著差异的物种，注意：Metastats分析是两两组间比较，LEfSe分析是多组间差异的比较。

6. 通常研究是的健康人和患病人群的粪便样本，解析患病人群肠道微生物的变化情况。其中最关键的是找到包括有益菌和有害菌在内的差异物种。如果样品的数目不够，会使得组间的统计学差异不够显著，这个时候看趋势也可以。

常见问题

每组测多少个样本够分析？

样本的重复对于测序的实验设计以及后续信息分析都非常重要，设置多个样品重复主要有以下几个作用：

1）能够消除组内误差、检测离群样本：异常样本的存在，会严重影响测序结果的准确性，通过后续信息分析可以发现异常样本，将其排除。

2）增强结果的可靠性。测序的样本越多，越能够降低背景差异，从而增加结果可信度；

一般情况下，个体之间差别较小的模式动物（比如大鼠，小鼠、鸡等）建议每组至少5个生物学重复，一般推荐10个生物学重复。

若是人类肠道粪便的样品，由于个体之间差别较大（比如环境、饮食、遗传条件、健康状态、年龄、性别等影响），建议加大取样量，每组不少于30个样品（取样太少没有说服力，统计学差异也可能不明显）。如果将多个样本混在一起建库测序，多个样合成了一个样，仍然需要多个重复样品。

一般建议多做几个生物学重复，若是遇到个体之间差别较大或者样品处理中出现意外等情况可以选择剔除个别样本，避免后期样本数量不够，如果补送则耽误时间，而且不是同一批的数据会导致重复性不好的情况发生。

16S rRNA基因测序还是16S rDNA测序

文献中存在16S rRNA基因测序和16S rDNA测序两种说法，实际上这两种说法都可以。

16S中的"S"是一个沉降系数，亦即反映生物大分子在离心场中向下沉降速度的一个指标，值越高，说明分子越大。rDNA 和rRNA中的小写字母"r"是ribosome（核糖体）的缩写。

rDNA指的是基因组中编码核糖体RNA（rRNA）分子的对应的DNA序列，也就是编码16S rRNA的基因。

rRNA指的是rDNA的转录产物，它是构成核糖体的重要成分，核糖体由许多小的rRNA分子组装而成，16S rRNA是其中一个元件。

所以我们研究16S菌群测序的对象都是DNA，而不是RNA。

绝对丰度和相对丰度

目前的16S rRNA扩增⼦测序只能检测微生物群落的相对丰度，不能得到绝对的丰度。

⽽只考虑相对丰度会带来很多问题，例如两个相对丰度相同的样本绝对丰度可能差别很⼤。

基于此问题，有方法是在DNA提取步骤之前加⼊一种群落中不存在的细菌作为内标 ，然后利用内标的绝对和相对丰度计算其他物种的绝对丰度。 或是向样品中加入合成的嵌合体DNA作为内参，可定量计算环境样品中微生物的绝对丰度，并在不同来源样品间进行比对。

不同统计方法的结果不一致怎么办？

由于不同的统计分析软件所使用的统计检验的方法有所不同，因此得出的结果也会存在差异。例如LEfSe使用了秩和检验和线性判别分析（LDA），如果某个菌的丰度通过Metastats算出来有差异，而LEfSe算出来没有差异，这时可以根据自己的研究背景选择最为符合预期的分析结果，这些统计分析方法筛选结果都是可信的。

选哪个16S数据库

注释数据库主要有SILVA、RDP、Greengenes。

因为SILVA数据库更新比较及时，因此是目前rRNA基因高通量测序后最常选用的参考数据库之一。此外，SILVA也可被用于平时菌种鉴定时，对rRNA基因sanger测序后的物种进行分类鉴定，此时主要用其SINA Alignment Service功能（https://www.arb-silva.de/aligner/），可非常方便地确定某条rRNA基因序列从门到属/种水平的分类信息并给出各分类水平相应的置信度。

RDP数据库全称“RibosomalDatabase Project”，该数据库提供质控、比对、注释的细菌、古菌16S rRNA基因和真菌28S rRNA基因序列。主要用其Classifier功能（http://rdp.cme.msu.edu/classifier/classifier.jsp），可非常方便地确定某条rRNA基因序列从门到属/种水平的分类信息并给出各水平相应的置信度。

Greengenes是数据库目前没有更新，停留在2013年5月更新的gg_13_5版本，而预测微生物群落功能的分析软件PICRUST，是基于gg_13_5数据库开发的，因此，想做PICRUST分析也必须依托Greengenes的gg_13_5数据库进行比对。

参考文献

[1] Tanja Magoč, Steven L. Salzberg. (2011).FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics (2011) 27 (21): 2957-2963.

[2] Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. NatureMethods, 7: 335–336.

[3] Edgar RC1, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics, 2011 Aug 15;27(16):2194-200.

[4] Edgar, R.C. (2010) Search and clustering orders of magnitude faster than BLAST, Bioinformatics, 26, 2460–2461.

[5]Wang, Q., Garrity, G. M., Tiedje, J. M., Cole, J. R. (2007). Naïve bayesian classifier for rapid assignment of rrnasequences into the new bacterial taxonomy. Applied and environmental microbiology, 73:5261-5267.

[6] Ondov BD, Bergman NH, and Phillippy AM. Interactive metagenomic visualization in a Web browser. BMC Bioinformatics. 2011 Sep 30; 12(1):385.

[7] Chao, A. 1984. Non-parametric estimation of the number of classes in a population. Scandinavian Journal of Statistics 11, 265-270.

[8] Shannon, C. E. and Weaver W. (1948) A mathematical theory of communication. The Bell System Technical Journal, 27, 379–423 and 623–656.

[9] Simpson, E. H. (1949). "Measurement of diversity". Nature 163: 688.

[10] Lozupone, C., Knight, R.(2005).UniFrac: a new phylogenetic method for comparing microbial communities. Applied and environmental microbiology, 71:8228-8235.

---

感谢您的支持，欢迎点击“好看"和转发！！长按下方二维码，即可关注 “基因的生物信息学分析”。