基于MicrobiotaProcess包进行PCoA分析！

 

点击上方

“科研后花园”

关注我们

具体过程如下：

1、加载所需R包：

rm(list = ls())setwd("D:/桌面/PCoA")#加载R包library(MicrobiotaProcess) # A comprehensive R package for managing and analyzing microbiome and other ecological data within the tidy framework # A comprehensive R package for managing and analyzing microbiome and other ecological data within the tidy frameworklibrary(microeco) # Microbial Community Ecology Data Analysislibrary(dplyr) # A Grammar of Data Manipulationlibrary(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics # Create Elegant Data Visualisations Using the Grammar of Graphics # Create Elegant Data Visualisations Using the Grammar of Graphics # Create Elegant Data Visualisations Using the Grammar of Graphics # Create Elegant Data Visualisations Using the Grammar of Graphics # Create Elegant Data Visualisations Using the Grammar of Graphicslibrary(phyloseq) # Handling and analysis of high-throughput microbiome census data

２、加载数据（这里以microeco自带数据集为例）：

data(sample_info_16S)data(otu_table_16S)data(taxonomy_table_16S)data(phylo_tree_16S)sample_info_16S <- as.data.frame(sample_info_16S)otu_table_16S <- as.data.frame(otu_table_16S)taxonomy_table_16S <- as.data.frame(taxonomy_table_16S)phylo_tree_16S <- as.data.frame(phylo_tree_16S)

３、数据清洗：

#统一分类信息taxonomy_table_16S %<>% tidy_taxonomy## 构造microtabledf <- microtable$new(sample_table = sample_info_16S, otu_table = otu_table_16S, tax_table = taxonomy_table_16S, phylo_tree = phylo_tree_16S, auto_tidy = F)##去除不属于非古菌和细菌的OTU# df$tax_table %<>% subset(Kingdom == "k__Archaea" | Kingdom == "k__Bacteria")df$tax_table %<>% .[grepl("Bacteria|Archaea", .$Kingdom), ]##去除“线粒体”和“叶绿体”污染df$filter_pollution(taxa = c("mitochondria", "chloroplast"))##统一各数据的样本和OTU信息df$tidy_dataset()##检查序列号df$sample_sums() %>% range##重采样以减少测序深度对多样性测量的影响，使每个样本的序列号相等。df$rarefy_samples(sample.size = 10000)df$sample_sums() %>% range

4、利用phyloseq包重新构造将数据转换为可用于MicrobiotaProcess包分析的数据格式：

ps <- phyloseq(sample_data(df$sample_table), otu_table(as.matrix(df$otu_table), taxa_are_rows=TRUE), tax_table(as.matrix(df$tax_table)))#转换数据格式df <- ps %>% as.MPSE()

5、PCoA分析：

###样本或组之间的距离计算##数据标准化df %<>% mp_decostand(.abundance=Abundance, method = "hellinger" #数据标准化方法, 可选择'total', 'max', 'frequency', 'normalize', 'range', 'rank', 'rrank', 'standardize' 'pa', 'chi.square', 'hellinger' and 'log', 具体可见??decostand               )
## 计算样本之间的距离，详细见??mp_cal_distdf %<>% mp_cal_dist(.abundance=hellinger, distmethod="bray")
###PCoA分析,详见??mp_cal_pcoadf %<>%   mp_cal_pcoa(.abundance=hellinger, distmethod="bray", action="add")
# 用adonis或anosim方法来检验各分组对群体的差异性是否显著df %<>% mp_adonis(.abundance=hellinger, .formula=~Group, distmethod="bray", permutations=9999, action="add")df %>% mp_extract_internal_attr(name=adonis)#可视化PCoA结果p1 <- df %>% mp_plot_ord( .ord = pcoa, .group = Group, .color = Group, .size = 1.2, .alpha = 1, ellipse = T, show.legend = T # 是否显示置信圈的图例 ) + scale_fill_manual(values=c("#00A087FF", "#3C5488FF","red")) +  scale_color_manual(values=c("#00A087FF", "#3C5488FF","red")) # point的大小也可以映射到其他变量，如Observe或Shannon# 然后alpha多样性和beta多样性将同时显示df %<>% mp_cal_alpha(.abundance=RareAbundance)df$Shannonp2 <- df %>% mp_plot_ord( .ord = pcoa, .group = Group, .color = Group, .size = Observe, .alpha = Shannon, ellipse = TRUE, show.legend = FALSE ) + scale_fill_manual( values = c("#00A087FF", "#3C5488FF","red"), guide = guide_legend(keywidth=0.6, keyheight=0.6, label.theme=element_text(size=6.5)) ) + scale_color_manual( values=c("#00A087FF", "#3C5488FF","red"), guide = guide_legend(keywidth=0.6, keyheight=0.6, label.theme=element_text(size=6.5)) ) + scale_size_continuous( range=c(0.5, 3), guide = guide_legend(keywidth=0.6, keyheight=0.6, label.theme=element_text(size=6.5))  )#更换其他分组p3 <- df %>% mp_plot_ord( .ord = pcoa, .group = Saline, .color = Saline, .size = 1.2, .alpha = 1, ellipse = T, show.legend = T # 是否显示置信圈的图例 ) + scale_fill_manual(values=c("#00A087FF", "#3C5488FF")) +  scale_color_manual(values=c("#00A087FF", "#3C5488FF")) p4 <- df %>% mp_plot_ord( .ord = pcoa, .group = Type, .color = Type, .size = 1.2, .alpha = 1, ellipse = F, show.legend = T # 是否显示置信圈的图例 ) +  scale_fill_manual(values=c("#ff0000", "#fbb034","#ffdd00","#c1d82f","#00a4e4","#8a7967")) #拼图p1 + p2 p3 + p4

PS: 以上内容是小编个人学习代码笔记分享，仅供参考学习，欢迎大家一起交流学习。

参考：

1）http://bioconductor.org/packages/release/bioc/vignettes/MicrobiotaProcess/inst/doc/MicrobiotaProcess.html

2）https://chiliubio.github.io/microeco_tutorial/

3）https://doi.org/10.1371/journal.pone.0061217

更多推荐