导读

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")
}