非度量多维排列 NMDS (Non-metric multidimensional scaling)分析

大家好，欢迎来到IT知识分享网。

多维排列 (Multidimensional scaling,MDS)是可视化多变量样品(如多个物种丰度、多个基因表达)相似性水平的一种方法。其基于距离矩阵进行一系列的排序分析。

经典的MDS (CMDS)分析就是前面提到的PCoA分析，也称为度量性MDS分析。而与之相对的是非度量多维排列 (Non-metric multidimensional scaling, NMDS)。

非度量多维排列 (NMDS)是基于相异矩阵或距离矩阵进行排序分析的间接梯度分析方法。NMDS的目标与PCA或PCoA类似 (一文读懂PCA分析 （原理、算法、解释和可视化）;一文学会PCA/PCoA相关统计检验(PERMANOVA)和可视化)，都是希望能在低维空间尽可能准确地展示样品在高维空间的关系。在低维空间，两个样品点距离越近，其相似度越高。NMDS的主要目的是识别和解释样品的分布模式，反应样品之间的顺序关系，找到能展示样品差异来源的梯度信息，如地理环境信息、生态信息等。

与MDS不同的是，NMDS分析将原始的距离矩阵转换为秩矩阵 (rank metric)再进行降维分析。NMDS弱化距离矩阵中具体值的大小，更关注其排序关系。假如样品A和样品B之间的距离是5，样品A和样品C之间的距离是10，转换后不再描述距离，而是说样品B是与样品A第1近的，样品C与样品A是第2近的，用排序的1,2代替原始的距离。所以称为“非参数”分析。

NMDS可以应用于 1) 存在配对距离缺失的数据中，2) 任何距离算法产生的矩阵，3) 定量、半定量、定性或混合变量的分析。

在多样本、物种数量多的情况下，NMDS模型能更准确地反映出距离矩阵的数值排序信息。因此当样本或者物种数量过多的时候使用NMDS会更加准确。

NMDS是不断迭代的过程，通过不断尝试找到样品在维度空间的最合适位置。其评估标准是stress值，表示 观察到的距离和拟合的距离的不一致性。Stress也可以理解为样品在降维后形成的空间的距离与其在原始多维空间的距离的差值。这个值越小越好，说明在低维空间更完整地捕获了高维空间的信息。对于NMDS二维分析，通常认为stress<0.2时有一定的解释意义；当stress<0.1时，可认为是一个好的排序；当 stress<0.05时，则具有很好的代表性。

NMDS在分析之前就会选择降维轴的数目并把数据拟合到所选的轴进行排序（轴越多，stress值就会越少；但轴越多，越难以解释）。NMDS不同于其它排序算法的地方在于其最终解不唯一，而是多次尝试获得一个能接受的解。NMDS算法不使用奇异值-奇异向量等因子分解技术，同时NMDS1、NMDS2也不一定是能解释最大差异的轴（不过后面分析时会使得第一轴解释的差异最大，以便更好可视化）。所以NMDS的轴可以按需转换。

在生物信息中，NMDS用于时间序列表达谱中鉴定基因的变化模式 (https://www.biorxiv.org/content/10.1101/v1.full)和宏基因组数据中分析微生物群落的差异。

metaMDS – NMDS分析的推荐步骤

NMDS通常被认为是生态学领域最稳定的非限制性排序方法 (Minchin 1987)。metaMDS是vagan中的一个组合了Minchin’s (1987) 推荐的分析步骤的NMDS分析函数，其包含下面这些步骤：

1. 数据转换：在参数autotransform=T(默认)时,如果输入的物种丰度矩阵（一般是抽平后的物种丰度矩阵）中最大的丰度值大于9则会进行Wisconsin双重标准化(每个物种丰度值除以该物种最大丰度，再在各个样品内计算相对丰度)；如果最大的丰度值大于50，则会对数据先进行开方处理再做Wisconsin转换。如果是自己转换过的数据，比如hellinger转换后的数据不想再被进一步转换，可以设置autotransform = FALSE。如果输入的是距离矩阵，这一步也会直接跳过。
2. 计算相异矩阵: 默认是Bray-Curtis，也通常是效果最好的。也可以选择vegdist输出的其它距离矩阵。对于非群体构成数据，可以用函数rankindex寻找自己的数据最合适的矩阵算法。如果输入的是距离矩阵，这一步也会直接跳过。
3. 最短路径插值 (Step-across dissimilarities): 如果大比例的样品没有共有物种，则排序将很难进行。在这一情况下，就需要通过计算这些不同样品间的最短路径代替样品之间的相异值或距离。默认的NMDS引擎函数monoMDS可以自动处理无共享物种的情况，不需要额外调用stepacross函数进行最短路径计算。而NMDS的另一个引擎函数isoMDS则不能自动处理这一情况，需要设置noshare=T触发stepacross步骤。如果输入的是距离矩阵，这一步也会直接跳过。
4. 多轮NMDS寻找最优解: NMDS会很容易陷入局部最优点，需要用不同的随机起始多运行几次才更可能获得全局最优解。metaMDS的策略是先运行PCoA分析并以其结果作为参考标准 (RUN 0)。如果设置了previous.best参数，则以该参数传入的NMDS结果作为参考。随后metaMDS会设定多个随机起始点运行NMDS分析 (参数try和trymax可以设置最小和最大尝试次数)。如果某一个NMDS的结果优于当前最优结果 (判断标准是：更低的stress值)，则该结果升级为当前最优结果，继续循环。可以设置trace = 2或更大的值跟踪这一优化过程。设置parallel=x设置进程数加速运算。
5. 结果优化: metaMDS在获得NMDS结果后，调用postMDS进一步优化结果：1) 结果整体向坐标轴中心聚拢; 2) 根据主成分旋转NMDS1使第一维解释的差异最大 (也可以调用函数MDSrotate旋转第一轴与指定的环境变量平行); 3) 群落相似性单位均分。
6. 物种得分: 在最终NMDS结果中用函数wascores计算物种的加权得分。

实战NMDS分析

继续使用之前的测试数据（如何读入自己的数据见前文和抄代码的时候总是遇到原始数据应该长什么样的问题）。

注意：测试数据已经做过了转置，每一行为一个样品，每一列为一个物种；自己的数据通常是每行为一个物种/OTU，每列为一个样品，读入后需要转置。

library(vegan) data("dune") data("dune.env") rownames(dune) <- paste0("Sample",1:nrow(dune)) # 已经做了转置，样品在行，物种在列 head(dune) Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere Bromhord Chenalbu Cirsarve Comapalu Eleopalu Elymrepe Empenigr Hyporadi Sample1 1 0 0 0 0 0 0 0 0 0 0 4 0 0 Sample2 3 0 0 2 0 3 4 0 0 0 0 4 0 0 Sample3 0 4 0 7 0 2 0 0 0 0 0 4 0 0 Sample4 0 8 0 2 0 2 3 0 2 0 0 4 0 0 Sample5 2 0 0 0 4 2 2 0 0 0 0 4 0 0 Sample6 2 0 0 0 3 0 0 0 0 0 0 0 0 0 Juncarti Juncbufo Lolipere Planlanc Poaprat Poatriv Ranuflam Rumeacet Sagiproc Salirepe Scorautu Trifprat Trifrepe Vicilath Sample1 0 0 7 0 4 2 0 0 0 0 0 0 0 0 Sample2 0 0 5 0 4 7 0 0 0 0 5 0 5 0 Sample3 0 0 6 0 5 6 0 0 0 0 2 0 2 0 Sample4 0 0 5 0 4 5 0 0 5 0 2 0 1 0 Sample5 0 0 2 5 2 6 0 5 0 0 3 2 2 0 Sample6 0 0 6 5 3 4 0 6 0 0 3 5 5 0 Bracruta Callcusp Sample1 0 0 Sample2 0 0 Sample3 2 0 Sample4 2 0 Sample5 2 0 Sample6 6 0 rownames(dune.env) <- paste0("Sample",1:nrow(dune.env)) # 通过metaMDS执行NMDS推荐分析步骤 # k=2; Number of dimensions. NB., the number of points n should be n > 2*k + 1, # and preferably higher in global non-metric MDS, and still higher in local NMDS. # try, tryamx: Minimum and maximum numbers of random starts in search of stable solution. # After try has been reached, the iteration will stop when two convergent solutions # were found or trymax was reached. # autotransform: Use simple heuristics for possible data transformation of typical # community data (see below). If you do not have community data, # you should probably set autotransform = FALSE. set.seed(1) dune_nmds <- metaMDS(dune, distance="bray", k=2, try=20, trymax=50, autotransform=T, model="global", stress=1, maxit=200, parallel=2, noshare=F) Run 0 stress 0. Run 1 stress 0. Run 2 stress 0. ... Procrustes: rmse 5.e-06 max resid 1.e-05 ... Similar to previous best Run 3 stress 0. ... Procrustes: rmse 6.e-06 max resid 2.061976e-05 ... Similar to previous best Run 4 stress 0. ... New best solution ... Procrustes: rmse 0.0 max resid 0.0 Run 5 stress 0. Run 6 stress 0. Run 7 stress 0. Run 8 stress 0. ... New best solution ... Procrustes: rmse 1.e-06 max resid 5.e-06 ... Similar to previous best Run 9 stress 0. Run 10 stress 0. ... Procrustes: rmse 1.e-06 max resid 4.e-06 ... Similar to previous best Run 11 stress 0. Run 12 stress 0. Run 13 stress 0. Run 14 stress 0. ... New best solution ... Procrustes: rmse 1.088269e-06 max resid 3.e-06 ... Similar to previous best Run 15 stress 0. Run 16 stress 0. ... New best solution ... Procrustes: rmse 1.e-06 max resid 3.e-06 ... Similar to previous best Run 17 stress 0. Run 18 stress 0. Run 19 stress 0. Run 20 stress 0. ... Procrustes: rmse 6.e-06 max resid 1.34356e-05 ... Similar to previous best * Solution reached

获得结果：stress=0.118结果还可以。(虽然指定了autotransform=T，但该套数据并没有触发数据转换)

dune_nmds Call: metaMDS(comm = dune, distance = "bray", k = 2, try = 20, trymax = 50, autotransform = T, noshare = F, model = "global", stress = 1, maxit = 200, parallel = 2) global Multidimensional Scaling using monoMDS Data: dune Distance: bray Dimensions: 2 Stress: 0. Stress type 1, weak ties Two convergent solutions found after 20 tries Scaling: centring, PC rotation, halfchange scaling Species: expanded scores based on 'dune'

如果做了数据转换，会有类似如下的新：

Square root transformation Wisconsin double standardization … global Multidimensional Scaling using monoMDS

Data: wisconsin(sqrt(varespec))

理想情况下，增加的排序距离等于增加的观测距离，所有的点都落在一条单调递增的线上，越一致越好，偏差大则结果差。Stress值就是下图回归线的残差。

stressplot(dune_nmds)

计算每个样品的拟合度，图中圈越小表示拟合度越高。

gof <- goodness(dune_nmds) plot(dune_nmds, type="t", main = "goodness of fit") points(dune_nmds, display="sites", cex=gof*100)

整理数据，绘制NMDS图。

dune_nmds_sample_result <- as.data.frame(scores(dune_nmds)) dune_nmds_sample_result$Sample <- rownames(dune_nmds_sample_result) dune_nmds_sample_result <- cbind(dune_nmds_sample_result, dune.env) dune_nmds_sample_result NMDS1 NMDS2 Sample A1 Moisture Management Use Manure Sample1 -0. -0. Sample1 2.8 1 SF Haypastu 4 Sample2 -0. -0. Sample2 3.5 1 BF Haypastu 2 Sample3 -0.0 -0. Sample3 4.3 2 SF Haypastu 4 Sample4 -0. -0. Sample4 4.2 2 SF Haypastu 4 Sample5 -0. -0.0 Sample5 6.3 1 HF Hayfield 2 Sample6 -0. 0. Sample6 4.3 1 HF Haypastu 2 Sample7 -0. -0.0 Sample7 2.8 1 HF Pasture 3 Sample8 0. -0. Sample8 4.2 5 HF Pasture 3 Sample9 0. -0. Sample9 3.7 4 HF Hayfield 1 Sample10 -0. 0.0 Sample10 3.3 2 BF Hayfield 1 Sample11 -0. 0. Sample11 3.5 1 BF Pasture 1 Sample12 0. -0. Sample12 5.8 4 SF Haypastu 2 Sample13 0. -0. Sample13 6.0 5 SF Haypastu 3 Sample14 0. 0. Sample14 9.3 5 NM Pasture 0 Sample15 0. 0. Sample15 11.5 5 NM Haypastu 0 Sample16 1.0 -0. Sample16 5.7 5 SF Pasture 3 Sample17 -0. 0. Sample17 4.0 2 NM Hayfield 0 Sample18 -0. 0. Sample18 4.6 1 NM Hayfield 0 Sample19 -0.0 1.0 Sample19 3.7 5 NM Hayfield 0 Sample20 1.0 0. Sample20 3.5 5 NM Hayfield 0 # Using the scores function from vegan to extract the species scores # and convert to a data.frame dune_nmds_species_result <- as.data.frame(scores(dune_nmds, "species")) dune_nmds_species_result$Species <- rownames(dune_nmds_species_result) dune_nmds_species_result NMDS1 NMDS2 Species Achimill -0. 0.0 Achimill Agrostol 0. -0. Agrostol Airaprae -0. 1. Airaprae Alopgeni 0. -0. Alopgeni Anthodor -0. 0. Anthodor Bellpere -0. -0. Bellpere Bromhord -0. -0. Bromhord Chenalbu 0. -0. Chenalbu Cirsarve -0. -0. Cirsarve Comapalu 1. 0. Comapalu Eleopalu 1. 0. Eleopalu Elymrepe -0. -0. Elymrepe Empenigr -0.0 1. Empenigr Hyporadi -0. 1. Hyporadi Juncarti 0. -0.0 Juncarti Juncbufo 0. -0. Juncbufo Lolipere -0. -0. Lolipere Planlanc -0. 0. Planlanc Poaprat -0. -0. Poaprat Poatriv -0. -0. Poatriv Ranuflam 1. 0.0 Ranuflam Rumeacet -0. -0. Rumeacet Sagiproc 0. -0. Sagiproc Salirepe 0. 0. Salirepe Scorautu -0. 0. Scorautu Trifprat -0. 0.0 Trifprat Trifrepe -0.0 0.0 Trifrepe Vicilath -0. 0. Vicilath Bracruta 0. 0. Bracruta Callcusp 1. 0. Callcusp

NMDS是一种散点图，图形中每个点通常代表样本，不同颜色/形状等代表样本所属的分组信息或其它关注的样本属性信息。同组内样品点距离远近说明了样本的重复性强弱，组间样本的远近则反应了组间样本在检测变量空间上的差异。通常需要标记stress信息，不标记轴的权重信息。

library(ggplot2) library(ggalt) p1 <- ggplot(data=dune_nmds_sample_result, aes(x=NMDS1,y=NMDS2,fill=Management, colour=Management,group=Management)) + geom_point(size=3) + # add the point markers geom_encircle(alpha=0.1, show.legend = F) + geom_hline(aes(yintercept=0),color="grey", linetype=3) + geom_vline(aes(xintercept=0),color="grey", linetype=3) + labs(title=paste0("Stress: ", round(dune_nmds$stress,3))) + theme_classic() p1

在NMDS二维空间同时绘制物种信息（物种较多时一团乱，可以筛选部分高峰度物种，再筛选dune_nmds_species_result进行绘制。或者也可以直接绘制物种信息。

p1 <- ggplot() + geom_point(data=dune_nmds_sample_result, aes(x=NMDS1,y=NMDS2, colour=Management,group=Management),size=3) + # add the point markers geom_encircle(data=dune_nmds_sample_result, aes(x=NMDS1,y=NMDS2,fill=Management,group=Management), alpha=0.1, show.legend = F) + geom_hline(aes(yintercept=0),color="grey", linetype=3) + geom_vline(aes(xintercept=0),color="grey", linetype=3) + geom_point(data=dune_nmds_species_result, aes(x=NMDS1,y=NMDS2)) + geom_text(data=dune_nmds_species_result, aes(x=NMDS1,y=NMDS2,label=Species)) + labs(title=paste0("Stress: ", round(dune_nmds$stress,3))) + theme_classic() p1

参考

1. NMDS https://archetypalecology.wordpress.com/2018/02/18/non-metric-multidimensional-scaling-nmds-what-how/
2. NMDS https://blog.csdn.net/weixin_/article/details/
3. https://www.cd-genomics.com/microbioseq/non-metric-multidimensional-scaling-nmds-in-microbial-sequencing-data-analysis-introduction-application-and-comparison.html
4. 很好的排序概念解释 http://albertsenlab.org/ampvis2-ordination/
5. https://blog.csdn.net/woodcorpse/article/details/?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2~aggregatepage~first_rank_ecpm_v1~rank_aggregation-4-.pc_agg_rank_aggregation&utm_term=beta%E5%A4%9A%E6%A0%B7%E6%80%A7&spm=1000.2123.3001.4430
6. https://mb3is.megx.net/gustame/dissimilarity-based-methods/nmds
7. https://jkzorz.github.io/2019/06/06/NMDS.html
8. https://mp.weixin..com/s/3NtVtqhEinqtbh2B-qnktg
9. https://mp.weixin..com/s/uWX9CUBNR2ybJ1LxFhKefQ