导读

非度量多维尺度分析（nonmetric multidimensional scaling, NMDS），是一种简介的梯度分析方法，也是基于距离或者相异性矩阵。与其它主要用于最大化变异和一致性的方法不一样，NMDS是一种排序方法。这对于距离缺失的数据有优势，只要先办法确定对象之间的位置关系，便可以进行NMDS分析。NMDS的计算过程会以代码的形式贴在下方，供大家参考

数据格式

NMDS与PCoA一样，NMDS可以基于任何类型距离、相异性矩阵对象（样方）进行排序。当然也可以是原始数据矩阵。这里我用的是weighted unifrac距离矩阵数据

图1 weighted 距离矩阵

分析代码

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比较好，但是也要根据所选择的的主成分数目来看

图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)###环境变量envfit###
#p<- p+ plot(ef1,p.max=0.05,col="black") #####plot(ef1)##所有环境因子图添加箭头###
#ef1
dev.off()

图3 常规方法作图

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] 中文版 《数量生态学-R语言的应用》（赖江山 译）高等教育出版社出版
[2] https://mb3is.megx.net/gustame/dissimilarity-based-methods/nmds