4.精美的Alpha多样性箱线图

本文详细介绍了如何使用R语言中的相关库,如readxl、vegan等,处理OTU表和分组表,计算Alpha多样性指数,并通过ggplot2绘制带有显著性标记的箱线图,展示了如何展示和比较不同组别的多样性差异。
摘要由CSDN通过智能技术生成

看了网上的很多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

搞定!

  • 3
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值