16S测序分析（二）菌群多样性分析
导读
多样性分析是16S测序分析中常见的分析方法。本文旨在向初学者介绍多样性分析中alpha多样性和beta多样性的由来、概念、计算，以及展示“如何用R语言实现alpha多样性和beta多样性的计算和可视化”的工作流程。
一、概念介绍
Alpha多样性、beta多样性和gamma多样性的概念由美国杰出的生态学家Robert Harding Whittaker提出。Whittaker将一个环境中总的物种多样性命名为gamma多样性。Gamma多样性由alpha多样性和beta多样性共同决定。Alpha多样性指在环境中的一个区域内的平均物种多样性。Beta多样性指环境中不同区域之间的差异。
二、Alpha多样性指数
1. Shannon指数：
Shannon指数最初是由Claude Shannon提出用来计算字符串文本信息熵的指标，后来逐渐发展成生态学研究中最常用的多样性指数。Shannon指数不只考虑物种丰富度（Richness，样本中物种数），而且同时考虑物种的均匀度（Evenness，平均程度），所以它是反应群落结构的综合指标。
计算方法：
公式
Pi：第i个物种的个体数占样本中总个体数的比例
R：样本中的总物种数
样本的Shannon指数越高，则其物种多样性也越高。
2. Simpson指数：
Simpson指数是由Edward H. Simpson在1949年提出评价生物多样性的另一种常用指标。它也是既考虑样本中的物种丰富度，又考虑物种均匀度的综合指标。
计算方法：
公式
Pi：第i个物种的个体数占样本中总个体数的比例
R：样本中的总物种数
Simpson指数的取值范围在0和1之间。样本的Simpson指数越高，则其物种多样性也越高。
三、Beta多样性的计算
Beta多样性专用于不同样本间的比较，它不能直接通过某一个样本的物种丰富度和均匀度计算出该样本的多样性度量值。Beta多样性是利用不同样本间的丰度变化或进化关系来计算样本间距离，从而反映样本间是否具有显著的微生物群落差异。
计算beta多样性的方法有很多：有最为常用的bray curtis距离、Jaccard距离还有欧式距离，他们考虑的是样本间物种丰度（有无）和均度（相对丰度）；另外还有Unifrac距离法，它是根据系统发生树进行比较。
四、准备工作
软件准备（window环境）：
1. R
地址：https://www.r-project.org/
版本：3.4.1
2. R包：openxlsx
功能：打开Excel文件
版本：4.1
3. R包：vegan
功能：包内含有多种生态学分析必备的函数
版本：2.4
4. Excel
功能：数据处理
数据准备：
标准化菌属相对丰度表。获取方法请参考16S测序分析系列（一）菌属丰度表获取。
五、Alpha多样性的计算和可视化
install.packages("openxlsx") ## 下载安装openxlsx包
install.packages("vegan") ## 下载安装vegan包
setwd("C:/mywd") ## 设置工作目录为C盘mywd文件夹
getwd() ## 进入工作目录
library(openxlsx) ## 调用openxlsx包
library(vegan) ## 调用vegan包
workbook <- "C:/mywd/data.xlsx"
## 读取并存储data.xlsx的内容到workbook
说明：data.xlsx必须是标准化后的菌属相对丰度表，格式如下：第一行是菌属名称，其余每行代表一个样本，每列代表一类菌属的在每一个样本中的相对丰度。
图片.png
mydataframe <- read.xlsx(workbook, 1)
## 读取并存储workbook中的第一个sheet的内容到mydataframe
alpha_shannon <- diversity(mydataframe, "shannon")
## 用vegan包中的diversity函数计算shannon多样性
write.xlsx(alpha_shannon, "data_diversity.xlsx")
## 获得shannon多样性指标的计算结果，即data_diversity.xlsx
## 保存alpha_shannon到当前的data_diversity.xlsx中（可自动创建）
数据处理：打开data_diversity.xlsx文件如“1”。其中第一个样本的shannon index是shannon index 1，以此类推。然后我们把这些shannon index加上变量名和分组信息（将其分为A，B两组），得到“2”。利用R或者Excel可以对“2”中的数据进行分组T检验，方法很简单，所以这里不再加以展示。接下来我将利用“2”制作boxplot来更直观的观察两组的shannon index是否有明显的差异。
图片.png
workbook2 <- "C:/mywd/data_diversity.xlsx"
## 读取并存储data_diversity.xlsx的内容到workbook2
mydataframe2 <- read.xlsx(workbook2, 1)
## 读取并存储workbook2中的第一个sheet的内容到mydataframe2
boxplot(mydataframe2$Diversity~ mydataframe2$Group, main="Shannon index", col=c("green", "red"))
## mian：命名；col：上色
## 可视化结果如下。
图片.png
由图可见，A组的多样性明显高于B组。
六、Beta多样性的计算和可视化
setwd("C:/mywd")
getwd()
library(openxlsx)
library(vegan)
workbook <- "C:/mywd/data.xlsx"
mydataframe <- read.xlsx(workbook, 1)
ord <- cmdscale(vegdist(mydataframe, method="bray"))
write.xlsx(ord, "data_ord.xlsx")
数据处理：data_ord.xlsx中保存的是每个样本的二维坐标信息，它是进行此分析的最关键的文件。用windows打开data_ord.xlsx文件，格式如“1”。第一行是变量名，其余每行是每个样本的坐标信息。添加分组信息，如下图“2”。为了方便作图，我这里直接使用颜色”Color”代替组别”Group”，即red代替Group A，blue代替Group B。接下来利用“2”中的样本坐标信息进行绘图。
图片.png
workbook3 <- "C:/mywd/data_ord.xlsx"
mydataframe3 <- read.xlsx(workbook3, 1)
plot(x=mydataframe3$V1,y=mydataframe3$V2,xlab="mds1",ylab="mds2",main="Beta diversity ( Bray Curtis )" , pch=19, col=mydataframe3$Color)
## x=mydataframe3$V1,y=mydataframe3$V2：以V1为横坐标，V2为纵坐标
## xlab：设置坐标名称
## main：设置图片名称
## pch=19：用点表示样本
## col=mydataframe3$Color：根据Color中的信息对点上色
legend("topright", pch=c(19, 19), col=c("red", "blue"), legend=c("Group A", "Group B"))
## "topright"：设置图注在右上角
## pch=c(19, 19)：绘制两个点代表两个组
## col：给点上色
## legend：给点命名
图片.png
结束语
到这里就完成了beta多样性的计算和可视化的所有工作。从图中可以看出A和B两个组的样本的肠道菌群有明显的差异，所以可以推测“分组变量”很可能和肠道菌群有联系。要想知道具体是“哪些菌”对此联系做了贡献，还需要进行更深层次的分析。本文到此就结束了，如有疑问欢迎留言交流。下次将为大家推出“如何从肠道菌群中寻找你感兴趣的细菌？”来介绍寻找“哪些菌”的方法。
相关阅读：
16S测序分析（一）菌属丰度表获取
16S测序分析（二）菌群多样性分析
16S测序分析（三）用LEfSe寻找组间差异细菌
16S测序分析（四）用MaAsLin寻找组间差异细菌
16S测序分析（五）用RandomForest寻找关键细菌
16S测序分析（六）用PICRUSt预测菌群KEGG代谢通路
