非度量多维尺度（NMDS）分析及R绘图

先上图

image.png
注：填充颜色由AI处理。
（样本和分组文件如下）（这里只是展示数据形式并不是原数据（分组名与代码中不一样））
image.png
image.png

vegan分析数据

行名样本，列名物种，内容丰度。

#加载R包
library(yyplot)
library(ggbiplot)
library(vegan)
library(ggplot2)
#导入样本和分组文件
sample <- read.table("sample.txt",sep = "\t",header = T)
group <- read.table("group1.txt",sep = "\t",header = T,row.names = 1)
group <- group[match(rownames(sample),rownames(group)),] #匹配行名，很重要
#nmds分析
nmds1 <- metaMDS(sample, distance = 'bray', k = 2)
summary(nmds1)
#提取数据
nmds1.stress <- nmds1$stress
nmds1.point <- data.frame(nmds1$point)
nmds1.species <- data.frame(nmds1$species)
sample_site <- nmds1.point[1:2]
sample_site$names <- rownames(sample_site)
colnames(sample_site)[1:2] <- c('NMDS1', 'NMDS2')
#合并分组数据
sample_site <- cbind(sample_site,group)

绘图

#NMDS图绘制
nmds_plot <- ggplot() +
  geom_point(data = sample_site, aes(NMDS1, NMDS2,color = site , shape = depth), size = 5, alpha = 0.8) + #可在这里修改点的透明度、大小
  #scale_shape_manual(values = c(17, 16)) + #可在这里修改点的形状
  #scale_color_manual(values = c('red', 'blue')) + #可在这里修改点的颜色
  geom_ord_ellipse(aes(sample_site$NMDS1,sample_site$NMDS2,group= sample_site$site2), ##添加0.8置信椭圆
                   ellipse_pro = 0.8,linetype=2,size=0.7,color='firebrick')+
  geom_ord_ellipse(aes(sample_site$NMDS1,sample_site$NMDS2,color= sample_site$site2,group= sample_site$site2),
                   ellipse_pro = 0.9,linetype=3,size=1)+ ##添加0.9置信椭圆
  theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) + #去掉背景  
  theme(legend.key = element_rect(fill = 'transparent'), legend.title = element_blank()) + #去掉图例标题及标签背景
  labs(x = 'NMDS axis1', y = 'NMDS axis2', title = paste('Stress =', round(nmds1$stress, 4))) +
  theme(plot.title = element_text(hjust = 0.5))+  #标题居中
  theme(panel.background = element_blank(),axis.line = element_line(color = "black"))#去上右边框
  
nmds_plot

注意！！geom_ord_ellipse()函数中的参数需要data$列名 而不是直接用列名

geom_ord_ellipse()函数参数中的fill好像只能填充一种颜色。默认无填充。

ggplot2自带的stat_ellipse() 同样0.8参数情况下会绘制不同大小的圈。

用法:

p+stat_ellipse(aes(fill=group),geom="polygon",level=0.8,alpha=0.1,show.lengend=F)
#或者
p+stat_ellipse(aes(color=group),level=0.8,linetype=2,show.lengend=F)

以下是NMDS结果的粗略评估

#stress记录了NMDS排序分析的应力函数值，points记录了各样本的排序坐标，
#species记录了各物种（OTU）的排序坐标。
#而在NMDS排序分析中，尽可能选择较低的应力函数值。一般情况下，应力函数值的值不要大于0.2。
write.csv(nmds1.point, 'nmds.sample.csv')
#简单绘图评估
nmds_plot <- nmds1
nmds_plot$species <- {nmds_plot$species}[1:10, ]
stressplot(nmds_plot, main = 'Shepard 图')
gof <- goodness(nmds_plot)
plot(nmds_plot,type = 't', main = '拟合度')
points(nmds_plot, display = 'sites', cex = gof * 200, col = 'red')

关键词：nmds，ggplot2，yyplot，gglayer，ggbiplot，github，geom_ord_ellipse
参考网站：
https://guangchuangyu.github.io/cn/2018/01/geom-ord-ellipse/
https://github.com/GuangchuangYu/gglayer#geom_ord_ellipse
https://mp.weixin.qq.com/s/4vw6chPvrz_tLmLXHQ7ThQ
https://mp.weixin.qq.com/s/_XH6u2hPv_JJi_MbGbeoHA
https://blog.csdn.net/qq_38854576/article/details/83024468/