導讀
非度量多維尺度分析（nonmetric multidimensional scaling, NMDS）统抬，是一種簡介的梯度分析方法冀宴，也是基于距離或者相異性矩陣灭贷。與其它主要用于最大化變異和一致性的方法不一樣，NMDS是一種排序方法花鹅。這對于距離缺失的數據有優(yōu)勢氧腰，只要先辦法確定對象之間的位置關系，便可以進行NMDS分析刨肃。NMDS的計算過程會以代碼的形式貼在下方古拴，供大家參考
數據格式
NMDS與PCoA一樣，NMDS可以基于任何類型距離真友、相異性矩陣對象（樣方）進行排序黄痪。當然也可以是原始數據矩陣。這里我用的是weighted unifrac距離矩陣數據
圖1 weighted 距離矩陣
分析代碼
NMDS排序分析可以通過生態(tài)學分析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比較好站绪，但是也要根據所選擇的的主成分數目來看
圖2 strees的scree圖遭铺，stress隨著主成分數的增加而減少
圖片來源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)###環(huán)境變量envfit###
#p<- p+ plot(ef1,p.max=0.05,col="black") #####plot(ef1)##所有環(huán)境因子圖添加箭頭###
#ef1
dev.off()
圖3 常規(guī)方法作圖
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")
圖4 ggplot2作圖
參考
[1] 中文版 《數量生態(tài)學-R語言的應用》（賴江山 譯）高等教育出版社出版
[2] https://mb3is.megx.net/gustame/dissimilarity-based-methods/nmds
