NMDS分析和绘图

简介：本节为大家介绍一种常用的排序方法NMDS分析，全称为非度量多维尺度分析 (non-metric multidimensional scaling)。NMDS是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。目前网络上关于NMDS分析的介绍和论述已经很多了，本公众号就不再赘述了。本节的关注点是绘制漂亮的NMDS散点图。

1、哪些数据适宜NMDS分析？

也许大家已经发现NMDS分析大多情况下是用来展示物种数据的一种分析方法，为什么呢？这是因为当前流行使用的物种数据大多数为OTU或ASV测序数据，这类数据包含丰富的0值。NMDS分析一般使用的是Bray-Curtis距离算法，该算法对0值不敏感，换句话说即使有很多的0值的情况下也能获得较为稳健的结果。因此，NMDS分析适宜于测序数据的分析。然而，因为该方法使用的是非参数的方法，所以不能给出每一轴对于数据分布的解释量（如有错误请指正），这是该方法存在的局限性。

加载NMDS分析所需要的包，如下：

library(tidyverse)

library(vegan)

library(MASS)

library(readxl)

数据导入

readxlsx <- function(file = "file.xlsx", n =3) {

require(readxl)

dat <- list()

i = 0

while (i < n) {

i = i+1

dat[[i]] <- read_excel(file, sheet = i, col_names = T)

}

return(dat)

}

otu <- readxlsx(file = "1-16S.xlsx", n =5)

数据处理

ITS <- otu[[4]] %>%

data.frame()#tibble data change to dataframe data

rownames(ITS) <- ITS$OTUid#defined row names

ITS2 <- ITS[,-1] %>%

t() %>% #transposition

data.frame()#row names are sample names, and colnum names are OTU id

head(ITS2[,1:6])

最终分析的数据格式如下：

ITS.nmds<-metaMDS(ITS2)

ITS.nmds#The smaller the value of stress, the better the goodness of fit

当stress > 0.2时表明使用该方法不合适，建议使用其它方法对数据进行分析

stressplot(ITS.nmds)#检查观测值非相似性与排序距离之间的关系

拟合结果显示，没有点出现在距离线段很远的位置，意味着该数据可以使用NMDS分析。

简单出图

1). 只显示样方点

ordiplot(ITS.nmds, type = "text",display = "sites")

2). 只显示物种信息

ordiplot(ITS.nmds, type = "text",display = "species")

2、ggplot2绘图

提取样方点

ITS.scores <- as.data.frame(scores(ITS.nmds)) #提取点

ITS.scores %>%

as_tibble(rownames = "sample") ->ITS_sites

ITS_sites

根据处理给数据分组

otu[[3]] %>%

dplyr::select(Code,Tdiff) %>%

mutate(group = if_else(Tdiff>0,"warmer",

if_else(Tdiff<0,"colder","in_situ"))) ->group

group

将分组信息添加到NMDS数据样点中

ITS_sites %>%

left_join(group, by = c("sample" = "Code")) %>%

filter(group!="NA")->ITS_sites2

ITS_sites2

因为我用的是已发表文章中的数据，数据给出的样方信息和分组信息数量不匹配所以我过滤掉了不匹配的部分，如果处理自己的数据则不必使用filter函数。

绘图背景等参数设置（直接粘贴并运行）

main_theme = theme(panel.background=element_blank(),

panel.grid=element_blank(),

axis.line.x=element_line(size=0.5, colour="black"),

axis.line.y=element_line(size=0.5, colour="black"),

axis.ticks=element_line(color="black"),

axis.text=element_text(color="black", size=12),

legend.position="right",

legend.background=element_blank(),

legend.key=element_blank(),

legend.text= element_text(size=12),

text=element_text(family="sans", size=12),

plot.title=element_text(hjust = 0.5,vjust=0.5,size=12),

plot.subtitle=element_text(size=12))

绘图

ggplot(data=ITS_sites2,aes(NMDS1,NMDS2)) +

geom_hline(aes(yintercept=0),colour="#d8d6d6",linetype=5)+

geom_vline(aes(xintercept=0),colour="#d8d6d6",linetype=5)+

geom_point(aes(color = group),shape = 19,size = 3.5)+

scale_color_manual(values = c("#2166ac","#f4a582","#e31a1c"))+

#scale_x_continuous(breaks = seq(-0.59,0.66,0.2),limits = c(-0.59,0.66))+

#scale_y_continuous(breaks = seq(-0.60,0.45,0.15),limits = c(-0.60,0.45))+

labs(x= "NMDS1", y = "NMDS2",color = "Treatments")+

theme_bw() +

main_theme

3、分组NMDS和“等温线”NMDS

分组NMDS

ggplot(data=ITS_sites2,aes(NMDS1,NMDS2)) +

geom_hline(aes(yintercept=0),colour="#d8d6d6",linetype=5)+

geom_vline(aes(xintercept=0),colour="#d8d6d6",linetype=5)+

geom_point(aes(color = group),shape = 19,size = 3.5)+

scale_color_manual(values = c("#2166ac","#f4a582","#e31a1c"))+

#scale_x_continuous(breaks = seq(-0.59,0.66,0.2),limits = c(-0.59,0.66))+

#scale_y_continuous(breaks = seq(-0.60,0.45,0.15),limits = c(-0.60,0.45))+

stat_ellipse(aes(fill=group),geom="polygon",level=0.95,alpha=0.15)+

labs(x= "NMDS1", y = "NMDS2",

color = "Treatments",fill = "Treatments")+

theme_bw() +

main_theme

“等温线”NMDS

“等温线”NMDS适用于处理比较多的情况，如梯度等

1）重新构建分组

ITS_sites2 %>%

mutate(group2=if_else(Tdiff< -5.7,"very cold",

if_else(Tdiff< 0,"cold",

if_else(Tdiff<5.7, "in situ",

if_else(Tdiff <9.6, "warm","hot"))))) ->ITS_sites3

ITS_sites3

ggplot(data=ITS_sites3,aes(NMDS1,NMDS2)) +

geom_hline(aes(yintercept=0),colour="#d8d6d6",linetype=5)+

geom_vline(aes(xintercept=0),colour="#d8d6d6",linetype=5)+

geom_point(aes(color = group2),shape = 19,size = 3.5)+

#scale_color_manual(values = c("#2166ac","#f4a582","#e31a1c"))+

#scale_x_continuous(breaks = seq(-0.59,0.66,0.2),limits = c(-0.59,0.66))+

#scale_y_continuous(breaks = seq(-0.60,0.45,0.15),limits = c(-0.60,0.45))+

stat_density2d(aes(color = group2),size = 0.6)+

labs(x= "NMDS1", y = "NMDS2",

color = "Treatments")+

theme_bw() +

main_theme