摘要：beta多样性分析网络上已经有很多，大部分出图都是基于R的基础包，很难看，而且很多图形不够完善，没有做方差分析，更没有将其直接输到图上，本文主要做NMDS分析并做一张完善的高质量图片，提取stress值，推荐适合NMDS结果的差异分析并通过命令展示在图形上，最后加上置信区间椭圆。

NMDS分析，网络上已近有很多相关教程分享其原理，与其他排序(PCA、PCoA、CCA、RDA)

方法的不同之处，简单来讲NMDS也是一种使用物种组成数据的排序称作非限制性排序；NMDS基于距离算法，优于PCA、PCoA、CCA、RDA的地方在于当样本或者物种数量过多的时候使用NMDS会更加准确；

基于R语言，开始作图：安装package，两个重要的包

install.packages("vegan")#微生态常用包，一定安装上

install.packages("ggplot2")

调用package

library("vegan")

#设定工作路径

setwd("E:/马兄作图")

#读取数据，一份otu.table文件和一份分组信息文件

design

otu=read.table("16s_otu_tablejidai.txt",row.names=1, header=T,sep="\t")

#如果数据调整为列名是OTU，行名是样本名

otu=t(otu)

#查看一下

head(otu)

#标准化，当然method有很多，可以跟换?decostand查看其它method

vare.hel

#计算矩阵

vare.dis

#使用NMDS的方法

vare.mds

到此分析过程就算完成了，使用基础包plot()作图当然已不是我们最佳的选择了，图形很基础，我就不展示了，下面基于ggplot做一张完善的图#首先提取前两轴坐标

point = scores(vare.mds)

#将分组文件和得分文件合并

index = merge(design, point,by="row.names",all=F)

计算Stress值

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

vare.mds

结果输出中寻找下面stress结果：

#显着性检验；anosim本质是基于排名的算法更加适合NMDSanosim.result

summary(anosim.result)

结果输出，得到下面这两个值：

#tiff输出图形，适合大部分出版刊物，入门级别分辩率300,18*14的长宽；

tiff(file="beta_bray_NMDS.tif", res = 300, compression ="none", width=180,height=140,units= "mm")

#开始出图，将上面得到的三个指标在图中更换stress，R，p，不多说，代码如下：

library("ggplot2")

p = ggplot(index, aes(x=NMDS1, y=NMDS2, color=SampleType)) +

geom_point(alpha=.7, size=2) +

labs(x=paste("NMDS1"),

y=paste("NMDS1"),

title="")#置信区间当然要加上，有三种方式，线条类型也可以更改

p+stat_ellipse(type = "t", linetype = 2)+

annotate("text",x=-1.07,y=-1,parse=TRUE,size=4,label=""R:"*0.7031",family="serif",fontface="italic",colour="darkred")+

annotate("text",x=-1.1,y=-0.9,parse=TRUE,size=4,label=""p:"*0.001",family="serif",fontface="italic",colour="darkred")+

annotate("text",x=-1,y=-0.78,parse=TRUE,size=4,label=""Stress:"*0.1208",family="serif",fontface="italic",colour="darkred")dev.off()

成图展示：

当然这还达不到对于拥有轻微强迫症的我对于图行要求：

请看下一篇文章

导读

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

数据格式

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

图1 weighted 距离矩阵

分析代码

NMDS排序分析可以通过生态学分析R包vegan中的metaMDS()函数实现。因为输入metaNMDS()的数据可以是原始数据矩阵，也可以是距离矩阵，这里拿上面列举的数据做示范。

运行NMDS分析

rm(list = ls())

path

setwd(path)

library(vegan) ## 加载包

weight_dm

sep = "\t",check.names = F) ## 加载weighted距离矩阵数据

meta_info

#group

set.seed(1234) ## 设置随机种子，以便结果可以重复

weight_nmds

weight_nmds$stress # 压力值，一般小于0.1比较好，但是也要根据所选择的的主成分数目来看

图2 strees的scree图，stress随着主成分数的增加而减少

图片来源GUSTA ME[2]

作图

tiff("NMDS.tiff",width=1000,height=1000)

p

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

#head(env)

#ef1

#p

#ef1

dev.off()

图3 常规方法作图

site.scores

site.scores

head(site.scores)

library("FactoMineR")

library("factoextra")

p

geom_point(aes(shape=group,color=group,fill=group),size=3)

p

p

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

ggsave("NMDS.tiff", height=8, width=8, units="in")

图4 ggplot2作图

参考

今天，看到赖江山老师在博客中分享了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。该函数将这些的过程封装到一个函数中：一般生态群落数据比较离散 ，用平方根转换数据 ，然后进行 Wisconsin双重 标准化，或物种除以它们的最大值 将数据 均一化 为相等的总数。这两个标 准化通常可以提高排序的质量，但是我们在最初分析中忘了考虑数据的转化问题 。

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

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

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

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

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

清理环境#清空内存

rm(list=ls())

准备主题和数据Mytheme

#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

提取数据 运算NMDSvegan_otu

OTU

if(taxa_are_rows(OTU)){

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

bray.mds

##jackard

x = decostand(x,"pa")

jaccard.mds

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

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

vare.mds0

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