######Fig 1, 2, 5######


library(tidyverse)
library(agricolae)
library(car)
library(reshape2)
library("ggplot2")
library("ggpubr")
data_wt = read.table("enzyme_data.txt", header=T, sep="\t")
head(data_wt)
data_wt = dcast(data_wt, ID + group ~ grou, value.var = "count")
head(data_wt)
plot_type = "box"
for (i in 3:ncol(data_wt))
  # i = 3
  ss = data_wt[i]
  colnames(ss) = c("count")
  ss$group = data_wt$group
  # ??????????????? Shapiro-Wilk normality test,??????p-value
  xx =  shapiro.test(ss$count)
  p1 = xx[[2]]
  # ?????????????????? Bartlett test of homogeneity of variances,??????p-value
  xc = bartlett.test(count ~ group, data = ss)
  p2 = xc[[3]]
  # ??????2. ???????????????
  if ( plot_type == "box")
    
    p1 = round(p1,3)
    p2 = round(p2,3)
    data_i = data_wt[i]
    ee	=as.matrix(data_i)
    dd = as.vector(ee)
    name_i = colnames(data_wt[i])
    model=aov(dd ~ group, data=data_wt)#????????????
    wtx1 = summary(model)
    wtx2 = wtx1[[1]]
    wtx3 = wtx2[5]
    out = LSD.test(model,"group", p.adj="none")#??????????????????,?????????P???
    aa = out$group#????????????:???????????????
    aa$group = row.names(aa)
    a = max(aa$dd)*1.2
    
    data_box = data_wt[c(1,2,i)]
    colnames(data_box) = c("ID" , "group","dd" )
    
    # out = LSD.test(model,"group", p.adj="none") # alternative fdr
    stat = out$groups
    data_box$stat=stat[as.character(data_box$group),]$groups
    max=max(data_box[,c("dd")])
    min=min(data_box[,c("dd")])
    x = data_box[,c("group","dd")]
    y = x %>% group_by(group) %>% summarise_(Max=paste('max(',"dd",')',sep=""))
    y=as.data.frame(y)
    rownames(y)=y$group
    data_box$y=y[as.character(data_box$group),]$Max + (max-min)*0.05
    
    
    
    # mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A")
    p = ggplot(data_box, aes(x=group, y=data_box[["dd"]], color=group)) +
      geom_boxplot(alpha=1, outlier.size=2, size=1.4, width=0.6, fill="transparent") +
      labs(x="Treatments",
           y=paste(name_i, sep = "_"),
           title = paste("",sep = ":"))+
      geom_text(data=data_box, aes(x=group, y=y, color=group, label= stat)) +
      geom_jitter(position=position_jitter(0.20),size=1.8, alpha=0.55)+theme(legend.position="none")
    p
    p=p+theme_bw()+
      # scale_fill_manual(values = mi, guide = guide_legend(title = NULL))+
      theme
        
        panel.grid.major=element_blank()
        panel.grid.minor=element_blank()
      
        plot.title = element_text(vjust = -8.5,hjust = 0.1)
        axis.title.y =element_text(size = 12,face = "bold",colour = "black")
        axis.title.x =element_text(size = 12,face = "bold",colour = "black")
        axis.text = element_text(size = 12,face = "bold")
        axis.text.x = element_text(colour = "black",size = 10)
        axis.text.y = element_text(colour = "black",size = 10)
        legend.text = element_text(size = 12,face = "bold")
        legend.position = "none"#??????????????????
        
  
    if (length(unique(data_box$group))>3){	p=p+theme(axis.text.x=element_text(angle=45,vjust=1, hjust=1))}
    FileName = paste(name_i,"_aov_LSD_box", ".pdf", sep = "")
    ggsave(FileName, p, width = 6, height = 6)
    
  
    ######Fig 7A#####
    
    
    library(vegan)
    library(ggplot2)
    library(ggrepel)
    fc=read.table("aha_soil.txt",header = T,row.names = 1)#????????????????????????
    sp=read.table("aha_bac_p.txt",header = T,row.names = 1)#????????????????????????
    spp=decostand(sp,method = "hellinger")#????????????????????????
    fcc=log10(fc)#????????????????????????
    uu=rda(spp~.,fcc)#RDA??????
    ii=summary(uu)  #??????????????????
    sp=as.data.frame(ii$species[,1:2])*2#?????????????????????,??????????????????????????????????????????,??????
    st=as.data.frame(ii$sites[,1:2])
    yz=as.data.frame(ii$biplot[,1:2])
    grp=as.data.frame(c(rep("CK",3),rep("FM1",3),rep("FM2",3),rep("FM2",3)))#????????????????????????
    colnames(grp)="group"
    ggplot() +
      #geom_text_repel(data = st,aes(RDA1,RDA2,label=row.names(st)),size=4)+#????????????
      geom_point(data = st,aes(RDA1,RDA2,shape=grp$group,fill=grp$group),size=4)+
      scale_shape_manual(values = c(21:25))+
      geom_segment(data = sp,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                                 type = "closed"),linetype=2, size=0,colour = "white")+
      geom_text_repel(data = sp,aes(RDA1,RDA2,label=row.names(sp)))+
      geom_segment(data = yz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                                 type = "closed"),linetype=1, size=1,colour = "black")+
      geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+
      labs(x=paste("RDA 1 (", format(100 *ii$cont[[1]][2,1], digits=4), "%)", sep=""),
           y=paste("RDA 2 (", format(100 *ii$cont[[1]][2,2], digits=4), "%)", sep=""))+
      geom_hline(yintercept=0,linetype=3,size=1) + 
      geom_vline(xintercept=0,linetype=3,size=1)+
      guides(shape=guide_legend(title=NULL,color="black"),
             fill=guide_legend(title=NULL))+
      theme_bw()+theme(panel.grid=element_blank())
    
    
    #####Fig 7B######
    
    library(vegan)
    library(ggplot2)
    library(ggrepel)
    fc=read.table("aha_soil.txt",header = T,row.names = 1)#????????????????????????
    sp=read.table("aha_fungi_p.txt",header = T,row.names = 1)#????????????????????????
    spp=decostand(sp,method = "hellinger")#????????????????????????
    fcc=log10(fc)#????????????????????????
    uu=rda(spp~.,fcc)#RDA??????
    ii=summary(uu)  #??????????????????
    sp=as.data.frame(ii$species[,1:2])*2#?????????????????????,??????????????????????????????????????????,??????
    st=as.data.frame(ii$sites[,1:2])
    yz=as.data.frame(ii$biplot[,1:2])
    grp=as.data.frame(c(rep("CK",3),rep("FM1",3),rep("FM2",3),rep("FM2",3)))#????????????????????????
    colnames(grp)="group"
    ggplot() +
      #geom_text_repel(data = st,aes(RDA1,RDA2,label=row.names(st)),size=4)+#????????????
      geom_point(data = st,aes(RDA1,RDA2,shape=grp$group,fill=grp$group),size=4)+
      scale_shape_manual(values = c(21:25))+
      geom_segment(data = sp,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                                 type = "closed"),linetype=2, size=0,colour = "white")+
      geom_text_repel(data = sp,aes(RDA1,RDA2,label=row.names(sp)))+
      geom_segment(data = yz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                                 type = "closed"),linetype=1, size=1,colour = "black")+
      geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+
      labs(x=paste("RDA 1 (", format(100 *ii$cont[[1]][2,1], digits=4), "%)", sep=""),
           y=paste("RDA 2 (", format(100 *ii$cont[[1]][2,2], digits=4), "%)", sep=""))+
      geom_hline(yintercept=0,linetype=3,size=1) + 
      geom_vline(xintercept=0,linetype=3,size=1)+
      guides(shape=guide_legend(title=NULL,color="black"),
             fill=guide_legend(title=NULL))+
      theme_bw()+theme(panel.grid=element_blank())
    
    
    #####Fig 6#####
    
    library(vegan)
    library(ape)
    library(ggplot2)
    library(ggrepel)
    setwd("G:\\Baiduyun\\CAFM\\picture")
    #?????????????????????????????????????????????????????????????????????
    data <- read.csv("NP_otu.txt", head=TRUE,sep="\t",row.names = 1)
    groups <- read.table("NP_group.txt",sep = "\t",header = F,colClasses = c("character"))
    groups <- as.list(groups)
    data <- t(data)
    data[is.na(data)] <- 0
    data <- vegdist(data,method = "bray")
    pcoa<- pcoa(data, correction = "none", rn = NULL)
    PC1 = pcoa$vectors[,1]
    PC2 = pcoa$vectors[,2]
    plotdata <- data.frame(rownames(pcoa$vectors),PC1,PC2,groups$V2)
    colnames(plotdata) <-c("sample","PC1","PC2","group")
    pich=c(21:24)
    cbbPalette <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442")
    Palette <- c("#000000", "#000000", "#000000", "#000000")
    pc1 <-floor(pcoa$values$Relative_eig[1]*100)
    pc2 <-floor(pcoa$values$Relative_eig[2]*100)
    
    otu.adonis=adonis(data~V2,data = groups,distance = "bray")
    pdf("NP_PCOA.pdf", height=6,width=6)
    
    
    #####Fig 3, 4#####
    
    
    setwd("C:\\Users\\Administrator\\Desktop\\Pediologia")
    library(ggplot2)
    data_06<-read.csv("Fild_CK_taxon.CSV",header = TRUE)
    p=ggplot(data_06,aes(x=taxon,y=value,fill=variable))
    p=p+geom_bar(stat="identity")+
      guides(fill=guide_legend(title = "Microbial phylum "))+
      xlab("")+ylab("Relative abundance (%)")+theme_bw()+
      scale_y_continuous(breaks = seq(0,100,25))+labs(title = "")
    p=p+theme(plot.title=element_text(hjust = 0.5,face = "bold"),
              axis.title=element_text(face = "bold",size = 12),
              legend.title=element_text(face = "bold",size = 10),
              axis.text.x=element_text(face = "bold",size = 10,angle=60,hjust=1))