install.packages("tidyselect")
install.packages("rstatix")
install.packages("ggpubr")
install.packages("rio")
install.packages("zip")
library(ggplot2)
library('ggpubr')
alpha <- read.table(file.choose(),header = TRUE,row.names = 1,sep = "\t") #各样本α多样性指数
group <- read.table(file.choose(),header = TRUE,sep = "\t")
rownames(group) = group[,1]
group <- group[rownames(alpha),]
group$group <- factor(group$group,labels = c('0m','25m','50m','117m','150m', '200m'))
group$group <- factor(group$group,labels = c('0.2-0.6μm','0.6-1.2μm','1.2-2μm','2-20μm','>20μm'))
data <- cbind(group,alpha)
##### 作图
library(reshape)
View(data)
windowsFonts(TNM = windowsFont("Times New Roman"))#设置字体
#install.packages("showtext")
#library(showtext)
#install.packages("grDevices")
#library(grDevices)
#install.packages("Cairo")
#library(Cairo)输出图片用的,还不太会用
plotdata = melt(data,id.vars = 1:2,measure.vars = 3:10,variable_name = "alpha",value.name = "value",family="TNM")
#CairoTIFF("p1.tiff", bg = "transparent",width = 800,height = 600)#,compression="lzw"
p = ggplot(data=plotdata, aes(x=group,y=value,family="TNM"))+geom_boxplot(aes(fill=alpha),width=0.3,)+
theme(axis.text.x =element_text(angle = 45,hjust = 0.5, vjust = 0.5, size= 8,family="TNM",face = "bold",color = "black"),axis.text.y =element_text(size = 7,color = "black",family="TNM",face = "bold"),legend.text = element_text(family="TNM"),
panel.grid.minor = element_blank() )
# 根据多样性指数进行分面
p = p + facet_wrap(~ alpha, scales="free") +
theme(strip.text = element_text(family="TNM"))+
scale_y_continuous(labels = scales::label_comma(accuracy =0.01))
p1 = p +stat_compare_means(method = "anova",aes(label = paste0(..p.signif..)), label.x = c(),hide.ns = TRUE, size=5,family="TNM",label.y.npc = 0.85)+#只显示p值#+label.x = 3,
labs(x = NULL,y="Species abundance (log to base 10)") +
#coord_flip()+#横纵坐标转置
#stat_compare_means(aes(label = ..p.signif..),method = "t.test",ref.group = ".all.", hide.ns = TRUE,)+
theme(panel.border = element_rect(fill = NA) ,panel.grid = element_blank(),panel.background = element_blank(),
axis.title.y = element_text(size = 13,family="TNM", face = "bold"))+
guides(fill = F)#panel.grid绘图区网格线
ggsave("p1.tiff",dpi=600)#输出图片