R:Bray-Curtis PCoA

导读

Bary-Curtis PCoA降维分析菌群结构。

一、输入数据

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 PCoA函数

library(vegan)
library(ggplot2)
library(ape)

pcoa_f = function(input, group)
{
    pcoa = pcoa(vegdist(input, method = "bray"))

    pcoa_point = data.frame(pcoa$vectors[, c(1, 2)], Category = group)
    pcoa_point$Group = substring(pcoa_point[, 3], 1, 1)
    pcoa_point$Mark = substring(pcoa_point[, 3], 2)
    write.csv(pcoa_point, file="bray_curtis_pcoa.csv")
    
    xylab = paste(c("Axis.1 [", "Axis.2 ["), round(pcoa$values$Relative_eig[1:2]*100, 2), "%]", sep = "")

    result = ggplot(pcoa_point, aes(x=Axis.1, y=Axis.2, color = Group)) +
      geom_point(aes(shape = Mark), size = 4) + 
      stat_ellipse(level = 0.95, show.legend = F, aes(fill = Category)) +
      theme_classic() +
      labs(x=xylab[1], y=xylab[2], title = paste('Bray-Curtis PCoA')) +
      scale_color_manual(values=c("A" = "#EE6363", "B" = "#1E90FF", "C" = "#00C5CD")) +
      guides(color = guide_legend(order=1), shape = guide_legend(order=2))
    ggsave(result, filename="bray_curtis_pcoa.pdf")
}

三、分析及结果

pcoa_f(df, df_group)


ellipse出错,圈圈没有出来,模拟数据效果不好

3.1 结果文件

3.2 结果表

3.3 结果图

四、实际分析举例图

一组bray curtis pcoa参数

ggplot(pcoa_point, aes(x=Axis.1, y=Axis.2, fill=source)) +
  geom_point(pch=21, size=7, color="white") + 
  stat_ellipse(level = 0.95, show.legend = F, 
               aes(color=source), size = 1) +
  theme_classic() +
  labs(x=xylab[1], y=xylab[2], 
       title = paste('Bray-Curtis PCoA'),
       fill = "Source") +
  theme(title = element_text(size = 15, face="bold")) +
  theme(legend.text=element_text(size=15)) +
  theme(legend.title=element_text(face='bold', size=20)) +
  theme(axis.title = element_text(size = 20)) +
  theme(axis.text = element_text(size = 20),
        axis.line = element_line(size = 1.2),
        axis.ticks = element_line(size = 1.2)) +
  theme(plot.title = element_text(face="bold", size=16)) +
  scale_fill_manual(
    values = c("Colon"="#00BA38",
               "Stool"="#F8766D",
               "Ileum"="#619CFF")) +
  scale_color_manual(
    values = c("Colon"="#00BA38",
               "Stool"="#F8766D",
               "Ileum"="#619CFF"))

更多参数

library("ggalt")
geom_encircle(aes(group = Group, colour = Group))  # 外点连线

更多阅读:
三文读懂PCA和PCoA
数量生态学笔记||非约束排序||PCoA
stat_ellipse参数
限制性/非限制性排序分析
画一个带统计检验的PCoA分析结果
Beta多样性PCoA和NMDS排序 三维展示技巧
带统计学的PCoA完美解决打样本量多组数据不好区分的问题!!
排序图(beta多样性)添加各种辅助线-无非是让不明显的规律看得到
RDA_环境因子群落结构统计检验_可视化 添加环境因子箭头

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容