微生物16S测序数据的正确打开方式

微生物16S测序数据的正确打开方式16S rRNA基因测序（也称16S rDNA测序）是最常用的菌群多样性分析的手段。对于新手，如果收到一份不讲“人话”的16S测序分析报告，很快就会被各种生态学术语、各种指数、各种分析方法弄晕。

7个问题串起16S测序的核心结果

怎么办？用你的研究逻辑来梳理16S测序数据（图1）。

简单地说，做16S测序是为了鉴定样本中的微生物（细菌）群组成，找微生物群与疾病或表型的相关性。

详细地说，

1）首先想了解在不同组样本中各有哪些微生物存在和丰富度（对应于菌群鉴定和α多样性分析）；

2）接着想看不同样本组间微生物群组成是否存在差异（对应于β多样性分析）；

3）如果是，那么就有必要找出引起不同组样本微生物群差异的关键菌。如果不是，那说明微生物群比如肠道菌群与疾病或表型可能并不相关（基于已有的研究，这种可能性

比较小）；

4）找到了关键菌，在临床上，很自然会想到，这些（个）关键菌是否可以作为Biomarker（对应于疾病诊断模型构建），比如用于区分糖尿病前期患者与健康组的标志物；

5）以及这些（个）菌是否与临床指标具有相关性（对应于菌群与临床指标的相关性分析）；也会进一步想到，既然不同组的微生物群落存在差异，又与疾病具有相关性，6）那么这些菌群是如何影响宿主的，可能参与了哪些代谢途径（对应于菌群基因功能预测）；

7）这些预测到的菌群功能是否与疾病有关，通常是肯定的。最后把这些结果整合起来分析，可以初步得出菌群组成的变化是如何与疾病或表型相关的。

顺着上述7个生物学问题来看16S测序结果，你会轻松拨开迷雾，直达核心结果。

图1 7个问题串起16S测序的核心结果

6张图就够发菌群与疾病相关性文章

编者对2019发表的数十篇以16S测序为主的肠道菌群与疾病关系研究文章（IF 5至10分）的内容进行了分析和归纳，发现大部分文章的Results部分都是由图1所列的核

心结果组成。以联川生物医学16S测序报告为例，具体讲解16S测序文章中的核心结果及其对应的图表。

（1）菌群鉴定与物种分布

采用最新的QIIME 2分析流程，并使用更严谨的DADA2算法对扩增子数据进行去噪，相当于以100%的相似性聚类（取代传统的OTU聚类），仅对低质量序列进行去除和校正等，获取扩增子序列变异，然后去冗余，即得到feature（特征）数据。将feature数据和16S数据库（如SILVA、NT-16S）进行序列比对，可以对样本中检测到的细菌从界（Kingdom）、门（Phylum）、纲（Class）、目（Order）、科（Family）、属（Genus）、种（Species）多个分类学层级进行物种鉴定和注释。然后根据各个分类层级上的物种相对丰度以物种分布堆叠图（图2A-C，以门水平为例）来直观展示。

图2A 物种分布堆叠图

图2B 物种聚类堆叠图

图2C 物种相对丰度Heatmap

在联川医学16S测序报告中，会提供上述三种主流的物种分布堆叠图，你可以选择其一使用。在图2A、B中，不同颜色的柱子对应不同的物种，柱子的长短代表该物种所占比

例的大小。图2B中左侧采用Bray-Curtis距离法分析样本间菌群组成的相似性并进行聚类。图2C中展示了不同细菌物种在不同样本中的相对丰度情况，颜色越红，丰度越高，颜色越蓝，丰度越低。

（2）菌群α多样性

α多样性是度量单个样本内有多少种微生物物种，以及每个物种所占比例的指标。在报告中，采用5种常用指数来度量α多样性：Observed species和Chao1反映样本中物种丰富度，但不考虑每个物种的占比情况（均匀度）；Shannon 和Simpson反映物种的丰富度和均匀度；Good’s Coverage 反映样本的测序深度。使用Wilcoxon秩和检验对上述各个指数的样本数据进行分析，筛选出各样本组比较中显著差异的α多样性指数并绘制小提琴图（图3）。

图3 小提琴图表示α多样性指数

小提琴图集合了箱形图和密度图的特征。上图以Good’s Coverage为例，左上角给出了差异分析使用的检验方法和计算得到的p值。当p<0.01，表示差异极显著；当p<0.05，表示差异显著；当p>0.05，则表示无显著性差异。

（3）菌群β多样性

β多样性是度量不同样本间菌群组成的相似度大小的指标，即关注各样本间的菌群组成差异。α多样性关注样本

自身的菌群丰富度和均匀度，而β多样性关注样本间的菌群组成与分布的差异。只有当样本（组）间菌群组成存在差异，才有可能进一步探讨菌群失调与疾病的关系。在报告中，采用主流的PCA、PCoA、NMDS、ANOSIM、Adonis、UPGMA等多种分析方法来考察和区分样本间的菌群组成差异（图4，以PCoA为例）。

图4 PCoA分析结果图

上图中每一个点代表一个样本，相同颜色的点来自同一个分组，两点之间距离越近表明两者的群落构成差异越小。左图是基于Unweighted UniFrac的PCoA分析结果，右图是基于Weighted UniFrac的PCoA分析结果。在这个例子中，采用Weighted UniFrac的PCoA分析更能把不同组的样本区分开来，且p值<0.01，具有显著统计学差异。需要说明的是，PCoA分析本身是没有p值计算的，p值来自于ANOSIM 分析的结果。在绘图时，把p值加入了PCoA图中。

由于每个项目的实验设计和样本菌群组成差异巨大，无法预先知道哪种β多样性分析方法是将样本间菌群差异区分开的最优方法。因此，在报告中提供了多种β多样性分析方法和产生的图片，在撰写文章时，你只需要从中选出最能解释生物学问题的图片用在文章中即可（通常展示一个或者两个β多样性分析结果）。

（4）显著差异菌群分析

通过β多样性分析，可以确定不同组间的微生物群落是存在差异的，接着需要进一步找出哪些菌（群）引起了组间的群落差异。只有找出核心菌（群），才能明确下一步的研究方向。在报告中，使用目前在文献中高频出现的方法——LEfSe（Linear discriminant analysis Effect Size），来做菌群差异分析，寻找生物标志物（Biomarker）。该方法综合了统计学上的差异分析和该差异物种对分组结果的影响力得分值，同时强调了统计学意义和生物相关性。LEfSe分析结果图，通常包括进化分支图（图5A）和LDA值分布柱状图（图5B）。需要说明的是，联川不仅提供LEfSe筛选差异菌群，还提供其他多种方法，如随机森林分析（图9）、秩和检验等。

图5A LDA值分布柱状图

上图主要展示了LDA score大于预设值的显著差异物种（less_strict设为2；more_strict 设为4），即具有统计学差异的Biomarker；柱状图的颜色代表各自的组别，长短代表的是LDA score，即不同组间显著差异物种的影响程度。

图5B 进化分支图

上图中，小圆圈: 图中由内至外辐射的圆圈代表了由界

（单个圆圈）至属（或种）的分类级别。不同分类级别上的每一个小圆圈代表该水平下的一个分类，小圆圈直径大小与相对丰度大小呈正比。颜色：无显著差异的物种统一上黄色，差异显著的物种Biomarker跟随组别进行上色，红色节点表示在红色组别中起到重要作用的微生物类群，绿色节点表示在绿色组别中起到重要作用的微生物类群。未能在图中显示的Biomarker对应的物种名会展示在右侧，字母编号与图中对应。

（5）菌群标志物预测能力评估

受试者工作特征（ROC）曲线分析是一种常用的统计学分析方法，在医学研究中主要用于评价诊断试验的效能。在报告中，通过绘制ROC曲线，并计算ROC曲线下面积（AUC），来确定哪种菌（群）具有最佳的诊断价值（图6）。

图6 菌群标志物ROC曲线分析

上图以灵敏度为纵坐标，特异度为横坐标绘制曲线。ROC 曲线越靠近左上角，试验的准确性就越高。若AUC值为1.0，反映出对两个群组的完美区分，且不存在预测误差。若AUC 值在1.0和0.5之间，在AUC>0.5的情况下，AUC越接近于1，说明诊断效果越好。AUC在0.5~0.7时有较低准确性，AUC 在0.7~0.9时有一定准确性，AUC在0.9以上时有较高准确性。AUC=0.5时，说明诊断方法完全不起作用，无诊断价值。

AUC<0.5不符合真实情况，在实际中极少出现。

（6）菌群基因功能预测

因为菌群功能预测软件PICRUSt（Phylogenetic Investigation of Communities by Reconstruction of Unobserved States）的出现，研究者能进一步基于16S测序数据预测菌群可能参与的代谢通路（尽管并没有测定菌群基因信息），以便能初步讨论菌群组成变化与疾病或表型是如何关联在一起的。在联川报告中，使用最新的PICRUSt 2，相比上一版，用于预测的参考基因组数据库已扩展超过10倍，可以获得包括COG，EC，KO，PFAM，TIGRFAM等数据库对菌群的基因功能注释结果。然后，再使用STAMP软件进行差异分析，得到在不同样本组中显著差异的菌群基因功能（图7，以pathway结果为例）。如果要系统研究菌群携带的基因及其功能，则应该做宏基因组测序。

图7 PICRUSt 2预测菌群基因功能

上图中比较了不同组菌群的KEGG pathway，并筛选出具有显著性组间差异的 pathway。左边柱状图代表某代谢通路的丰度分别占两组样本中所有代谢通路的百分比，右边为corrected p值。

至此，一篇医学微生物组16S测序文章的主要结果和图表就基本齐备了。当然，完整的医学16S测序报告还包括更

多内容（图8），而且16S测序数据还有许多扩展性以及个性化的分析图表（图9），联川会根据研究者的具体需求来提供。

图8 联川生物医学16S测序报告内容

图9 16S测序数据个性化分析图表（部分）

值得参考的菌群与疾病关系的研究思路

甲状腺癌患者的肠道菌群和代谢谱的变化（IF4.982）

发表期刊：International Journal of Cancer

影响因子：4.982

发表时间：2019

研究内容：肠道菌群、代谢谱与甲状腺癌之间的关系

样本数量：16S测序：30例甲状腺癌（TC）vs 35例健康对照（HCs）;代谢组：15例TC vs 15例HCs

样本类型：粪便样本

实验方法：16S rRNA基因测序+非靶向代谢组检测

实验设计：

所有受试者均为汉族，出生于中国东北地区，饮食结构相似

研究思路：

采用目前主流的肠道菌群与疾病关系的研究策略，联合肠道菌群16S测序和粪便代谢组（极具性价比的检测组合），并结合临床指标来一起讨论，给出甲状腺癌患者的肠道菌群特征和代谢谱，潜在的疾病标志物，以及肠道菌群影响肿瘤发生发展的潜在途径。

研究结果：

本研究对TC患者的肠道菌群及其代谢产物进行了全面的研究。鉴定出21个细菌属和72个代谢产物发生显著变化，并找到8种代谢产物结合5个细菌属对TC与HCs的鉴别更为有效（AUC=0.97）。