R包vegan的非度量多维标度（NMDS）分析

非度量多维标度（Non-metric Multidimensional Scaling，NMDS）分析，是主坐标分析（PCoA）的非度量替代方法。NMDS基于对象之间给定的距离测度（可以为任意类型的距离），尝试在预先设定数量的排序轴中排序对象，排序目的不是最大限度地保留对象之间给定的原始距离，而是反映对象之间的顺序关系。

NMDS的基本概念，可参考前文。

本篇以某16S扩增子测序所得细菌群落数据为例，简介R包vegan的NMDS过程。

示例文件、R脚本等，已上传至百度盘：

https://pan.baidu.com/s/1g_EjuWxtZokLHftlx_YRLQ

物种数据的NMDS分析

网盘文件“otu_table.txt”为通过16S测序所得的细菌群落丰度表格，OTU水平。

接下来期望通过NMDS，评估3组群落的物种组成差异水平。

NMDS执行

vegan中，metaMDS()函数用于执行NMDS。metaMDS()随机构建对象的初始结构，并不断尝试以获得最佳的NMDS排序结果。

对于输入metaMDS()的数据，可以为两种形式：

（1）已经计算好的距离矩阵，可以为任意距离测度；

（2）样方-物种丰度矩阵，即原初的物种丰度表，同时在metaMDS()函数中通过distance参数指定需要的距离类型。

注意：两种模式看起来似乎是相同的，但实际上计算后，结果并不一致。详见下文。

##第 1 种模式，输入距离矩阵排序
#读取 OTU 丰度表
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))
 
#根据物种组成计算样方距离，如 Bray-curtis 距离，详情 ?vegdist
bray_dis <- vegdist(otu, method = 'bray')      #结果以 dist 数据类型存储
 
#输出距离矩阵
#write.table(as.matrix(bray_dis), 'bray_distance.txt', sep = '\t', col.names = NA, quote = FALSE)
 
#或者读取已经准备好的距离矩阵文件，如 Bray-curtis 距离，排序
dis <- read.delim('bray_distance.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
bray_dis <- as.dist(dis)   #转为 dist 数据类型
 
#NMDS 排序，定义 2 个维度，详情 ?metaMDS
nmds_dis <- metaMDS(bray_dis, k = 2)
 
##第 2 种模式，直接输入 OTU 丰度表，在函数中指定距离参数排序
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))
 
#NMDS 排序，定义 2 个维度，详情 ?metaMDS
nmds_otu <- metaMDS(otu, distance = 'bray', k = 2)

查看主要结果内容

将主要的信息，如应力函数值、样方得分等提取出来。

这时通过比较结果，就能看到两种模式存在区别了。别问我这是为什么，我不知道……问原作者去。

##第 1 种模式
#应力函数值，一般不大于 0.2 为合理
nmds_dis$stress
 
#样方得分
nmds_dis_site <- data.frame(nmds_dis$points)
#write.table(nmds_dis_site, 'nmds_dis_site.txt', sep = '\t', col.names = NA, quote = FALSE)
 
#由于物种变量在距离矩阵的计算过程中丢失，因此若想补充物种变量
#物种变量可通过丰度加权平方被动添加至排序图中，详情 ?wascores
nmds_dis_species <- wascores(nmds_dis$points, otu)
#write.table(nmds_dis_species, 'nmds_dis_species.txt', sep = '\t', col.names = NA, quote = FALSE)
 
##第 2 种模式
#应力函数值，一般不大于 0.2 为合理
nmds_otu$stress
 
#样方得分
nmds_otu_site <- data.frame(nmds_otu$points)
#write.table(nmds_otu_site, 'nmds_otu_site.txt', sep = '\t', col.names = NA, quote = FALSE)
 
#物种得分，这种模式下可直接计算出物种得分，具体怎么算出来的，问作者吧……
nmds_otu_species <- data.frame(nmds_otu$species)
#write.table(nmds_otu_species, 'nmds_otu_species.txt', sep = '\t', col.names = NA, quote = FALSE)

别问我哪种模式更好，根据本菜鸡的经验：

很多情况中，第2种模式一般更能区分组间差异，且stress值也小。但仅能指定特定的距离类型（metaMDS()函数的distance参数），对于其它vegan中没有的距离类型，应用起来倒不是不行，但必须修改原函数，就比较复杂了。

然而就算法而言，第1种模式更符合NMDS定义，而且距离测度适用广泛，只要提供现有的距离矩阵（任意来源，不受metaMDS()函数的distance限制）即可。

所以......还请大家看着自己琢磨吧.......

NMDS排序图

然后接上文，再通过排序图比较一下两种模式的差别。

#ordiplot() 作图展示，比较两种方法的不同
par(mfrow = c(2, 2))

#输入距离矩阵直接排序，仅样方
ordiplot(nmds_dis, type = 'none', main = paste('第1种模式，仅样方，Stress =', round(nmds_dis$stress, 4)))
points(nmds_dis, pch = 19, cex = 0.7, col = c(rep('red', 12), rep('orange', 12), rep('green3', 12)))

#输入距离矩阵直接排序，样方 + 被动物种投影
ordiplot(nmds_dis, type = 'none', main = paste('第1种模式，样方+物种，Stress =', round(nmds_dis$stress, 4)))
points(nmds_dis_species, pch = 3, cex = 0.5, col = 'gray')
points(nmds_dis, pch = 19, cex = 0.7, col = c(rep('red', 12), rep('orange', 12), rep('green3', 12)))

#输入原始物种丰度表并指定距离类型的排序，仅样方
ordiplot(nmds_otu, type = 'none', display = 'site', main = paste('第2种模式，仅样方，Stress =', round(nmds_otu$stress, 4)))
points(nmds_otu, display = 'site', pch = 19, cex = 0.7, col = c(rep('red', 12), rep('orange', 12), rep('green3', 12)))

#输入原始物种丰度表并指定距离类型的排序，样方 +物种
ordiplot(nmds_otu, type = 'none', main = paste('第2种模式，样方+物种，Stress =', round(nmds_otu$stress, 4)))
points(nmds_otu, display = 'sp', pch = 3, cex = 0.5, col = 'gray')
points(nmds_otu, display = 'site', pch = 19, cex = 0.7, col = c(rep('red', 12), rep('orange', 12), rep('green3', 12)))
#怎么能少得了 ggplot2 的份？再基于导出的坐标作个图吧，按第 1 种模式的结果来
library(ggplot2)

#物种太多，就选代表性的展示，比方说 top10 丰度物种
abundance <- apply(otu, 2, sum)
abundance_top10 <- names(abundance[order(abundance, decreasing = TRUE)][1:10])

species_top10 <- data.frame(nmds_dis_species[abundance_top10,1:2])
species_top10$name <- rownames(species_top10)

#添加分组信息
nmds_dis_site$name <- rownames(nmds_dis_site)
nmds_dis_site$group <- c(rep('A', 12), rep('B', 12), rep('C', 12))

p <- ggplot(data = nmds_dis_site, aes(MDS1, MDS2)) +
geom_point(aes(color = group)) +
stat_ellipse(aes(fill = group), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) + #添加置信椭圆，注意不是聚类
scale_color_manual(values = c('red3', 'orange3', 'green3')) +
scale_fill_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid.major = element_line(color = 'gray', size = 0.2), panel.background = element_rect(color = 'black', fill = 'transparent'),
plot.title = element_text(hjust = 0.5), legend.position = 'none') +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_text(data = species_top10, aes(label = name), color = 'blue', size = 3) + #添加 top10 丰度物种标签
labs(x = 'NMDS1', y = 'NMDS1', title = 'NMDS（带top10丰度物种投影）') +
annotate('text', label = paste('Stress =', round(nmds_dis$stress, 4)), x = 0.35, y = 0.3, size = 4, colour = 'black') + #标注应力函数值
annotate('text', label = 'A', x = -0.27, y = 0.08, size = 5, colour = 'red3') +
annotate('text', label = 'B', x = -0.03, y = 0.11, size = 5, colour = 'orange3') +
annotate('text', label = 'C', x = 0.13, y = -0.15, size = 5, colour = 'green3')

p

#ggsave('NMDS_dis.pdf', p, width = 5.5, height = 5.5)
ggsave('NMDS_dis.png', p, width = 5.5, height = 5.5)

NMDS评估

如前所述，NMDS不是基于特征向量提取的排序方法，因此与PCA、CA、PCoA等不同，NMDS轴不存在解释量一说（文献中有时看到NMDS轴标注了解释量，鬼知道怎么算的，问作者去）。通常情况下根据应力函数值来判断排序模型的合理性，值不大于0.2是合理的。

此外，还可通过比较NMDS排序图内对象的（欧式）距离与原始距离矩阵中的值，进行线性或非线性回归的R²来评估NMDS的拟合度。

#第 1 种模式，NMDS 评估
par(mfrow = c(1, 2))
stressplot(nmds_dis, main = '第1种模式，Shepard 图')
gof <- goodness(nmds_dis)
plot(nmds_dis,type = 'text', main = '第1种模式，拟合度')
points(nmds_dis, display = 'site', cex = gof * 100, col = 'red')

#第 2 种模式，NMDS 评估
par(mfrow = c(1, 2))
stressplot(nmds_otu, main = '第2种模式，Shepard 图')
gof <- goodness(nmds_otu)
plot(nmds_otu, type = 'text', display = 'sites', main = '第2种模式，拟合度')
points(nmds_otu, display = 'site', cex = gof * 100, col = 'red')

#注：拟合R2越大越合理；气泡图越小越合理

友情链接