NMDS分析
哪些数据适宜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 idhead(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")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分组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_sites3ggplot(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链接：https://www.aliyundrive.com/s/CFvkjJ3nECi
如有问题，可以加入我们的群聊一起讨论，如下：
赞 (0)