一篇医学微生物组16S测序文章必备的六张图 | 微生物专题

上期介绍了《一篇医学微生物组16S测序文章必备的五个内容 | 微生物专题》，接下来我们逐一介绍这些必备内容需要用哪些图片来直观呈现，以及这些图片的含义。

什么是Feature

在开始前，我们先要引出一个重要的新名词——feature。说到feature，就必须提到它的前辈OTU了。OTU（Operational Taxonomic Units，操作分类单元）是过去在微生物群落多样性研究中非常重要的一个概念，是一种将复杂测序数据降维成简单数据的做法。通常按照97%的相似性阈值将序列聚类为不同的 OTU，每一个OTU通常被视为一个微生物物种。然而，这种聚类方法已显露出许多不足，比如偏好性，序列可变性等，也使得不同研究的数据无法合并或难以合并起来做大数据的整合分析。近来研究大多倾向于采用100%相似度聚类取代传统的97%的相似性聚类，以达到更准确的菌群鉴定和多样性分析。在最新的QIIME 2分析流程中，采用DADA2（或Deblur）算法对扩增子数据进行去噪，相当于以100%的相似度聚类，仅对低质量序列进行去除和校正等，获取扩增子序列变异，然后去冗余，即得到feature（特征）数据。后续扩增子测序数据分析都是基于feature数据展开。如果你对OTU很熟悉，其实把feature看做OTU来理解也无妨。我们的扩增子测序分析已全面基于feature方法进行（发微生物组16S测序高分文章，四招让你事半功倍 | 微生物专题）。

图1 feature取代OTU成为扩增子测序数据分析的趋势

言归正传，我们来介绍与医学微生物组16S测序文章必备内容对应的图片。

1. 菌群鉴定与物种分布

使用feature数据，去和16S数据库SILVA和NT-16S进行序列比对，就可以对样本中检测到的16S序列从界（Kingdom）、门（Phylum）、纲（Class）、目（Order）、科（Family）、属（Genus）、种（Species）多个分类学层级进行物种鉴定和注释。然后根据各个分类层级上的物种相对丰度来绘制第一张必备的图片——物种分布堆叠图

Figure 1-1物种分布堆叠图

Figure 1-2 物种聚类堆叠图

Figure 1-3物种相对丰度Heatmap

在医学16S测序报告中，我们会提供三种主流的物种分布堆叠图（图2-1、2-2、2-3，以门水平为例），你可以选择其一使用。在图2-1、2-2中，不同颜色的柱子对应不同的物种，柱子的长短代表该物种所占比例的大小。图2-2中左侧采用Bray-Curtis距离法分析样本间菌群组成的相似性并进行聚类。图2-3中展示了不同细菌物种在不同样本中的相对丰度情况，颜色越红，丰度越高，颜色越蓝，丰度越低。

2. 菌群α多样性

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

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

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

3. 菌群β多样性

β多样性是度量不同样本间菌群组成的相似度大小，即关注各样本间的菌群组成差异。α多样性关注样本自身的菌群丰富度和均匀度，而β多样性关注样本间的菌群组成与分布的差异。只有当样本（或组）间菌群组成存在差异，我们才有可能进一步探讨菌群失调与疾病的关系。在医学16S测序报告中，我们采用主流的PCA、PCoA、NMDS、ANOSIM、Adonis、UPGMA等多种分析方法来考察和区分样本间的菌群组成差异（图4，以较常使用的PCoA为例）。

Figure 3 PCoA分析结果图

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

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

4. 显著差异菌群分析

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

Figure 4-1 LDA值分布柱状图

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

Figure 4-2 进化分支图

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

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

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

Figure 5 菌群标志物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的出现，我们能进一步对16S测序数据进行挖掘，预测菌群基因可能携带的功能（尽管并没有测定菌群基因组信息），以便我们能初步讨论菌群失调与疾病是如何关联在一起的。在医学16S测序报告中，我们使用最新的PICRUSt 2，相比上一版，用于预测的参考基因组数据库已扩展超过10倍，可以获得包括COG，EC，KO，PFAM，TIGRFAM等数据库对菌群的基因功能注释结果。然后，再使用STAMP进行差异分析，得到在不同样本组中显著差异的菌群基因功能（图7，以pathway结果为例）。

Figure 6 PICRUSt 2预测菌群基因功能

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

有了这6张必备的figures，一篇医学微生物组16S测序文章的主要结果就齐备了。不仅如此，我们还为大家准备了更多为文章加分添彩的图片，部分展示如下。

关于联川医学16S测序报告介绍，请查看《医学16S测序报告：为忙碌的医生而准备 | 微生物》 。

相关阅读