导读
非度量多维尺度分析(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