library(vegan)
library(ggplot2)
library(RColorBrewer)
library(cowplot)
############first step: import data##############################################
#OTU 
#otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- read.csv('per-i1j3.csv',row.names = 1)
otu <- data.frame(t(otu))
#samples group # 
group <- data.frame(rownames(otu))
group$site <- c(rep('A', 36), rep('B', 30), rep('C', 36), rep('D', 27))
group$site
#####################second step，PERMANOVA analysis###########################
# Bray-Curtis distance measure
adonis_result <- adonis2(otu~site, group, distance = 'bray', permutations = 999)
adonis_result
#output 
otuput <- data.frame(adonis_result, check.names = FALSE, stringsAsFactors = FALSE)
otuput <- cbind(rownames(otuput), otuput)
names(otuput) <- c('', 'Df', 'Sums of squares', 'Variation (R2)','F.Model',  'Pr (>F)')
write.table(otuput, file = 'PERMANOVA.result_all.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')
# get R2 and Pvalue
#otu_adonis_R2 <- paste("R^2: ",round(adonis_result$aov.tab$R2[1],2))
#otu_adonis_P <- paste("P-value: ", round(adonis_result$aov.tab$`Pr(>F)`[1],2))
#PCoA plot
bray_dis <- vegdist(otu, method = 'bray')    
pcoa <- cmdscale(bray_dis, k = 2, eig = TRUE)
pcoa_eig <- pcoa$eig
pcoa_exp <- pcoa$eig/sum(pcoa$eig)
site <- scores(pcoa)

pcoa1 <- paste('PC1（', round(100*pcoa_exp[1], 2), '%）')
pcoa2 <- paste('PC2（', round(100*pcoa_exp[2], 2), '%）')
site <- data.frame(pcoa$point)[1:2]
site$name <- rownames(site)
site$Site <- c(rep('A', 36), rep('B', 30), rep('C', 36), rep('D', 27))
p4 <- ggplot(data = site, aes(X1, X2)) +geom_point(aes(shape=Site,color=Site))+
  labs(x = pcoa1, y = pcoa2, title = 'PCoA (i=1，j=3)')+
  stat_ellipse(aes(fill = Site), geom = 'polygon', level = 0.95, alpha = 0.5, show.legend = T) +   #，it is not clustering
  scale_color_manual(values = c("#404D5B","#3CA48D","#5492C7","#99B86B"))+
  scale_fill_manual(values = c("#404D5B","#3CA48D","#5492C7","#99B86B")) +
theme(axis.text = element_text(size=10),
        #panel.border = element_rect(color = "black"),
        #axis.line = element_line(linewidth = 0.5),
        plot.title=element_text(hjust=0,size=12),
        legend.title = element_text(size=16),
        legend.text = element_text(size=12),
        legend.position = "right") +
  geom_vline(xintercept = 0, color = 'gray', size = .1) +
  geom_hline(yintercept = 0, color = 'gray', size = .1) +
  labs(x = pcoa1, y = pcoa2, size=3)+
  annotate('text', label = expression(R^2:0.55), x = 0.5, y = 0.5, size = 5, colour = 'black')+
  annotate('text', label = expression(italic(P):  0.001), x = 0.5, y = 0.45, size = 5, colour = 'black')+
  ylim(-0.4,0.52)
p4
#tiff(file=paste("fig.4.minor",".tiff",sep = ""),
     #width =8,height=7.5,units = "in", res=400)
#p4
#dev.off()
#getwd()
####################third step，The differences among the small groups in the whole ##################################
#devtools::install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
#install.packages("pairwiseAdonis")
library(pairwiseAdonis)
otu.pairwise.adonis <- pairwise.adonis(x=otu, factors=group$site, sim.function = "vegdist",
                                       sim.method = "bray",
                                       p.adjust.m = "BH",
                                       reduce = NULL,
                                       perm = 999)
otu.pairwise.adonis
otu.pairwise.adonis <- data.frame(otu.pairwise.adonis, stringsAsFactors = FALSE)

###output
write.table(otu.pairwise.adonis, 'otu.pairwise.adonis.txt', row.names = FALSE, sep = '\t', quote = FALSE, na = '')
##
#install.packages("ggpubr")
#install.packages("pathchwork")
library(ggpubr)
library(patchwork)
table <- ggtexttable(otu.pairwise.adonis[,c("pairs","p.adjusted","sig")], rows = NULL, 
                     theme = ttheme("blank")) %>% 
  tab_add_hline(at.row = 1:2, row.side = "top", linewidth = 1)  %>% 
  tab_add_hline(at.row = nrow(otu.pairwise.adonis)+1, row.side = "bottom", linewidth = 1)  
table 
qqq=ggdraw(p4)+draw_plot(table,x=0.25,y=-0.43,scale = 0.4,width = 1,height = 1.3)
qqq
tiff(file=paste("per-i1j3",".tiff",sep = ""),
     width =8,height=7.5,units = "in", res=300)
qqq
dev.off()