看了网上的很多A多样性箱线图,但是都很一般,而且带显著性标记的也很少,那么今天就出一期如何做出好看的A多样性箱线图。
如何从上图画出下图的感觉!
先说下需要的两个文件:抽平过的OTU表和分组表
OTU表:flat_ASV.xlsx
分组表:group.xlsx
话不多说,上代码!
library(readxl);library(vegan)
ASV <- read_xlsx("flat_ASV.xlsx",sheet="flat_ASV")
ASV=as.data.frame(ASV)
# 将ASV矩阵的第一列设置为行名
row.names(ASV) <- ASV[, 1]
# 从ASV矩阵中移除第一列
ASV <- ASV[, -1]
# 从ASV矩阵中移除第1和第2列
ASV <- ASV[, -c(1, 2)]
ASV <- t(ASV)
#构建Alpha多样性指数计算函数#####
alpha <- function(x,tree = NULL,base = exp(1)) {
est <- estimateR(x)
Richness <- est[1, ]
Chao1 <- est[2, ]
ACE <- est[4, ]
Shannon <- diversity(x, index = "shannon", base = base)
Simpson <- diversity(x, index = "simpson") #Gini-Simpson 指数
Pielou <- Shannon / log(Richness, base)
goods_coverage <- 1 - rowSums(x == 1) / rowSums(x)
result <- data.frame(Richness, Chao1, ACE,Shannon, Simpson, Pielou, goods_coverage)
if (!is.null(tree)) {
PD_whole_tree <- pd(x, tree, include.root = FALSE)[1]
names(PD_whole_tree) <- 'PD_whole_tree'
result <- cbind(result, PD_whole_tree)
}
result
}
alpha_index <- alpha(ASV)
#write.csv(alpha_index,file = "alpha_index.csv", quote = FALSE)
group <- read_xlsx("group.xlsx")
group=as.data.frame(group)
# 将ASV矩阵的第一列设置为行名
row.names(group) <-group[, 1]
group[, "sample"] <- NULL # 替换'ColumnName'为实际列名
# 按行名合并alpha_index和group矩阵
dat <- merge(alpha_index, group, by = "row.names", all = TRUE)
library(reshape2)
dat <- dat[,-c(2,4,8)]
dat <- melt(dat,id.vars = c(1,6),variable.name = "Alpha")#短数据变长数据
head(dat, n = 3)
rk <- c("Model","TMZ")
dat$group <- factor(dat$group,levels = rk) #设定作图顺序
names(dat)[4] <- 'V' #重新命名value列名
names(dat)[2] <- 'Group'
#计算均值标准差#####
Mean <- aggregate(dat$V,list(dat$Group,dat$Alpha),mean)
Sd <- aggregate(dat$V,list(dat$Group,dat$Alpha),sd)
Max <- aggregate(dat$V,list(dat$Group,dat$Alpha),max)
colnames(Mean) <- c("Group","Alpha","mean")
colnames(Sd) <- c("Group","Alpha","sd")
colnames(Max) <- c("Group","Alpha","max")
df3 <- as.data.frame(Reduce(function(x,y) merge(x,y,all = TRUE),list(Mean,Sd,Max)))
head(df3)
#write.csv(df3,file = "alpha_meansd.csv",row.names = FALSE) #保存文件
#LSD方差分析#####
ONE_LSD <- function(data,group,compare,value){
library(agricolae)
a <- data.frame(stringsAsFactors = F)
type <- unique(data[,group])
for (i in type)
{
# sub_dat <- subset(data,group == i)
sub_dat <- data[data[,group]==i,]
# fit <- aov(value ~ compare,sub_dat)
fit <- aov(sub_dat[,value] ~ sub_dat[,compare] )
out <- LSD.test(fit,'sub_dat[, compare]',p.adj='BH')
out$groups$type <- i
out$groups$compare <- rownames(out$groups)
a <- rbind(a,merge(out$means[,1:2], out$groups,by='sub_dat[, value]'))
}
names(a) <- c('mean','std','lsd',group,compare)
return(a)
}
df4 <- ONE_LSD(data=dat,group='Alpha',compare='Group',value='V')
head(df4)
lab <- rep(c("a","b"),each =4)#根据def4lsd的结果手动输入差异符号
data <- data.frame(df3,lab = lab)
#作图#####
library(ggplot2);library(ggsci)
col<-c( "#E81D22","#20e81d")
#箱线图#####
#windowsFonts(myFont=windowsFont("Arial"))
data$Alpha <- factor(data$Alpha,levels = c("Chao1","Shannon","Simpson","Pielou"))
P1 <- ggplot(data = dat,aes(x=Group,y=V))+
stat_boxplot(linewidth=1,width=0.4,geom = "errorbar",alpha=0.9,aes(color=Group))+ ##由于自带的箱形图没有胡须末端没有短横线,使用误差条的方式补上
geom_boxplot(linewidth=0.8,aes(fill=Group,color=Group),alpha=0.2,outlier.fill="white",outlier.shape=7,outlier.color="red")+ #outlier设置异常值各种参数
geom_jitter(aes(fill=Group),colour="white",shape = 21,size=2,width=0.3,alpha=0.6)+ #在箱线图上添加点图(抖点图),更好反应数据分布于
scale_fill_manual(values = col)+
scale_color_manual(values = col)+
facet_wrap(.~Alpha,scales = "free_y")+
labs(x=NULL,y="Alpha Diversity",fill = NULL,color = NULL)+ #这里如果不加color = NULL,就会出现两个图例
scale_x_discrete(breaks=dat$Group,labels=dat$Group)+
geom_text(data = data,aes(x=Group,y=max,label=lab),family = "serif",vjust = -0.5,show.legend = FALSE)+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 90))+
theme_bw()+
theme(legend.background=element_rect(color="white"),panel.grid.major = element_blank(), #去除网格线
panel.grid.minor = element_blank(),strip.text.x = element_text(size = 16, color = "black") ,strip.background = element_rect(
color="black", fill="white", linewidth=0.5, linetype="solid"),panel.border = element_rect(color = "black"),axis.title = element_text(size = 15),
axis.text = element_text(color = "black",size = 12),
legend.text = element_text(size = 12,color = "black"),
text = element_text(family = "serif"))
P1
搞定!