基于Bray-Curtis距离的nMDS分组绘图问题求助

我来帮你排查下这个nMDS绘图的问题，看了你的代码，有几个关键问题需要调整：

1. 语法错误直接导致绘图失败 #

你最后一行代码的ylab="axis 2)少了一个闭合的双引号，这会让R直接报错，根本跑不出图。先把这个引号补上，改成：

plot(nMDS$points, col=co[data_2$Treatment], pch=shape[data_2$Treatment], cex=1.2, main="", xlab="axis 1", ylab="axis 2")

2. 分组索引逻辑错误 #

你的data_2$Treatment大概率是字符型的处理名称（比如"Control"、"Treat1"这类），但co[data_2$Treatment]这种用法是错误的——向量索引需要是数字/逻辑值，字符型索引得要求你的颜色/形状向量的名字和Treatment的取值完全对应。

解决方法有两种：

方法一：将Treatment转为因子，用水平索引取值 #

# 先把Treatment转成因子，确保处理顺序符合你的预期
data_2$Treatment <- factor(data_2$Treatment)
# 用因子的水平序号来取颜色和形状
plot(nMDS$points, 
     col=co[as.numeric(data_2$Treatment)], 
     pch=shape[as.numeric(data_2$Treatment)], 
     cex=1.2, main="", xlab="axis 1", ylab="axis 2")

方法二：给颜色/形状向量命名，匹配Treatment取值 #

# 假设你的Treatment取值是"Trt1","Trt2","Trt3","Trt4"，对应命名颜色和形状
co <- c("Trt1"="red", "Trt2"="blue", "Trt3"="green", "Trt4"="black")
shape <- c("Trt1"=1, "Trt2"=2, "Trt3"=3, "Trt4"=4)
# 此时直接用co[data_2$Treatment]就能正确匹配
plot(nMDS$points, col=co[data_2$Treatment], pch=shape[data_2$Treatment], cex=1.2, main="", xlab="axis 1", ylab="axis 2")

3. 更稳妥的vegan专属绘图方式 #

其实vegan自带的绘图函数更适配nMDS结果，你可以先初始化空图，再用points()添加分组点，这样不容易出错，还能顺便加图例：

library("vegan")
data_1 <- Spp_Abund[,4:10]
data_2 <- Spp_Abund[,1:3]
nMDS <- metaMDS(data_1, distance = "bray", k=2)

# 先画空图，不显示任何点
plot(nMDS, type = "n", main="nMDS Plot", xlab="axis 1", ylab="axis 2")
# 处理分组变量
data_2$Treatment <- factor(data_2$Treatment)
co <- c("red", "blue", "green", "black")
shape <- c(1,2,3,4)
# 添加分组点
points(nMDS, col=co[as.numeric(data_2$Treatment)], pch=shape[as.numeric(data_2$Treatment)], cex=1.2)
# 添加图例，方便识别分组
legend("topright", legend=levels(data_2$Treatment), col=co, pch=shape, cex=1)

4. 用ggplot2绘制更美观的nMDS图 #

如果你想要更灵活的样式，推荐用ggplot2来画：

library(ggplot2)
# 提取nMDS坐标并合并处理信息
nmds_df <- as.data.frame(nMDS$points)
nmds_df$Treatment <- data_2$Treatment

# 绘图
ggplot(nmds_df, aes(x=MDS1, y=MDS2, color=Treatment, shape=Treatment)) +
  geom_point(size=3) +
  scale_color_manual(values=c("red", "blue", "green", "black")) +
  scale_shape_manual(values=c(1,2,3,4)) +
  labs(x="axis 1", y="axis 2") +
  theme_bw()

内容的提问来源于stack exchange，提问作者mpearson