########################     α-diversity
H<-vegan::diversity(bac_asv, "shannon")
S<-specnumber(bac_asv)
J<- H/log(S)
simp<-vegan::diversity(bac_asv, "simpson")


FT=as.factor(FT)
cycle=as.factor(cycle)
H_aov<-aov(H~straw*FT*cycle)
S_aov<-aov(S~straw*FT*cycle)
J_aov<-aov(J~straw*FT*cycle)
simp_aov<-aov(simp~straw*FT*cycle)
versatility_aov<-aov(env$versatility~straw*FT*cycle)
summary(H_aov)
summary(S_aov)
summary(J_aov)
summary(simp_aov)
summary(versatility_aov)
##random forest
library(rfPermute)
library(randomForest)
library(labeling)
library(rfUtilities)
library(ggplot2)
F.ran<-cbind(treatment,env,H,S,J,simp)
set.seed(123)
pp<-rfPermute(H~AK+AP+PH+LRNH4_N+LRN03_N+JTNH4_N+JTN03_N+N+C+NH4_loss+NO3_loss, data=env, ntree = 500, na.action = na.omit, nrep = 1000)

a=importance(pp, sort.by = NULL, decreasing = TRUE)[,1:2]
b=rp.importance(pp, sort.by = NULL, decreasing = TRUE)[,c(2,4)]
richness_dat2 <-merge( a,b,by = 'row.names', all = F)
row.names(richness_dat2) <- richness_dat2$Row.names#格式变了转回去
richness_dat2<-richness_dat2[,-1] 
library(tibble)
library(dplyr)
library(forcats)
library(ggplot2)
library(ggsci)
pdf("MSE_C_H.pdf",width=9.5,heigh=6.5)
richness_dat2 %>%
  as_tibble(rownames = "names") %>%
  data.frame() %>%
  mutate(label = if_else(X.IncMSE.pval < 0.001,"***",
                         if_else(X.IncMSE.pval <0.01,"**",
                                 if_else(X.IncMSE.pval<0.05,"*",""))),
         X.IncMSE = as.numeric(X.IncMSE)) %>%
  arrange(X.IncMSE) %>%
  mutate(group = if_else(label=="","In_sig","Sig"),
         names = forcats::fct_inorder(names)) %>%
  ggplot(aes(x = names, y = X.IncMSE))+
  geom_bar(aes(fill = label),stat = "identity")+
  scale_fill_npg()+
  theme_bw()+
  theme(legend.key.size = unit(35,"pt"))+
  theme(axis.text.x = element_text(size = 25,color="black"),axis.text.y = element_text(size = 20,color="black"))+
  theme(axis.title.x = element_text(size = 25),axis.title.y = element_text(size = 25))+
  theme(panel.grid.major = element_blank(),panel.grid.minor=element_blank(),legend.title =element_blank() )+
  theme(legend.position = 'none') +
  geom_text(aes(y = X.IncMSE+1 ,label = label),size=9)+
  labs(x = "Variables", y = "Increase in MSE_M (%)")+
  coord_flip()
dev.off()


################################# β-diversity
#Mantel test   

AK.dist<-distance(AK)
AP.dist<-distance(AP)
PH.dist<-distance(PH)
LRNH4_N.dist<-distance(LRNH4_N)
LRN03_N.dist<-distance(LRN03_N)
JTNH4_N.dist<-distance(JTNH4_N)
JTN03_N.dist<-distance(JTN03_N)
N.dist<-distance(N)
C.dist<-distance(C)
C.N.dist<-distance(C.N)
NH4_loss.dist<-distance(NH4_loss)
NO3_loss.dist<-distance(NO3_loss)
bac_asv.dist<-distance(bac_asv, "bray")


bac_asv_AK<-mantel(bac_asv.dist~AK.dist)
bac_asv_AP<-mantel(bac_asv.dist~AP.dist)
bac_asv_PH<-mantel(bac_asv.dist~PH.dist)

bac_asv_LRNH4_N<-mantel(bac_asv.dist~LRNH4_N.dist)
bac_asv_LRN03_N<-mantel(bac_asv.dist~LRN03_N.dist)
bac_asv_JTNH4_N<-mantel(bac_asv.dist~JTNH4_N.dist)
bac_asv_JTN03_N<-mantel(bac_asv.dist~JTN03_N.dist)
bac_asv_N<-mantel(bac_asv.dist~N.dist)
bac_asv_C<-mantel(bac_asv.dist~C.dist)
bac_asv_NH4_loss<-mantel(bac_asv.dist~NH4_loss.dist)
bac_asv_NO3_loss<-mantel(bac_asv.dist~NO3_loss.dist)


##PCoA
library(vegan)
library(ggplot2)
adonis(bac_asv ~ straw*FT*cycle,  perm=999)


#PCOA
distance <- vegdist(bac_asv, method = 'bray')
pcoa <- cmdscale(distance, k = (nrow(bac_asv) - 1), eig = TRUE)
ordiplot(scores(pcoa)[ ,c(1, 2)], type = 't')
pcoa$eig
point <- data.frame(pcoa$point)
pcoa_eig <- (pcoa$eig)[1:2] / sum(pcoa$eig)
sample_site <- data.frame({pcoa$point})[1:2]
names(sample_site)[1:2] <- c('PCoA1', 'PCoA2')
pcoa_eig
sample_site=read.csv("sample_site_C.csv")
sample_site$cycle<-factor(cycle)
sample_site$treat<-factor(treatment)
sample_site$FT<-factor(FT)
sample_site$straw<-factor(straw)
library(ggforce)
pdf("****.pdf", width=14, height=12)
ggplot(data = sample_site, aes(fill = factor(straw), x= PCoA1, y=PCoA2))+
  geom_point(aes(colour = factor(straw), shape=factor(cycle)), size=4.5)+
  scale_shape_manual(values = c(17, 15, 16, 18))+
  coord_cartesian(xlim = c(-0.4,0.3),ylim = c(-0.4,0.4))+
  ylab("PCo2 (1.82%)") + 
  xlab("PCo1 (8.15%)")+
  theme_bw()+
  theme(axis.text.x = element_text(size = 25,color="black"),axis.text.y = element_text(size = 25,color="black"))+
  theme(axis.title.x = element_text(size = 25),axis.title.y = element_text(size = 25))+
  theme(legend.key.size = unit(35,"pt"))+
  #geom_mark_ellipse(aes(fill=straw,colour=straw),alpha=0)
  #   geom_encircle(aes(fill=straw),expand=0,spread=0.5,s_shape=1,size=3,linetype = 1,alpha=0.2)+
  # stat_ellipse(data=sample_site,
  #             geom = "polygon",level=1,
  #             linetype = 1,size=0.8,
  #             aes(fill=straw),
  #             alpha=0.1,
  #             show.legend = F)
  #geom_mark_ellipse(data=sample_site,
  #aes(fill=straw,colour=straw),
  #alpha=0)+
scale_color_manual(values = c("#E69F00","#56B4E9", "#CC79A7","#F8766D","#00BFC4"))+
  guides(fill=F,color=F,shape=F) 

dev.off()
##NMDS
bac_asv.nmds<-metaMDS(bac_asv,trace = FALSE)
x <- bac_asv.nmds$points[,1]
y <- bac_asv.nmds$points[,2]

sample_site <- data.frame({bac_asv.nmds$point})[1:2]
aa=cbind(sample_site,all[,c(1:4)])

names(sample_site)[1:2] <- c('NMDS1', 'NMDS2')

sample_site$Cycle<-factor(Cycle)
sample_site$FT<-factor(FT)
sample_site$Straw<-factor(Straw)
pdf("nmds_CAZy.pdf", width=7, height=6)
ggplot(data = sample_site, aes(fill = factor(Straw), x= NMDS1, y=NMDS2))+
  geom_point(aes(colour = factor(Straw), shape=factor(Cycle)), size=2.5)+
  scale_color_manual(values = c("red", "orange", "brown", "green", "blue", "purple")) +
  scale_shape_manual(values = c(11, 15, 16, 17, 18))+
  ylab("NMDS2") + 
  xlab("NMDS1")+
  theme_bw()+
  theme(axis.text.x = element_text(size = 25,color="black"),axis.text.y = element_text(size = 25,color="black"))+
  theme(axis.title.x = element_text(size = 25),axis.title.y = element_text(size = 25))+
  theme(legend.key.size = unit(35,"pt"))
dev.off()



##################  network


bac1<-read.csv("allbac_asv_network.csv",head=T, row.names=1)
bac <- t(bac1[,-c(1:7)])

occor = corr.test( bac,use="pairwise",method="spearman",adjust="fdr",alpha=.05,ci=FALSE)
occor.r = occor$r 
occor.p = occor$p 

occor.r[occor.p>0.01|abs(occor.r)<0.75] = 0 

igraph = graph_from_adjacency_matrix(occor.r,mode="upper",weighted=TRUE,diag=FALSE)
igraph

bad.vs = V(igraph)[igraph::degree(igraph) == 0]
igraph = igraph::delete.vertices(igraph, bad.vs)
igraph

igraph.weight = E(igraph)$weight

E(igraph)$weight = NA

set.seed(123)
plot(igraph,main="allbac_asv network",vertex.frame.color=NA,vertex.label=NA,edge.width=1,
     vertex.size=5,edge.lty=1,edge.curved=TRUE,margin=c(0,0,0,0))


fc <- cluster_fast_greedy(igraph, weights = NULL)
modularity <- modularity(igraph, membership(fc))


comps <- membership(fc)


module_sizes2 <- table(comps)


module_order <- order(module_sizes2,decreasing = T)  
new_module_names <- c(1:length(module_sizes2)) 
comps <- factor(comps, levels = module_order, labels = new_module_names)  
module_sizes <- table(comps)




top_modules <- names(head(sort(module_sizes, decreasing = TRUE), N))
colbar = rainbow(nrow(as.data.frame(table(top_modules)))+1)


comps_o <- ifelse(comps %in% top_modules, comps, 1000)
V(igraph)$color <- ifelse(comps %in% top_modules, colbar[comps_o], "gray")



ggplot(data,aes(x=factor(treatment),y=abundance,fill=phylum))+
  geom_bar(stat="identity",position="fill",width=0.9)+
  #facet_wrap(~Cycle,nrow = 1)+###cycle合起来，cycle分开就不##
  scale_fill_manual(values=c("gray","red","orange" ,"yellow", "green","cyan", "blue", "purple","pink","brown","wheat","#4B0082"))+
  scale_x_discrete(limits=c("1CK", "2CK", "3CK","1RR","2RR","3RR"))+
  labs(y="Fun Relative Abundance",x="Treatment")+
  #coord_cartesian(ylim = c(0,100))+
  theme_bw()+
  theme(legend.key.size = unit(35,"pt"))+
  theme(axis.text.x = element_text(size = 25,color="black",angle = 45, hjust = 1, vjust = 1),axis.text.y = element_text(size = 25,color="black"))+
  theme(axis.title.x = element_text(size = 25),axis.title.y = element_text(size = 25))+
  theme(legend.position = "bottom")+
  guides(fill=F,color=F,shape=F)