使用R语言,进行NMDS统计和绘图
#!/usr/bin/env Rscript
#install.packages("vegan") #安装vegan包
library(vegan)
sample <- read.table(file.choose(),sep = "\t",row.names = 1)#第一行为列名
group <- read.table(file.choose(),sep = "\t",header = T,row.names = 1)#其中header=T表示将文件中第一行设为列名字。row.names= 1表示第一列设为行名
group1 <- group[match(rownames(sample),rownames(group)),] #匹配行名,很重要group = group[rownames(sample),]
#nmds分析
nmds1 <- metaMDS(sample, distance = 'bray', k = 3)
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,group1)
#给分组排序
sample_site$Depth <- factor(sample_site$Depth,levels =c('0m','25m','50m','117m','150m', '200m',
ordered = FALSE ))
sample_site$Particlesize <- factor(sample_site$Particlesize, levels = c("0.2-0.6μm","0.6-1.2μm","1.2-2μm","2-20μm",">20μm",
ordered = FALSE))
View(sample_site)
#NMDS图绘制
windowsFonts(TNM = windowsFont("Times New Roman"))#设置字体
nmds_plot <- ggplot() +
geom_point(data = sample_site, aes(NMDS1, NMDS2,color = Particlesize, shape = Depth), size = 5, alpha = 0.8) + #可在这里修改点的透明度、大小
scale_shape_manual(values = c(15,17,18,19,20,79)) + #可在这里修改点的形状
scale_color_manual(values = c("#99CCFF","#FF9999","#FFCC33","#99CC99",
"#996699" )) + #可在这里修改点的颜色
scale_fill_manual(values = c("black" ),aesthetics = c("fill")) +#填充颜色
#geom_ord_ellipse(aes(sample_site$NMDS1,sample_site$NMDS2,group= sample_site$Particlesize),ellipse_pro = 0.8,linetype=2,size=0.7,color='firebrick')+##添加0.8置信椭圆
#geom_ord_ellipse(aes(sample_site$NMDS1,sample_site$NMDS2,group= sample_site$Particlesize,color= sample_site$Particlesize),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.text =element_text(family="TNM"), legend.title = element_blank()) + #去掉图例标题及标签背景
labs(x = 'Dimension1', y = 'Dimension2', title = paste('Stress =', round(nmds1$stress, 2))) +
theme(plot.title = element_text(hjust = 0.8,family="TNM",face = "bold"))+ #标题
theme(axis.text.x =element_text(family="TNM",face = "bold"),axis.text.y =element_text(family="TNM",face = "bold"),
legend.text = element_text(family="TNM"),panel.grid.minor = element_blank())+
theme(axis.title.y = element_text(family="TNM",face = "bold"),axis.title.x = element_text(family="TNM",face = "bold"))
#theme(panel.background = element_blank(),axis.line = element_line(color = "black"))#去上右边框
nmds_plot