R-16s分析作图之NMDS非度量多维尺度分析

摘要：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<-read.table('design16s_otu_tablejidai.txt',header=T,sep="\t",row.names=1)
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<-decostand(otu,method="hellinger")
#计算矩阵
vare.dis <- vegdist(vare.hel,method="bray")
#使用NMDS的方法
vare.mds <- metaMDS(vare.hel,distance = "bray")

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

计算Stress值

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

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

#显著性检验；anosim本质是基于排名的算法更加适合NMDS

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

#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="")

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")

成图展示：

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

请看下一篇文章

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。