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做一张完善的图

  • #首先提取前两轴坐标

  • 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本质是基于排名的算法更加适合NMDS

  • anosim.result<-anosim(vare.dis, design$SampleType,permutations =999)

  • 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()

成图展示：

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

请看下一篇文章

赞 (0)

相关推荐