导读
箱型图（Boxplot）或者盒图是一种能同时展示一组或多组数据的极值、四分位数、中位数和离群值，显示数据离散情况的统计图。下面介绍一下如何在R软件中利用箱形图可视化两组微生物群的Alpha多样性指数。过程分为以下几步：1）模拟原始数据；2）模拟分组；3）标准化；4）计算多样性；5）T检验；6）ggplot绘制Boxplot，保存结果。完成以上步骤将得到如下结果：
1.png
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=" ")
# 给样本命名，结果如下：
Species 1 Species 2 ... Species 10
Sample 1 36.70358 1255.7217
Sample 2 2082.09677 1022.0450
...
Sample 20
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)
# 做分组文件，如下：
sample_id group
1 Sample 1 A
...
10 Sample 10 A
11 Sample 11 B
...
20 Sample 20 B
3 标准化丰度
data_norm=data
# 新建data_norm数据框，赋初值
sample_sum=apply(data, 1, sum)
# 计算每个样本的菌种总丰度
for(i in 1:20)
{
for(j in 1:10){
data_norm[i,j]=data[i,j]/sample_sum[i]
# [每个样本的每个物种的丰度] / [该样本物种总丰度] = 相对丰度
}
}
apply(data_norm, 1, sum)
# 检查每个样本物种总丰度是否成功标准化到1，如下：
Sample 1 Sample 2 Sample 3 Sample 4 Sample 5 Sample 6 Sample 7
1 1 1 1 1 1 1
Sample 8 Sample 9 Sample 10 Sample 11 Sample 12 Sample 13 Sample 14
1 1 1 1 1 1 1
Sample 15 Sample 16 Sample 17 Sample 18 Sample 19 Sample 20
1 1 1 1 1 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_norm_shannon
Sample 1 2.209006
...
Sample 20 2.198923
data_ggplot=data.frame(data_ggplot, data_group["group"])
# 添加分组信息，如下：
data_norm_shannon group
Sample 1 2.209006 A
...
Sample 10 2.241821 A
Sample 11 2.100319 B
...
Sample 20 2.198923 B
5 T 检验
成组设计T检验需要两个条件：1）个体间相互独立；2）两组样本均来自正态分布总体；3）方差齐。我们可以利用直方图和QQ图观察样本是否呈正态分布，当然也可以直接用夏皮罗-威尔克（Shapiro-Wilk）假设检验的方法判断样本数据的分布。如果样本数据符合正态分布，那么可以接着用Bartlett检验不同组样本的方差齐性。两个条件都满足的话就可以做T检验了。我这里使用的是随机产生的数据，所以分析结果就不要太在意了。
5.1 直方图
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)
# 绘制多样性指数分布直方图
# 观察形状若为倒钟形那便是接近正态分布的
2.png
5.2 QQ图
qqnorm(group_A)
qqnorm(group_B)
# 绘制多样性指数QQ图
# 观察形状是一条连接主对角线的线那便是接近正态分布
3.png
5.3 Shapiro-Wilk正态性检验
shapiro.test(data_ggplot$data_norm_shannon)
# 夏皮罗-威尔克（Shapiro-Wilk）检验正态性，p>0.05，接受原假设，符合正态分布
4.png
5.4 Bartlett方差齐性检验
bartlett.test(data_norm_shannon~group,data=data_ggplot)
# 巴特利特（Bartlett）检验方差齐性，p>0.05，接受原假设，即两样本数据方差齐
5.png
5.5 T 检验
with(data_ggplot,t.test(formula=data_norm_shannon~group,conf.level=0.95))
# T检验，p>0.05，没法拒绝原假设，两组Shannon多样性指数差异不显著
6.png
6 ggplot绘制箱形图
library(ggplot2)
# 加载R包ggplot2
alpha_boxplot=ggplot(data_ggplot, aes(x=group, y=data_norm_shannon, fill=group))+
# 添加数据、xy值、 颜色参数给画图函数ggplot
geom_boxplot()+
# 盒图
labs(title="Alpha diversity", x="Group", y="Shannon index")+
# 标题
theme(plot.title=element_text(hjust=0.5), legend.title=element_blank())
# 标题居中
pdf('result.pdf')
alpha_boxplot
dev.off()
# 保存结果，打开result.pdf文件，结果如下：
1.png
参考
1）R语言-t检验
https://blog.csdn.net/weixin_38322363/article/details/89811964
2）Shapiro-Wilk检验
https://www.jianshu.com/p/e202069489a6
3）方差齐性检验的原理
https://zhuanlan.zhihu.com/p/26943946
同步发布于微信公众号：微生态
