技术贴 | R语言菌群Alpha多样性分析和绘图 / 开普饭

本文由阿童木根据实践经验而整理，希望对大家有帮助。

原创微文，欢迎转发转载。

导读

箱型图（Boxplot）或者盒图是一种能同时展示一组或多组数据的极值、四分位数、中位数和离群值，显示数据离散情况的统计图。下面介绍一下如何在R软件中利用箱形图可视化两组微生物群的Alpha多样性指数。过程分为以下几步：1）模拟原始数据；2）模拟分组；3）标准化；4）计算多样性；5）T检验；6）ggplot绘制Boxplot，保存结果。

完成以上步骤将得到如下结果：

1 模拟丰度矩阵

set.seed(1995)

# 设置随机种子：为了能多次获取同一组随机数

data=matrix(abs(round(rnorm(200, mean=1000, sd=500))), 20, 10)

# 获取随机正整数矩阵：20个样本，10个菌种

colnames(data)=paste("Species", 1:10, sep=" ")

# 给细菌命名

rownames(data)=paste("Sample", 1:20, sep=" ")

# 给样本命名，结果如下：

2 模拟分组文件

group=c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B")

# 把样本分为A、B两组

sample_id=rownames(data)

# 取样本名

data_group=data.frame(sample_id, group)

# 做分组文件，如下：

3 标准化丰度

data_norm=data

# 新建data_norm数据框，赋初值

for(i in 1:20){

sample_sum=apply(data, 1, sum)

# 计算每个样本的菌种总丰度

for(j in 1:10){

data_norm[i,j]=data[i,j]/sample_sum[i]

# [每个样本的每个物种的丰度] / [该样本物种总丰度] = 相对丰度

}

apply(data_norm, 1, sum)

# 检查每个样本物种总丰度是否成功标准化到1，如下：

4 计算Alpha多样性

library(vegan)

# 加载VEGAN包

data_norm_shannon=diversity(data_norm, "shannon")

# 使用VEGAN包中的diversity函数计算每个样品的Shannon Alpha多样性指数

data_ggplot=data.frame(data_norm_shannon)

# 多样性计算结果如下：

data_ggplot=data.frame(data_ggplot, data_group["group"])

# 添加分组信息，如下：

5 T检验

成组设计T检验需要两个条件：1）个体间相互独立；2）两组样本均来自正态分布总体；3）方差齐。我们可以利用直方图和QQ图观察样本是否呈正态分布，当然也可以直接用夏皮罗-威尔克（Shapiro-Wilk）假设检验的方法判断样本数据的分布。如果样本数据符合正态分布，那么可以接着用Bartlett检验不同组样本的方差齐性。两个条件都满足的话就可以做T检验了。我这里使用的是随机产生的数据，所以分析结果就不要太在意了。

group_A=data_ggplot$data_norm_shannon[data_ggplot$group=='A']

# 获取A组的多样性数据

group_B=data_ggplot$data_norm_shannon[data_ggplot$group=='B']

# 获取B组的多样性数据

hist(group_A)

hist(group_B)

# 绘制多样性指数分布直方图

# 观察形状若为倒钟形那便是接近正态分布的