# dune是一套包含了20个样品和30个物种丰度数据的统计表。
# 其格式是常见OTU表转置后的格式，
# 每一行是一个样品，每一列是一个物种 (检测指标)。
# dune.env是元数据信息，包含数据的分子信息、生存环境信息等，记录了5个因素 (同时包含连续变量信息和分组变量信息)：
# A1: 土壤厚度信息 a numeric vector of thickness of soil A1 horizon.
# Moisture: 湿度等级信息，分4个等级，1 < 2 < 4 < 5.
# Management: 分组信息，不同的管理方式 a factor with levels: BF (Biological farming), HF (Hobby farming), NM (Nature Conservation Management), and SF (Standard Farming).
# Use: 一个分组信息 an ordered factor of land-use with levels: Hayfield < Haypastu < Pasture.
# Manure: 一个分组信息，0 < 1 < 2 < 3 < 4
library(vegan)
data(dune)
data(dune.env)
## [1] 20 30
head(dune)
head(dune.env)
#绘制一个PcOA的图看一下
# 计算加权bray-curtis距离
dune_dist <- vegdist(dune, method="bray", binary=F)
dune_pcoa <- cmdscale(dune_dist, k=3, eig=T)
dune_pcoa_points <- as.data.frame(dune_pcoa$points)
sum_eig <- sum(dune_pcoa$eig)
eig_percent <- round(dune_pcoa$eig/sum_eig*100,1)
colnames(dune_pcoa_points) <- paste0("PCoA", 1:3)
dune_pcoa_result <- cbind(dune_pcoa_points, dune.env)
head(dune_pcoa_result)
library(ggplot2)
ggplot(dune_pcoa_result, aes(x=PCoA1, y=PCoA2, color=Management)) +
  labs(x=paste("PCoA 1 (", eig_percent[1], "%)", sep=""),
       y=paste("PCoA 2 (", eig_percent[2], "%)", sep="")) +
  geom_point(size=4
  ) + stat_ellipse(level=0.6) +
  theme_classic()
#样品中重复太少了，做不出置信椭圆。换个方式，用ggalt包中的geom_encircle把样品包起来。
install.packages("ggalt")
library(ggalt)
ggplot(dune_pcoa_result, aes(x=PCoA1, y=PCoA2, color=Management, group = Management)) +
  labs(x=paste("PCoA 1 (", eig_percent[1], "%)", sep=""),
       y=paste("PCoA 2 (", eig_percent[2], "%)", sep="")) +
  geom_point(size=5) + 
  geom_encircle(aes(fill=Management), alpha = 0.1, show.legend = F) +
  theme_classic() + coord_fixed(1)
# 基于bray-curtis距离进行计算
set.seed(1)
dune.div <- adonis2(dune ~ Management, data = dune.env, permutations = 999, method="bray")
dune.div
dune_adonis <- paste0("adonis R2: ",round(dune.div$R2,2), "; P-value: ", dune.div$`Pr(>F)`)
library(ggalt)
p<-ggplot(dune_pcoa_result, aes(x=PCoA1, y=PCoA2, color=Management, group = Management)) +
  labs(x=paste("PCoA 1 (", eig_percent[1], "%)", sep=""),
       y=paste("PCoA 2 (", eig_percent[2], "%)", sep=""),
       title=dune_adonis) +
  geom_point(size=5) + 
  geom_encircle(aes(fill=Management), alpha = 0.1, show.legend = F) +
  theme_classic() + coord_fixed(1)
#整体有差异了，后面就看看具体那两组之间有差异，哪两组之间无差异~~~
#这时就需要pairwise.adonis来进行配对检验了。
head(dune)
head(dune.env)
library(pairwiseAdonis)
dune.pairwise.adonis <- pairwise.adonis(x=dune, factors=dune.env$Management, sim.function = "vegdist",
                                        sim.method = "bray",
                                        p.adjust.m = "BH",
                                        reduce = NULL,
                                        perm = 999)
dune.pairwise.adonis
library(ggpubr)
library(patchwork)
tab2 <- ggtexttable(dune.pairwise.adonis[,c("pairs","R2","p.value","p.adjusted")], rows = NULL, 
                    theme = ttheme("blank")) %>% 
  tab_add_hline(at.row = 1:2, row.side = "top", linewidth = 1)  %>% 
  tab_add_hline(at.row = nrow(dune.pairwise.adonis)+1, row.side = "bottom", linewidth = 1)  

p2 = p + tab2 
p2
# 调布局
p2 + plot_layout(design=c(area(1,1), area(2,1)))
# p / tab2