# 引入包
install_if_missing <- function(pkg , BiocM=FALSE ) {
  if (!require(pkg, character.only = TRUE)) {
    if( BiocM ){
      BiocManager::install(pkg ,update=FALSE )
    }else{
      install.packages(pkg, dependencies = TRUE)
    }
    suppressPackageStartupMessages( library(pkg, character.only = TRUE) )
  }
}
required_packages <- c("vegan", "ggpubr")
invisible(lapply(required_packages, install_if_missing))

sangon <- c('#CA0000','#0087FF','#00BA1D','#CF00FF','#00DBE2','#FFAF00','#0017F4','#006012','#6F1A1A','#B5CF00','#009E73','#F0E442','#D55E00','#CC79A7','#FF8A8A','#AA6400','#50008A','#00FF58','#E175FF','#FFCC99','#009999','#CC0066','#C0C0C0','#666600','#505050','#56B4E9','#0072B2','#33FF33','#99004c','#CCFF99','#660066','#9370DB','#D8BFD8','#BC8F8F','#2F4F4F','#FF6347','#877878','#CD5C5C','#FF0000','#00FF00','#000080','#0D7339','#3D90D9','#D98236','#DF514F','#8CCDA3','#F5E2D7','#FACBBE','#EA4885','#E69F00')

name <- 'genus_bray' 
data <- read.table('genus_bray_distance.xls', header=T, sep='\t', check.names=F, row.names=1)
rownames(data) <- colnames(data) 
sample_list <- read.table('Sample_list.txt', sep='\t', header=F, check.names=F, colClasses='character', na.strings='')
data <- data[match(sample_list[, 2], rownames(data)), match(sample_list[, 2], colnames(data)), drop=F]
if (nrow(data) < 3){
  stop('have only less than three sample. will not do NMDS.\n', call.=FALSE)
}
k <- ifelse(nrow(data) == 3, 2, 3)
sol <- metaMDS(as.dist(data), try=20, trymax=1000, k=k)
write.table(sol$points, paste0(name, '_nmds_axis.xls'), sep='\t', col.names=NA, quote=F)
df <- data.frame(Dissimilarity=sol$diss, Distance=sol$dist)
write.table(df, paste0(name, '_stressplot.xls'), sep='\t', row.names=FALSE, quote=F)
pdf(paste0(name, '_stressplot.pdf'))
stressplot(sol)
dev.off()

si <- as.data.frame(sol$points)
stress <- round(as.numeric(sol$stress),5)

legend <- 'right'
color <- 'Samples'
si$Samples <- rownames(si)
n <- nrow(si)
if (length(unique(sample_list[,1])) == nrow(sample_list)) {
  si$Samples <- sample_list[match(rownames(si), sample_list[,2]), 1]
} else if (length(unique(sample_list[,1])) > 1) {
  si$Groups <- sample_list[match(rownames(si), sample_list[,2]), 1]
  si$Groups <- factor(si$Groups, levels=unique(si$Groups))
  color <- 'Groups'
  legend <- 'right'
  n <- nlevels(factor(si$Groups))
}
legend.title <- NULL
mean.point <-  TRUE

ellipse <- TRUE
ellipse.border.remove <- FALSE
label <- 'Samples'
pointch <- rep(15:18, 15)
if (n >= 50) {
  sangon <- rep('#CA0000', n)
  pointch <- rep(15, n)
}

p <- ggscatter(data.frame(si), x='MDS1', y='MDS2', color=color, shape=color, sangon=sangon, ellipse=ellipse, ellipse.level=0.8, ellipse.type='confidence', ellipse.alpha=0.1, mean.point=FALSE, star.plot=FALSE, ellipse.border.remove=ellipse.border.remove, repel=TRUE, label=label, legend.title=legend.title, legend=legend, xlab='NMDS1', ylab='NMDS2', title=paste('Stress :', stress), show.legend.text = FALSE) + scale_shape_manual(values=pointch[1:n])
p <- p +  border() + geom_vline(xintercept=0, linetype=2, colour='grey') + geom_hline(yintercept=0, linetype=2, colour='grey') + theme(plot.title = element_text(hjust=0.5))
p