微生物组学研究-16s测序从建库到数据分析全流程解析

2022-12-12 12:12 作者:请叫我小鸡崽 0人读过 | 我要投稿

微生物组学研究-16s测序从建库到数据分析全流程解析

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

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

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

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

5）以及这些（个）菌是否与临床指标具有相关性（对应于菌群与临床指标的相关性分析）；也会进一步想到，既然不同组的微生物群落存在差异，又与疾病具有相关性，

6）那么这些菌群是如何影响宿主的，可能参与了哪些代谢途径（对应于菌群基因功能预测）；

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

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

ITS: 鉴定种以下的水平

18S：鉴定种和种以上的水平

16S：鉴定属水平，少量可以达到种水平

测序物种的差异：16S测序可以对目标区域（可变区域）进行扩增，可以鉴定低丰度物种种类，结果高于宏基因组的结果。
功能分析的差异：16S 只能预测代谢通路，得不到物种对代谢通路的信息；宏基因组是全基因组的测序，可以得到具体的代谢通路的情况

特定时空下，并不是所有的基因都表达
宏转录组测序：特定时空下，有活性的微生物群落的组成

16S 测序测的是编码16S rRNA的DNA序列

5S: 碱基少，不足以反映菌群的差异性
23S：碱基数多，变异多，不容易区分亲缘关系
16S：对保守区域设计探针，然后扩增可变区域；碱基数量少，变异少

测序结果的reads直接进行97%相似程度归为划分小组，一个小组为一个OTU
一个OUT可以注释到一个物种，一个物种可以被多个OUT注释

#### OUT解释

OTU(operational taxonomic units) 是在系统发生学研究或群体遗传学研究中，为了便于进行分析，人为给某一个分类单元（品系，种，属，分组等）设置的同一标志。通常按照 97% 的相似性阈值将序列划分为不同的 OTU，每一个 OTU 通常被视为一个微生物物种。相似性小于97%就可以认为属于不同的种，相似性小于93%-95%，可以认为属于不同的属。样品中的微生物多样性和不同微生物的丰度都是基于对OTU的分析。

下表是对每个样本在分类字水平上的数量进行统计，并且在表栺中列出了在每个分类字水平上的物种数目

其中SampleName表示样本名称；Phylum表示分类到门的OTU数量；Class表示分类到纲的OTU数量；Order表示分类到目的OTU数量；Family表示分类到科的OTU数量；Genus表示分类到属的OTU数量；Species表示分类到种的OTU数量。

α多样性关注【样本自身】的菌群丰富度和均匀度，而β多样性关注【样本间】的菌群组成与分布的差异。只有当样本（组）间菌群组成存在差异，才有可能进一步探讨菌群失调与疾病的关系。
α多样性是度量【单个样本内】有多少种微生物物种，以及每个物种所占比例的指标。
β多样性是度量【不同样本间】菌群组成的相似度大小的指标，即关注各样本间的菌群组成差异
a多样性：样本内的比较（单个样本的比较）
β多样性：样本间的比较（样本间的比较）

Alpha多样性（样本内多样性）

Alpha多样性是指一个特定区域或者生态系统内的多样性，常用的度量指标有Chao1 丰富度估计量（Chao1 richness estimator）、香农 - 威纳多样性指数（Shannon-wiener diversity index）、辛普森多样性指数（Simpson diversity index）等。

计算菌群丰度：Chao、ace；

计算菌群多样性：Shannon、Simpson。

Simpson指数值越大，说明群落多样性越高；Shannon指数越大，说明群落多样性越高。

Chao1：是用chao1 算法计算群落中只检测到1次和2次的OTU数估计群落中实际存在的物种数。Chao1 在生态学中常用来估计物种总数，由Chao (1984) 最早提出。Chao1值越大代表物种总数越多。

Schao1=Sobs+n1(n1-1)/2(n2+1)

其中Schao1为估计的OTU数，Sobs为观测到的OTU数，n1为只有一条序列的OTU数目，n2为只有两条序列的OTU数目。

Shannon：用来估算样品中微生物的多样性指数之一。它与 Simpson 多样性指数均为常用的反映 alpha 多样性的指数。Shannon值越大，说明群落多样性越高。

Ace：用来估计群落中含有OTU 数目的指数，由Chao 提出，是生态学中估计物种总数的常用指数之一，与Chao1 的算法不同。

Simpson：用来估算样品中微生物的多样性指数之一，由Edward Hugh Simpson ( 1949) 提出，在生态学中常用来定量的描述一个区域的生物多样性。Simpson 指数值越大，说明群落多样性越高。

Beta多样性分析（样品间差异分析）

也许我们有听说Beta多样性在最近10年间成为生物多样性研究的热点问题之一。具体解释下：

Beta多样性度量时空尺度上物种组成的变化, 是生物多样性的重要组成部分, 与许多生态学和进化生物学问题密切相关！

PCoA分析

PCoA（principal co-ordinates analysis）是一种研究数据相似性或差异性的可视化方法，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，PCoA 可以找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样品点之间的相互位置关系，只是改变了坐标系统。

重要的是，它是可以用来观察个体或群体间的差异的。

PCA分析

主成分分析（Principal component analysis）PCA 是一种研究数据相似性或差异性的可视化方法，通过一系列的特征值和特征向量进行排序后，选择主要的前几位特征值，采取降维的思想，PCA 可以找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样品点之间的相互位置关系，只是改变了坐标系统。

显著差异菌群分析

通过β多样性分析，可以确定不同组间的微生物群落是存在差异的，接着需要进一步找出哪些菌（群）引起了组间的群落差异。只有找出核心菌（群），才能明确下一步的研究方向。在报告中，使用目前在文献中高频出现的方法——LEfSe（Linear discriminant analysis Effect Size），来做菌群差异分析，寻找生物标志物（Biomarker）。该方法综合了统计学上的差异分析和该差异物种对分组结果的影响力得分值，同时强调了统计学意义和生物相关性。

LEFse原理

Step1. 首先在多组样本中采用 非参数因子Kruskal-Wallis秩和检验 检测不同分组间丰度差异显著的物种；也就是图中按class1 和class2两个大的分组，每一行都进行检验，初步得到差异物种，通过检验的打钩进入step2检验；

Step2. 再利用Wilcoxon秩和检验，对每一组中的亚组进行两两检验，具有显著差异的再进行下一轮检验。

Step3. 最后用线性判别分析（LDA）对数据进行降维并评估差异显著的物种的影响力（即LDA score）。

前两步的Kruskal-Wallis秩和检验、Wilcoxon秩和检验 比较简单，类似T检验或者方差检验等，只不过T检验和方差分析为参数检验（要求数据符合方差齐性、正态分布），而在微生物多样性分析中，样品物种丰度分布不确定，多采用非参数检验，所以采用非参数的Kruskal-Wallis秩和检验、Wilcoxon秩和检验。比较复杂一点的就是最后的LDA分析。

LDA是一种监督学习的降维技术，也就是说其数据集中的每个样本是有类别输出的。是在目前机器学习、数据挖掘领域经典且热门的一个算法这点和PCA不同。PCA是不考虑样本类别输出的无监督降维技术。LDA是有监督的，所以LDA算法可以很好的利用样本的分组信息，得到的结果更可靠，这就是LDA分析优势。理解了LDA分析的原理，就不难理解LEfSe的分析结果了。

LDA的全称是Linear Discriminant Analysis （线性判别分析），是一种supervised learning （有监督学习）。有些资料上也称为是Fisher' s Linear Discriminant，由Ronald Fisher发明自1936年，是在目前机器学习、数据挖掘领域经典且热的一个算法。

LDA的思想可以用一句话概括，就是”投影后类内方差最小，类间方差最大”。简单来说就是一种投影，是将一个高维的点投影到一个低维空间，我们希望映射之后，不同类别之间的距离越远越好，同类别之中的距离越近越好。

是不是很抽象哇，来举个栗子吧。假设我们有两类数据：分别为红色和蓝色，如下图所际，这些数据特征是二维的，我们希望将这些数据投影到一维的一条直线，让每一种类别数据的投影点尽可能的接近，而红色和蓝色数据中心之间的距离尽可能的大。

从直观上可以看出，右图要比左图的投影效果好，因为右图的红色数据和色数据各个较为集中，类别之间的距离明显。左图则在边界处数据混杂。当然在实际应用中，我们的数据多个类别的，我们的原始数据一般也是超过二维的，投影后的也一般不是直线，而是一-个低维的超平面。

菌群标志物预测能力评估

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

上图以灵敏度为纵坐标，特异度为横坐标绘制曲线。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不符合真实情况，在实际中极少出现。

菌群基因功能预测

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