導讀

非度量多維尺度分析（nonmetric multidimensional scaling, NMDS），是一種簡介的梯度分析方法，也是基於距離或者相異性矩陣。與其它主要用於最大化變異和一致性的方法不一樣，NMDS是一種排序方法。這對於距離缺失的資料有優勢，只要先辦法確定物件之間的位置關係，便可以進行NMDS分析。NMDS的計算過程會以程式碼的形式貼在下方，供大家參考

資料格式

NMDS與PCoA一樣，NMDS可以基於任何型別距離、相異性矩陣物件（樣方）進行排序。當然也可以是原始資料矩陣。這裡我用的是weighted unifrac距離矩陣資料

分析程式碼

NMDS排序分析可以通過生態學分析R包vegan中的metaMDS()函式實現。因為輸入metaNMDS()的資料可以是原始資料矩陣，也可以是距離矩陣，這裡拿上面列舉的資料做示範。

執行NMDS分析

rm(list = ls())
path <- "D:/元基因組/3. 16S 擴增子下游測序分析&其他分析方法集錦/土壤菌群與水稻土溶解有機質/16S"
setwd(path)
library(vegan) ## 載入包
weight_dm <- read.table("weighted_unifrac_otu_table_css.txt",header = T,row.names = 1,
                        sep = "\t",check.names = F) ## 載入weighted距離矩陣資料
meta_info <- read.csv("meta_info.csv",header = T,row.names = 1,check.names = F) ## 載入樣品分組資訊
#group <- meta_info$group 
set.seed(1234) ## 設定隨機種子，以便結果可以重複
weight_nmds <- metaMDS(weight_dm,trace = F) # trace = F 表示的是不要在螢幕上輸出執行的過程和結果
weight_nmds$stress # 壓力值，一般小於0.1比較好，但是也要根據所選擇的的主成分數目來看

圖片來源GUSTA ME[2]

作圖

tiff("NMDS.tiff",width=1000,height=1000)
p<- plot(scores(weight_nmds, choice=c(1, 2)), col=c("purple","green","blue","red")[meta_info$group], 
         cex=1.5, font=1,  pch=5) ##NMDS作圖###
legend('topleft', legend=levels(meta_info$group), col=c("purple","green","blue","red"),pch=5,cex=1.5,box.lty =1)###NMDS新增Legend###
with(weight_nmds,ordiellipse(weight_nmds, meta_info$group,display = "sites", kind = "se", conf = 0.95, lwd=2,cex=0.8,lty=1,col=c("purple","green","blue","red"),font=2,label = FALSE))
#dev.off()
#env<-meta_table[,c(5:10)]
#head(env)
#ef1<-envfit(otu.nmds,env,na.rm=TRUE)###環境變數envfit###
#p<- p+ plot(ef1,p.max=0.05,col="black") #####plot(ef1)##所有環境因子圖新增箭頭###
#ef1
dev.off()

site.scores <- as.data.frame(scores(weight_nmds)) ## 抽提出樣本和NMDS1，NMDS2資訊
site.scores <- cbind(site.scores,group)
head(site.scores)
library("FactoMineR")
library("factoextra")
p <- ggplot(data = site.scores,aes(x=NMDS1,y=NMDS2)) + 
  geom_point(aes(shape=group,color=group,fill=group),size=3) 
  
p <- p + theme_bw()
p <- p + theme(panel.grid.major = element_blank(), plot.background=element_blank(),
         panel.grid.minor = element_blank(), 
         legend.position="top", 
         axis.title.x = element_text(size = 16), 
         axis.text.x = element_text(angle=0,color="black",vjust = 0.95,hjust = 0.95,size=12), 
         axis.title.y = element_text(size=16),
         axis.text.y = element_text(size=14,color="black"), 
         strip.text.x = element_text(size=18), 
         legend.text = element_text(size=14), 
         legend.title = element_text(size=16))
p<-p+geom_vline(xintercept=0,linetype="dashed",color="blue")+geom_hline(yintercept = 0,linetype="dashed",color="blue")
p <- p + geom_text(x=0.12,y=-0.2,label=paste("stress = ",round(weight_nmds$stress,3),sep = ""),color="blue",size=6)
ggsave("NMDS.tiff", height=8, width=8, units="in")

參考

[1] 中文版 《數量生態學-R語言的應用》（賴江山 譯）高等教育出版社出版
[2] https://mb3is.megx.net/gustame/dissimilarity-based-methods/nmds