先上图
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/
网友评论