搜索

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

生信宝典

共 913字，需浏览 2分钟

 · 2022-04-26

多维排列 (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/538918v1.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.1192678 
## Run 1 stress 0.1886532 
## Run 2 stress 0.1192678 
## ... Procrustes: rmse 5.822837e-06  max resid 1.845818e-05 
## ... Similar to previous best
## Run 3 stress 0.1192678 
## ... Procrustes: rmse 6.697234e-06  max resid 2.061976e-05 
## ... Similar to previous best
## Run 4 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 0.02027078  max resid 0.06496407 
## Run 5 stress 0.1192679 
## Run 6 stress 0.1939202 
## Run 7 stress 0.1808911 
## Run 8 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 1.651369e-06  max resid 5.655239e-06 
## ... Similar to previous best
## Run 9 stress 0.1192678 
## Run 10 stress 0.1183186 
## ... Procrustes: rmse 1.505928e-06  max resid 4.480433e-06 
## ... Similar to previous best
## Run 11 stress 0.1192679 
## Run 12 stress 0.1192678 
## Run 13 stress 0.1886532 
## Run 14 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 1.088269e-06  max resid 3.328672e-06 
## ... Similar to previous best
## Run 15 stress 0.1192678 
## Run 16 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 1.797827e-06  max resid 3.862713e-06 
## ... Similar to previous best
## Run 17 stress 0.2075713 
## Run 18 stress 0.1192678 
## Run 19 stress 0.1192679 
## Run 20 stress 0.1183186 
## ... Procrustes: rmse 6.241246e-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.1183186 
## 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.84054551 -0.71582776  Sample1  2.8        1         SF Haypastu      4
## Sample2  -0.50485425 -0.40893664  Sample2  3.5        1         BF Haypastu      2
## Sample3  -0.08267021 -0.43667531  Sample3  4.3        2         SF Haypastu      4
## Sample4  -0.11562467 -0.52223859  Sample4  4.2        2         SF Haypastu      4
## Sample5  -0.62654363 -0.08669513  Sample5  6.3        1         HF Hayfield      2
## Sample6  -0.54269773  0.11315559  Sample6  4.3        1         HF Haypastu      2
## Sample7  -0.54030010 -0.05820584  Sample7  2.8        1         HF  Pasture      3
## Sample8   0.28115518 -0.16683738  Sample8  4.2        5         HF  Pasture      3
## Sample9   0.11057585 -0.44257696  Sample9  3.7        4         HF Hayfield      1
## Sample10 -0.51697230  0.02738090 Sample10  3.3        2         BF Hayfield      1
## Sample11 -0.33831007  0.35081631 Sample11  3.5        1         BF  Pasture      1
## Sample12  0.44246800 -0.36410934 Sample12  5.8        4         SF Haypastu      2
## Sample13  0.41863051 -0.58334738 Sample13  6.0        5         SF Haypastu      3
## Sample14  0.94331987  0.47606816 Sample14  9.3        5         NM  Pasture      0
## Sample15  0.89599692  0.22235599 Sample15 11.5        5         NM Haypastu      0
## Sample16  1.08109124 -0.17966386 Sample16  5.7        5         SF  Pasture      3
## Sample17 -0.85990338  0.98711836 Sample17  4.0        2         NM Hayfield      0
## Sample18 -0.17720324  0.52341328 Sample18  4.6        1         NM Hayfield      0
## Sample19 -0.07001774  1.01213863 Sample19  3.7        5         NM Hayfield      0
## Sample20  1.04240524  0.25266697 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.82280076  0.04328704 Achimill
## Agrostol  0.71096003 -0.28924414 Agrostol
## Airaprae -0.52817837  1.67987383 Airaprae
## Alopgeni  0.39096747 -0.58596256 Alopgeni
## Anthodor -0.72021737  0.65914560 Anthodor
## Bellpere -0.47838087 -0.24446722 Bellpere
## Bromhord -0.61896693 -0.33475908 Bromhord
## Chenalbu  0.59187098 -0.92197922 Chenalbu
## Cirsarve -0.15184056 -0.82169646 Cirsarve
## Comapalu  1.28932838  0.60840015 Comapalu
## Eleopalu  1.24504861  0.16149275 Eleopalu
## Elymrepe -0.42014574 -0.68023468 Elymrepe
## Empenigr -0.08835329  1.69629767 Empenigr
## Hyporadi -0.41569548  1.44600334 Hyporadi
## Juncarti  0.91146010 -0.08309672 Juncarti
## Juncbufo  0.26476745 -0.60761447 Juncbufo
## Lolipere -0.51199146 -0.24807077 Lolipere
## Planlanc -0.70645212  0.32063694 Planlanc
## Poaprat  -0.38844035 -0.25091394  Poaprat
## Poatriv  -0.15906762 -0.47836751  Poatriv
## Ranuflam  1.14363918  0.09953894 Ranuflam
## Rumeacet -0.52479268 -0.10530738 Rumeacet
## Sagiproc  0.14315579 -0.18745689 Sagiproc
## Salirepe  0.57484082  0.91104875 Salirepe
## Scorautu -0.13956707  0.25000522 Scorautu
## Trifprat -0.77154453  0.08564591 Trifprat
## Trifrepe -0.07526697  0.04516674 Trifrepe
## Vicilath -0.46793995  0.54914931 Vicilath
## Bracruta  0.15071778  0.18979787 Bracruta
## Callcusp  1.42118155  0.38376646 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_53819139/article/details/114133818

  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/107032032?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2~aggregatepage~first_rank_ecpm_v1~rank_aggregation-4-107032032.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.qq.com/s/3NtVtqhEinqtbh2B-qnktg

  9. https://mp.weixin.qq.com/s/uWX9CUBNR2ybJ1LxFhKefQ

往期精品(点击图片直达文字对应教程)

机器学习

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集


浏览 39
点赞
评论
收藏
分享

手机扫一扫分享

举报
评论
图片
表情
推荐
为啥我敢说Python是数据分析界的扛把子语言？
首先声明下这篇文字不是卖课的，也不是无脑吹Python，咱只讲事实，认认真真讨论下Python是不是数据分析领域最好的语言。因为我在知乎上看到非常多人在问这个问题，想必大家是关心的。我的观点是，目前来看所有编程语言里，做数据分析Python是最好的选择，没有之一。列举几个事实：1、Kaggle、天池
Python大数据分析
9
文献分享 | 体外冲击波+本体感觉训练治疗「运动创伤性踝关节炎」患者的效果分析
■  近年来，国内患有运动创伤性踝关节炎患者的数量在逐年增加，患者主要为青壮年。在日常生活与工作中，不健康运动方式或者运动过度，都会对踝关节健康造成危害。■  运动创伤性踝关节炎在临床上还被称为足球踝，该病症大多数是运动员踝关节出现急性损伤后，过早负重，韧带因为尚未修复
乐普医疗AI
0
AKShare-股票数据-赚钱效应分析
作者寄语本次更新股票数据-赚钱效应分析接口。本次主要修复该接口，目前该接口可以获取 上涨，涨停，真实涨停，st st*涨停，下跌，跌停 等字段的数据。欢迎加入专注于财经数据和量化投资的【数据科学实战】社区，加入方式参见文末！更新接口"stock_market_activity_legu"  
数据科学实战
275
马上要结婚了，这次又因为买房吵起来了，目前手里有80万，想买长沙全款就能搞定，女朋友非要让我贷款50万去成都买，怎么办才好？
婚姻是一座围城，城外的人想进去，城里的人想出来。这是钱锺书在1947年发表的经典名句，穿越时代历久弥新。有人说恋爱时的核心利益是欣赏，是情感交换，彼此凭感觉就够了，互相让双方开心就可以。而婚姻时的核心利益是成全，是价值交换，一旦结婚，就看双方能否长期为彼此提供“价值”。有位马上要结婚的小哥，如今即将
TTTEED
3
基于源码的 Webpack 结构分析
分享背景即使目前优秀的构建工具层出不穷，Webpack 还是保持着其在现代前端开发工具链中不可替代的地位。这主要得益于其优秀的灵活性以及强大的生态系统。然而，随着版本更替，Webpack 的功能越来越庞大，整体的代码量日渐夸张，大大提高了学习难度。与此同时，大多数人对 Webpack 的使用都停留在
全栈前端精选
1
一个超强的Python库！HTTP请求性能分析工具推荐：httpstat！
什么是Python httpstat？httpstat是一个基于命令行的工具，用于在终端中展示HTTP请求的详细统计信息。它以可视化和易读的方式显示了HTTP请求的各个阶段的性能数据，如DNS解析、TCP连接、TLS握手、发送请求、服务器处理、接收响应等。使用httpstat可以帮助开发人员更好地理
测试开发技术
10
100种分析思维模型之：本质思维
你好，我是林骥。北京有个出租车司机，原来住在四合院，房改之后，他就把四合院卖了，买了一辆桑塔纳，当时感觉风光无限，现如今则是后悔不已。假如这个出租车司机当时具有本质思维，也许就不会做出卖房买车的决定。下面介绍 100 种分析思维模型的第 97 种：本质思维，它能帮助我们更好地看清事物的本质，进而做出
林骥
18
英伟达NVLink技术深度分析
NVIDIA NVLink已成为高性能计算（HPC）和人工智能领域的关键技术。本文深入探讨了NVLink的复杂性，并介绍NVSwitch芯片、NVLink服务器和NVLink交换机，以揭示其在不断演进的先进计算领域中的重要性。关于CPU、服务器和存储详细技术，请参考“服务器基础知识全解(终极版)”，
架构师技术联盟
10
LIST | 文本分析代码列表
个人感觉 textdata.cn 文本分析代码案例都集中在这里了，我将内容按大类分成- Python语法- 数据采集- 数据处理&Pandas- 文本分析- 数据标注&机器学习- 可视化- R语言-&nb
大邓和他的Python
254
Python实证指标构建与文本分析
概览为何要学Python？在科学研究中，数据的获取及分析是最重要的也是最棘手的两个环节！在前大数据时代，一般使用实验法、调查问卷、访谈或者二手数据等方式，将数据整理为结构化的表格数据，之后再使用各种计量分析方法，对这些表格数据进行分析。但大数据时代，网络数据成为各方学者亟待挖掘的潜在宝藏，大量商业信
大邓和他的Python
457
点赞
评论
收藏
分享

手机扫一扫分享

举报