导读
Bray-Curtis NMDS降维分析菌群结构。
一、输入数据
1.1 菌属丰度矩阵
df = data.frame(abs(round(matrix(rnorm(729, 100, 50), 27, 27))))
rownames(df) = paste(rep(paste(rep(c("A", "B", "C"), each=3), rep(c("I","II","III")), sep=""), each=3), 1:3, sep="")
colnames(df) = paste("genus", 1:27, sep="")
1.2 样品分组
df_group = paste(rep(c("A", "B", "C"), each=9), rep(c("I","II","III"), each=3), sep="")
data.frame(rownames(df), df_group)
二、BC NMDS函数
library(vegan)
library(ggplot2)
library(ape)
# nmds
nmds_f = function(input, group)
{
nmds = metaMDS(input, distance="bray", k=2)
nmds_point = data.frame(nmds$points)
nmds_point$Category = group
nmds_point$Group = substring(nmds_point[, 3], 1, 1)
nmds_point$Mark = substring(nmds_point[, 3], 2)
# substring(x, start, end) no end is all
write.csv(nmds_point, file="bray_curtis_nmds.csv")
result = ggplot(nmds_point, aes(x=MDS1, y=MDS2, color = Group)) +
geom_point(aes(shape = Mark), size = 4) +
stat_ellipse(level = 0.95, show.legend = F, aes(fill = Category)) +
theme_classic() + # 经典，只有xy轴
labs(title = paste('Bray-Curtis NMDS Stress =', round(nmds$stress, 4))) +
scale_color_manual(values=c("A" = "#EE6363", "B" = "#1E90FF", "C" = "#00C5CD")) +
scale_shape_manual(values=c(1, 2, 0)) +
guides(color = guide_legend(order=1), shape = guide_legend(order=2))
# guides(color=guide_legend(title = NULL)) +
# theme(legend.title=element_blank()) +
# geom_hline(yintercept = 0, size=0.1, linetype="dashed") +
# geom_vline(xintercept = 0, size=0.1, linetype="dashed") +
# geom_text(aes(label = rownames(nmds_point))) +
# theme(panel.grid=element_blank(), panel.background=element_rect(color='black', fill='transparent')) +
# theme_bw() + # 黑框，无背景
# theme(legend.title = element_text(size = 13, face = "bold"), plot.title = element_text(hjust = 0.5))
ggsave(result, filename="bray_curtis_nmds.pdf")
}
三、分析及结果
nmds_f(df, df_group)
ellipse圈圈没出来，数据是模拟的，没事
3.1 结果文件
3.2 结果表
3.3 结果图
四、分析举例
