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

1、微生物16S测序数据的正确打开方式 16S rRNA基因测序（也称16S rDNA测序）是最常用的菌群多样性分析的手段。对于新手，如果收到一份不讲“人话”的16S测序分析报告，很快就会被各种生态学术语、各种指数、各种分析方法弄晕。 7个问题串起16S测序的核心结果 怎么办？用你的研究逻辑来梳理16S测序数据（图1）。 简单地说，做16S测序是为了鉴定样本中的微生物（细菌）群组成，找微生物群与疾病或表型的相关性。 详细地说， 1）首先想了解在不同组样本中各有哪些微生物存在和丰富度（对应于菌群鉴定和多样性分析）； 2）接着想看不同样本组间微生物群组成是否存在差异（对应于多样性分析）； 3）如果是，

2、那么就有必要找出引起不同组样本微生物群差异的关键菌。如果不是，那说明微生物群比如肠道菌群与疾病或表型可能并不相关（基于已有的研究，这种可能性比较小）； 4）找到了关键菌，在临床上，很自然会想到，这些（个）关键菌是否可以作为Biomarker（对应于疾病诊断模型构建），比如用于区分糖尿病前期患者与健康组的标志物； 5）以及这些（个）菌是否与临床指标具有相关性（对应于菌群与临床指标的相关性分析）；也会进一步想到，既然不同组的微生物群落存在差异，又与疾病具有相关性， 6）那么这些菌群是如何影响宿主的，可能参与了哪些代谢途径（对应于菌群基因功能预测）； 7）这些预测到的菌群功能是否与疾病有关，通常是肯

3、定的。最后把这些结果整合起来分析，可以初步得出菌群组成的变化是如何与疾病或表型相关的。 顺着上述7个生物学问题来看16S测序结果，你会轻松拨开迷雾，直达核心结果。 图1 7个问题串起16S测序的核心结果 6张图就够发菌群与疾病相关性文章 编者对2019发表的数十篇以16S测序为主的肠道菌群与疾病关系研究文章（IF 5至10分）的内容进行了分析和所列的核1部分都是由图Results归纳，发现大部分文章的心结果组成。以联川生物医学16S测序报告为例，具体讲解16S测序文章中的核心结果及其对应的图表。 （1）菌群鉴定与物种分布 采用最新的QIIME 2分析流程，并使用更严谨的DADA2算法对扩增子数

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

5、Heatmap 在联川医学16S测序报告中，会提供上述三种主流的物种分布堆叠图，你可以选择其一使用。在图2A、B中，不同颜色的柱子对应不同的物种，柱子的长短代表该物种所占比例的大小。图2B中左侧采用Bray-Curtis距离法分析样本间菌群组成的相似性并进行聚类。图2C中展示了不同细菌物种在不同样本中的相对丰度情况，颜色越红，丰度越高，颜色越蓝，丰度越低。 （2）菌群多样性 多样性是度量单个样本内有多少种微生物物种，以及每个物种所占比例的指标。在报告中，采用5种常用指数来度量多样性：Observed species和Chao1反映样本中物种丰富度，但不考虑每个物种的占比情况（均匀度）；Shan

6、non和Simpson反映物种的丰富度和均匀度；Goods Coverage反映样本的测序深度。使用Wilcoxon秩和检验对上述各个指数的样本数据进行分析，筛选出各样本组比较中显著差异的多样性指数并绘制小提琴图（图3）。 图3 小提琴图表示多样性指数 小提琴图集合了箱形图和密度图的特征。上图以Goods Coverage为例，左上角给出了差异分析使用的检验方法和计算得到的p值。当p0.01，表示差异极显著；当p0.05，则表示无显著性差异。 （3）菌群多样性 多样性是度量不同样本间菌群组成的相似度大小的指标，即关注各样本间的菌群组成差异。多样性关注样本自身的菌群丰富度和均匀度，而多样性关注样

7、本间的菌群组成与分布的差异。只有当样本（组）间菌群组成存在差异，才有可能进一步探讨菌群失调与疾病的关系。在报告中，采用主流的PCA、PCoA、NMDS、ANOSIM、Adonis、UPGMA等多种分析方法来考察和区分样本间的菌群组成差异（图4，以PCoA为例）。 图4 PCoA分析结果图 上图中每一个点代表一个样本，相同颜色的点来自同一个分组，两点之间距离越近表明两者的群落构成差异越小。左图是基于Unweighted UniFrac的PCoA分析结果，右图是基于Weighted UniFrac的PCoA分析结果。在这个例子中，采用Weighted UniFrac的PCoA分析更能把不同组的样本

8、区分开来，且p值0.5的情况下，AUC越接近于1，说明诊断效果越好。AUC在0.50.7时有较低准确性，AUC在0.70.9时有一定准确性，AUC在0.9以上时有较高准确时，说明诊断方法完全不起作用，无诊断价值。AUC=0.5性。AUC0.5不符合真实情况，在实际中极少出现。 （6）菌群基因功能预测 因为菌群功能预测软件PICRUSt（Phylogenetic Investigation of Communities by Reconstruction of Unobserved States）的出现，研究者能进一步基于16S测序数据预测菌群可能参与的代谢通路（尽管并没有测定菌群基因信息），以

9、便能初步讨论菌群组成变化与疾病或表型是如何关联在一起的。在联川报告中，使用最新的PICRUSt 2，相比上一版，用于预测的参考基因组数据库已扩展超过10倍，可以获得包括COG，EC，KO，PFAM，TIGRFAM等数据库对菌群的基因功能注释结果。然后，再使用STAMP软件进行差异分析，得到在不同样本组中显著差异的菌群基因功能（图7，以pathway结果为例）。如果要系统研究菌群携带的基因及其功能，则应该做宏基因组测序。 图7 PICRUSt 2预测菌群基因功能 上图中比较了不同组菌群的KEGG pathway，并筛选出具有显著性组间差异的 pathway。左边柱状图代表某代谢通路的丰度分别占两

10、组样本中所有代谢通路的百分比，右边为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测序和粪便代谢组（极具性价比的检测