基于R语言的微生物群落组成多样性分析——NMDS分析

       NMDS分析，即非度量多维尺度分析（NMDS analysis and plot NMDS），是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法，一般用组间样本的秩次（数据排名rank order）上的差异来定义距离，常用作微生物群落研究的β分析。今天，小编就带大家看一下如何使用vegan包做NMDS分析并进行可视化！

安装、加载R包

rm(list=ls())#clear Global Environmentsetwd('D:\\NMDS')#设置工作路径#安装所需R包install.packages("vegan")install.packages("ggplot2")#加载包library(vegan)#计算距离时需要的包library(ggplot2)#绘图包

加载数据

主要为OTU表格：

#读取数据，一般所需是数据行名为样本名、列名为OTUxxx的数据表otu_raw <- read.table(file="otu.txt",sep="\t",header=T,check.names=FALSE ,row.names=1)#由于排序分析函数所需数据格式原因，需要对数据进行转置otu <- t(otu_raw)

NMDS分析

1、计算bray_curtis距离

otu.distance <- vegdist(otu, method = 'bray')

2、NMDS排序分析

#NMDS排序分析——vegan包中的metaMDS函数df_nmds <- metaMDS(otu.distance, k = 2)#结果查看——关注stress、points及species三个指标summary(df_nmds)

3、适用性检验——基于stress结果

#应力函数值（<=0.2合理）df_nmds_stress <- df_nmds$stressdf_nmds_stress#检查观测值非相似性与排序距离之间的关系——没有点分布在线段较远位置表示该数据可以使用NMDS分析stressplot(df_nmds)

4、提取作图数据

#提取作图数据df_points <- as.data.frame(df_nmds$points)#添加samp1es变量df_points$samples <- row.names(df_points)#修改列名names(df_points)[1:2] <- c('NMDS1', 'NMDS2')head(df_points)

可视化

1、绘制基础散点图

p <- ggplot(df_points,aes(x=NMDS1, y=NMDS2))+#指定数据、X轴、Y轴 geom_point(size=3)+#绘制点图并设定大小 theme_bw()#主题p

2、添加分组数据并进行个性化展示

#读入分组文件group <- read.table("group.txt", sep='\t', header=T)#修改列名colnames(group) <- c("samples","group")#将绘图数据和分组合并df <- merge(df_points,group,by="samples")head(df)#使用ggplot2包绘图color=c("#1597A5","#FFC24B","#FEB3AE")#颜色变量p1<-ggplot(data=df,aes(x=NMDS1,y=NMDS2))+#指定数据、X轴、Y轴，颜色 theme_bw()+#主题设置 geom_point(aes(color = group), shape = 19, size=3)+#绘制点图并设定大小 theme(panel.grid = element_blank())+ geom_vline(xintercept = 0,lty="dashed", size = 1, color = 'grey50')+ geom_hline(yintercept = 0,lty="dashed", size = 1, color = 'grey50')+#图中虚线 geom_text(aes(label=samples, y=NMDS2+0.03,x=NMDS1+0.03, vjust=0, color = group),size=3.5, show.legend = F)+#添加数据点的标签 stat_ellipse(data=df, geom = "polygon",level=0.95, linetype = 2,size=0.5, aes(fill=group), alpha=0.2)+ scale_color_manual(values = color) +#点的颜色设置 scale_fill_manual(values = color)+#椭圆颜色 theme(axis.title.x=element_text(size=12),#修改X轴标题文本 axis.title.y=element_text(size=12,angle=90),#修改y轴标题文本 axis.text.y=element_text(size=10),#修改x轴刻度标签文本 axis.text.x=element_text(size=10),#修改y轴刻度标签文本 panel.grid=element_blank())+#隐藏网格线 ggtitle(paste('Stress=',round(df_nmds_stress, 3)))#添加应力函数值p1

3、也可通过添加边际箱线图展示组间差异性

#加载包，对组间进行统计检验以及组合图的拼接library(ggpubr)library(ggsignif)# 绘制y轴为PC2值的分组箱线图p2 <- ggplot(df,aes(x=group,y=NMDS2))+ stat_boxplot(geom = "errorbar", width=0.1,size=0.5)+ geom_boxplot(aes(fill=group), outlier.colour="white",size=0.5)+  theme(panel.background =element_blank(),  axis.line=element_line(color = "white"), axis.text.y = element_blank(), panel.grid = element_blank(), axis.ticks = element_blank(), legend.position = 'none')+ xlab("") + ylab("")+ scale_fill_manual(values=c("#1597A5","#FFC24B","#FEB3AE"))+ geom_signif(comparisons = list(c("A","B"), c("A","C"), c("B","C")),              map_signif_level = T,               test = t.test, y_position = c(0.3,0.4,0.35), tip_length = c(c(0,0), c(0,0), c(0,0)), size=0.8,color="black")p2# 绘制y轴为PC1值的分组箱线图p3 <- ggplot(df,aes(x=group,y=NMDS1))+ stat_boxplot(geom = "errorbar", width=0.1,size=0.5)+ coord_flip()+  geom_boxplot(aes(fill=group),  outlier.colour="white",size=0.5)+ theme(panel.background =element_blank(), axis.line=element_line(color = "white"), axis.text.x = element_blank(), panel.grid = element_blank(), axis.ticks = element_blank(), legend.position = 'none')+ xlab("") + ylab("")+ scale_fill_manual(values=c("#1597A5","#FFC24B","#FEB3AE"))+ geom_signif(comparisons = list(c("A","B"), c("A","C"), c("B","C")),              map_signif_level = T, test = t.test, y_position = c(0.48,0.62,0.55), tip_length = c(c(0,0), c(0,0), c(0,0)), size=0.8,color="black")p3
# ggpubr::ggarrange()函数对图进行拼接ggarrange(p3, NULL, p1, p2, widths = c(5,2), heights = c(2,4), align = "hv")

4、通过AI进行微调：

源码及作图数据可在后台回复NMDS获取！

我就知道你“在看”