# rm(list=ls())
# # If devtools package is not installed, first install it
# #install.packages("devtools")
# #devtools::install_github("ChiLiubio/microeco")
# pacman::p_load(microeco,tidytree,treeio,ggplot2,ggtree,microeco)
# # Let's create a microtable object with more information
# dataset <- microtable$new(sample_table = sample_info_16S,
#                           otu_table = otu_table_16S,
#                           tax_table = taxonomy_table_16S,
#                           phylo_tree = phylo_tree_16S)
# sample_table<-dataset$sample_table
# otu_table<-dataset$otu_table
# tax_table<-dataset$tax_table
# phylo_tree<-dataset$phylo_tree
# write.csv(otu_table,"feature_table.csv")
# write.csv(sample_table,"sample_table.csv")
# write.csv(tax_table,"tax_table.csv")
# write.tree(as.phylo(phylo_tree),"phylo_tree_16S.nwk")
#参考链接：https://mp.weixin.qq.com/s/xbuycJN52fZhP8zNG1eJ9A
##Step1：创建microtable对象
rm(list=ls())
pacman::p_load(microeco,tidytree,treeio,ggplot2,ggtree,microeco,aplot)
feature_table <- read.csv('feature_table.csv', row.names = 1)
sample_table <- read.csv('sample_table.csv', row.names = 1)
tax_table <- read.csv('tax_table.csv', row.names = 1)
# 创建microtable对象
dataset <- microtable$new(sample_table = sample_table,
                          otu_table = feature_table, 
                          tax_table = tax_table)
dataset
head(sample_table)
head(feature_table)
head(tax_table)

##Step2：PCoA分析
# PCoA
dataset$cal_betadiv()
t1 <- trans_beta$new(dataset = dataset, group = "Group", measure = "bray")
t1$cal_ordination(ordination = "PCoA")
tmp <- t1$res_ordination$scores
# differential test with trans_env class
t2 <- trans_env$new(dataset = dataset, add_data = tmp[, 1:2])
# 'KW_dunn' for non-parametric test
t2$cal_diff(group = "Group", method = "anova")
t2
#主图
p1 <- t1$plot_ordination(plot_color = "Group",
                         plot_shape = "Group",
                         plot_type = c("point", "ellipse"))
###PCo1 箱型图1  
p2 <- 
  t2$plot_diff(measure = "PCo1", add_sig = T) + 
  theme_bw() + 
  coord_flip() + 
  theme(legend.position = "none", 
        axis.title.x = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank())
###PCo2 箱型图2
p3 <- t2$plot_diff(measure = "PCo2", add_sig = T) + 
  theme_bw() + 
  theme(legend.position = "none", 
        axis.title.y = element_blank(), 
        axis.text.x = element_blank(), 
        axis.ticks.x = element_blank())
### 组合图
g <- p1 %>% 
  insert_top(p2, height = 0.2) %>% 
  insert_right(p3, width = 0.2)
g

#上图的水平轴和右侧图的垂直轴与主图的轴线并不完全对应
####手动控制坐标范围
#对PCOA的主成分坐标缩放，统一坐标
x_lim <- range(tmp[, 1]) * 1.4
y_lim <- range(tmp[, 2]) * 1.4
p1 <- p1 + scale_y_continuous(limits = y_lim, expand = c(0, 0)) + 
  scale_x_continuous(limits = x_lim, expand = c(0, 0))
p2 <- p2 + scale_y_continuous(limits = x_lim, expand = c(0, 0))
p3 <- p3 + scale_y_continuous(limits = y_lim, expand = c(0, 0))

g1 <- p1 %>% 
  insert_top(p2, height = 0.2) %>% 
  insert_right(p3, width = 0.2)
g1

#切换主题
p2 <- p2 + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
p3 <- p3 + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
g2 <- (p1 + theme_classic()) %>% 
  insert_top(p2, height = 0.15) %>% 
  insert_right(p3, width = 0.15)
g2