NMDS非度量多维尺度分析—基于微生物群落

今天，看到赖江山老师在博客中分享了vegan中的一些函数的中文帮助文件，翻译专业，可读性强，这本材料是我们熟悉vegan原理和提高内涵的有力学习途径。（末尾有彩蛋）

本文主要做NMDS分析并做一张完善的高质量图片，提取stress值，推荐适合NMDS结果的差异分析并通过命令展示在图形上，最后加上置信区间椭圆。

非度量多维尺分析（ NMDS）是一种很好的排序方法，因为它可以使用 具有生态学意义的 方法来度量群落差异 。一个好的 相异性测度与环境梯距离具有很好的秩 关系。因为 NMDS只使用秩信息，并且映射的在有序空间 上是非线性的， 故它能处理任意 类型 的非线性物种 矩阵 ，并能有效、稳健地找 到潜在梯度。

NMDS分析，网络上已近有很多相关教程分享其原理，与其他排序（PCA、PCoA、CCA、RDA） 方法的不同之处，简单来讲NMDS也是一种使用物种组成数据的排序称作非限制性排序；NMDS基于距离算法，优于PCA、PCoA、CCA、RDA的地方在于当样本或者物种数量过多的时候使用NMDS会更加准确；

vegan 的ordiplot()函数可以用来绘制NMDS 的结果：

plot(vare.mds, type = "t")

vegan 包中的metaMDS()函数不需要单独计算相异矩阵，直接 将原始数据矩阵作为输入。结果比以前更丰富 ，除了奥杜尔包中isoMDS()结果中 的成分外还有很多其他结果输出：nobj, nfix, ndim, ndis, ngrp, diss, iidx, jidx, xinit, istart, isform, ities, iregn, iscal, maxits, sratmx, strmin, sfgrmn, dist, dhat, points, stress, grstress, iters, icause, call,model, distmethod, distcall, data, distance, converged, tries,engine, species。该函数将这些的过程封装到一个函数中：

  1. 一般生态群落数据比较离散 ，用平方根转换数据 ，然后进行 Wisconsin双重 标准化，或物种除以它们的最大值 将数据 均一化 为相等的总数。这两个标 准化通常可以提高排序的质量，但是我们在最初分析中忘了考虑数据的转化问题 。

  2. 默认使用 Bray-Curtis相异系数。

  3. 运行多次独立的isoMDS()，并在一定次数的尝试之后停止，或者找到两个具有最小应力 函数之后停止 ，返回了最佳的排序结果。

  4. 旋转排序图，使样方坐标的最大差位于第一轴上。

  5. 对排序结果进行标准化，使一个单元应于将群落相似性从重复相似性减半。

  6. 函数发现物种 排序轴为样方排序轴 的加权平均值 并将其扩大，使物种和样 方排序轴 具有相等的方差，可以使用shrink = TRUE撤消。metaMDS()的帮助页面将提供更多细节，并解释函数使用 的帮助页面将提供更多细节，并解释函数使用 过程 。

清理环境

#清空内存
rm(list=ls())

准备主题和数据

Mytheme <- theme_bw() +
#scale_fill_brewer(palette = "YIOrRd", guide = guide_legend(title = NULL), limits = c("CK1","CK3","CK5","CK7","CK9","CK11","CK13","CK15","CK17","CK19"))+
theme(

panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),

plot.title = element_text(vjust = -8.5,hjust = 0.1),
axis.title.y =element_text(size = 24,face = "bold",colour = "black"),
axis.title.x =element_text(size = 24,face = "bold",colour = "black"),
axis.text = element_text(size = 20,face = "bold"),
axis.text.x = element_text(colour = "black",size = 14),
axis.text.y = element_text(colour = "black",size = 14),
legend.text = element_text(size = 15,face = "bold")
#legend.position = "none"#是否删除图例

)

#设定路径
path = getwd()

# 导入包
library(phyloseq)
library(ggplot2)
suppressMessages(library(vegan))

# 使用示例数据，注意是phyloseq封装好的
data("GlobalPatterns")
ps = GlobalPatterns

提取数据 运算NMDS

vegan_otu <- function(physeq){
OTU <- otu_table(physeq)
if(taxa_are_rows(OTU)){
OTU <- t(OTU)
}
return(as(OTU,"matrix"))
}
x = as.data.frame(t(vegan_otu(GlobalPatterns)))
head(x)
x = as.matrix(x)

x = t(t(x)/colSums(x,na=T))* 100 # normalization to total 100
head(x)

##bray
bray.mds<-metaMDS(t(x), distance="bray", k=2, trymax=100) #maximum numbers of random starts in search of stable solution
bray.mds

##jackard
x = decostand(x,"pa")
jaccard.mds<-metaMDS(t(x), distance="jaccard", k=2, trymax=100)
jaccard.mds

str(bray.mds) #structure
# ## 输出坐标
bray_axis = bray.mds$points
jaccard_axis = jaccard.mds$point

计算Stress值

Stress值是反映模型合适程度的指标，NMDS会多次打乱数据计算Stress值，直到找到最合适的模型，也就是最低的Stress值；理想状况下，Stress值为0，一般Stress值低于0.1较为合理（本数据这个值偏高一些）

# 读入实验设计和Alpha多样性值
design = as.data.frame(sample_data(ps))
head(design)

########outbray出图坐标准备
outbray = as.data.frame(bray_axis)
index = merge(outbray,design, by="row.names",all=F)
head(index)
dim(index)

stress = paste("bray ","stress: ",round(bray.mds$stress,3), sep = "")
stress

使用坐标和stress出图

mi = c("#FFF5EB" ,"#FEE6CE" ,"#FDD0A2", "#FDAE6B", "#FD8D3C", "#F16913", "#D94801", "#A63603", "#7F2704","black")
# mi=c("#1B9E77" ,"#D95F02", "#7570B3","#E7298A")
p <-ggplot(index, aes(x=MDS1, y=MDS2, fill = SampleType)) +
geom_point(alpha=.7, size=5, pch = 21) +
labs(x=paste("NMDS1",sep=""),
y=paste("NMDS2" ,sep=""),
title=stress)+
#stat_ellipse( linetype = 2,level = 0.65,aes(group =group, colour = group))+
#stat_ellipse( linetype = 1,level = 0.8)+
#geom_text_repel(aes(label=points$id),size=4)+
scale_fill_manual(values = mi)+
#labs(title = "toamto hea and dis")+
guides(color=guide_legend(title = NULL),shape=guide_legend(title = NULL))+
#scale_y_continuous(expand = c(0,0))+
geom_hline(aes(yintercept=0), colour="black", linetype=2) +
geom_vline(aes(xintercept=max(index$MDS2/2)), colour="black", linetype="dashed")
p
# points$id=row.names(points)
# p+geom_text(aes(label=points$id),size=4)#?stat_ellipse
p = p + Mytheme
p

plot_name = paste(path,"/a3_NMDS.pdf",sep = "")
ggsave(plot_name, p, width = 8, height = 6)

成图展示：

基于ade4 的NMDS分析（此中文来自赖江山老师团队翻译结果）

非度量多维尺分析（ Non-metric multidimensional scaling, NMDS）可以用 ）可以用 MASS包中的 isoMDS()函数实现，输入相异矩阵即可。函数实现，输入相异矩阵即可。vegan的 vegdist()函数可 以计算群落 的相异矩阵 ，默认是 Bray-Curtis相异系数，现在通常称为 Steinhaus 相似指数，在芬兰称为 Sorensen指数。基本步骤如下 指

library(vegan)
library(MASS)

data(varespec)
vare.dis <- vegdist(varespec)
vare.mds0 <- isoMDS(vare.dis)

NMDS排序结构通过迭代来不断最小化应力函数（stress function），默认情 况下是找到两个维度并使用度量尺度分析(cmdscale)作为初始结构进行调整。从 跟踪(trace)信息中可以看出迭代过程(通过设置参数trace = F来隐藏迭代过程)。

isoMDS()返回一个排序构建过程和应力函 数的列表(item points, stress)。S应力函数是一个拟合度统计量，是排序空 间内对象结构与原始距离矩阵之间的相异程度 的度量。NMDS将观察到的群落差异非线性地映射 到排序空间上，可以基于任何类型距离矩阵对 对象进行排序。可以用MASS包的函数 Shepard()或者vegan包的stressplot()函数来评估 NMDS的结果（Shepard图）。

stressplot(vare.mds0, vare.dis)

stressplot函数绘制了一个Shepard图，其中 横坐标为原始距离，纵坐标为排序距离，用单 调的折线拟合。此外，stressplot()显示了这两者 距离相关性，如拟合度（goodness of fit）与应 力函数的关系是R²= 1 - S²。“fit-based R²”是拟合 值θ(d)和运算出的排序图上距离d之间的相关 性，或者是折线和点之间的相关性。它应该是 线性的，即使拟合有点弯曲，通常仍被称为“线性拟合”。这两个相关性都是基 于Shepard图中的残差，但是它们的零模型有所不同。在线性拟合中，零模型是 所有排序距离相等，拟合为一条水平直线。这听起来很合理，但是需要N-1维的 N个点的零模型，而这个零模型在排序空间中是没有几何意义的。基本应力采 用零模型，所有的观测都放在同一点上，这在几何上是可能的。注意，有时人 们使用群落差异和排序距离之间的相关性。但是由于NMDS是一种非线性方 法，因此这样做既危险又具有误导性：使用该准则，具有更多非线性关系的分类将会出现更多错误。

reference

http://wap.sciencenet.cn/home.php?mod=space&uid=267448&do=blog&id=1194081

赞 (0)

相关推荐