【R语言】非度量多维标度分析法(Non-metric multidimensional scaling，NMDS)

Python016

【R语言】非度量多维标度分析法(Non-metric multidimensional scaling，NMDS),第1张

非度量多维尺度分析（NMDS 分析）是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法，能够反映对象间的顺序关系。与PCoA类似，NMDS可以基于任何类型距离矩阵对对象（样方）进行排序；但也有不同之处，NMDS不在基于距离矩阵数值，而是根据排位顺序进行计算。对于存在距离缺失的数据而言有优势，只要想办法确定对象间的位置关系，即可进行NMDS分析。由于NMDS不是特征根排序技术，也不再让排序轴载更多的变差目的；因此NMDS排序图可以任意旋转、中心化和倒置（赖江山, 数量生态学）。 预先设定排序轴的数量m； 在m维空间内构建对象的初始结构，初始结构是调整对象之间位置关系的起点； 在m维空间内，用一个迭代程序不断调整对象位置，目标是不断最小化应力函数（Stress function，其值被转化为0~1间的数值，可以检验 NMDS 分析结果的优劣。通常认为 stress<0.2 时可用 NMDS 的二维点图表示，其图形有一定的解释意义；当 stress<0.1 时，可认为是一个好的排序；当 stress<0.05 时，则具有很好的代表性）不断调整对象位置，直至应力函数值不再减少，或已达到预定的值； 大部分NMDS会根据PCA结果旋转最终的排序图，使结果更容易解读。 Shepard 图：能够比较NMDS中对象间距离与原始距离测度矩阵中的值或者通过二者间进行线性（Linear fit）或非线性（Non-metric fit）拟合的R^2 评估，拟合R^2 越大越好/两个R^2 值越一致越好。 1、Shepard图的R^2 ：拟合R^2 越大越好/两个R^2 值越一致越好； 2、Stress值：Stress<0.2（有一定可靠性）；Stress<0.05（结果较好）；Stress<0.02（结果很好）；Stress<0.01（结果极好）。 注：stress=0.12，说明结果不好，但有一定的可靠性。 图形中的点代表样本，不同颜色/形状代表样本所属的分组信息。同组样本点距离远近说明了样本的重复性强弱，不同组样本的远近则反应了组间样本距离在秩次（数据排名）上的差异。样本相似性距离计算方式对结果有影响，选择输入不同相似性距离值的矩阵，得到的结果存在着不同程度差异。 横纵坐标轴含义：NMDS是距离值的秩次（数据排名）信息的评估，图上样本信息仅反映样本间数据秩次信息的远近，而不反映真实的数值差异，横纵坐标轴并无权重意义，横轴不一定比纵轴更加重要。 NMDS是非参数的方法，不应该有解释量，不知道他们这些NMDS的解释率是怎么算出来的？我估计是用主坐标分析(PCoA)的解释量来作为NMDS的解释量，但这种张冠李戴的做法可能是不合适的（引自赖江山老师： http://blog.sciencenet.cn/blog-267448-1146112.html ）。

非度量多维标度（NMDS）分析

非度量多维标度（Non-metric Multidimensional Scaling，NMDS）是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。

与PCoA相似，NMDS都使用样本相似性距离矩阵进行降维排序分析，但NMDS侧重反映距离矩阵中数值的排序关系，NMDS图形常用作微生物群落研究的β分析。

来看看R的vegan包和ggplot2包是如何绘制NMDS图的。

使用数据文件如下：

①otu_table.txt

②group.txt

方法如下：

1.调用vegan包metaMDS()函数执行NMDS，读取数据；

2.根据物种组成计算样方距离，如 Bray-curtis 距离；

3.NMDS 排序，定义 2 个维度；

4.提取应力函数值、样方得分，查看主要结果；

5.物种变量可通过丰度加权平方被动添加至排序图中；

6.使用ggplot2包作图。

计算β多样性指数需要用到phyloseq包。它的安装方式不同于简单的install.packages（“phyloseq”）

有两种方法可以安装

1.先安装BiocManager

install.packages("BiocManager")

library("BiocManager")

BiocManager::install("phyloseq")

library("phyloseq")

2.source("https://bioconductor.org/biocLite.R")

biocLite("phyloseq")

#安装phyloseq

library("phyloseq")

安装并加载了phyloseq包后，开始读取数据，前面计算α多样性，用到的是read.table……

qiimedata <- import_qiime(otufilename = "feature-table.taxonomy.txt", mapfilename = "mapping_file.txt", treefilename = "tree.rooted.nwk", refseqfilename = "dna-sequences.fasta")

#读取数据，参数都是文件名，注意加后缀

#otufilename指定out表格，mapfilename指定map文件（分组数据)

#treefilename指定有根进化树文件

#refseqfilename指定代表序列文件

otu<-qiimedata@[email protected]

#从qiimedata里面提取otu

sum_of_otus<-colSums(t(otu))

#t_转置,colsums计算列的和,即计算各个otu检测到的总序列数，为了筛掉一些总序列数过低的otu（可能是测序错误）

sum_of_otus

#查看otu总序列数

selected_otu<-names(sum_of_otus)[sum_of_otus>10]

#获取总序列数大于10的otu id

sub_qiimedata <- prune_taxa(selected_otu, qiimedata)

#筛选总序列数大于10的otu的phyloseq数据

weighted_unifrac<-distance(sub_qiimedata,method = 'wunifrac')

#计算样本间加权unifrac

unweighted_unifrac<-distance(sub_qiimedata,method = 'unifrac')

#计算样本间非加权unifrac

bray_curtis <- distance(sub_qiimedata, method='bray')

write.table(as.matrix(bray_curtis),"bray_curtis.txt",sep = '\t',quote = FALSE,col.names = NA)

#保存距离矩阵

#计算样本间Bray-Curtis距离矩阵，method 可选" wunifrac ", " unifrac " ，"jaccard"等

pcoa_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "PCoA")

#基于Bray-Curtis距离矩阵的PCoA排序分析

p<-plot_ordination(sub_qiimedata, pcoa_of_bray_curtis, type="samples", color="Group1",shape = "Group1")

#将PCoA排序分析结果可视化

library("ggplot2")

p<-p+ scale_colour_manual(values=c("#DC143C","#808000","#00CED1")) + geom_point(size=2) +ggtitle("PCoA of Bray-Curtis distance")+theme(text = element_text(size = 15))

#修改图形大小,ggtitle加标题,stat_ellipse加椭圆

#用scale_colour_manual(values=c())自定义颜色，可查颜色的16进制对照表

p

nmds_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "NMDS")

#基于Bray-Curtis距离矩阵的NMDS排序分析

p1<-plot_ordination(qiimedata, nmds_of_bray_curtis, type="samples", color="Group1")

#将NMDS排序分析结果可视化

# color=“Group1”指定不同分组的点染不同颜色

p1

p1<-p1+ geom_point(size=3) +ggtitle("NMDS of Bray-Curtis distance") + stat_ellipse()+theme(text = element_text(size = 15))

#对图片进行适当修饰， stat_ellipse()加椭圆， ggtitle()加标题

ggsave(plot = p1,“nmds_of_bary_curtis.pdf",dpi = 300,width

PCoA中的两个点距离，接近β多样性指数

PCA(Principal Components Analysis)即主成分分析，也称主分量分析或主成分回归分析法，首先利用线性变换，将数据变换到一个新的坐标系统中然后再利用降维的思想，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上。这种降维的思想首先减少数据集的维数，同时还保持数据集的对方差贡献最大的特征，最终使数据直观呈现在二维坐标系。

PCoA(Principal Co-ordinates Analysis)分析即主坐标分析，可呈现研究数据相似性或差异性的可视化坐标，是一种非约束性的数据降维分析方法，可用来研究样本群落组成的相似性或相异性。它与PCA类似，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样本点之间的相互位置关系，只是改变了坐标系统。两者的区别为PCA是基于样本的相似系数矩阵(如欧式距离)来寻找主成分，而PCoA是基于距离矩阵(欧式距离以外的其他距离)来寻找主坐标。

NMDS图中两个点的距离的排序，接近β多样性指数的排序