R可视化——基于MicrobiotaProcess包进行物种差异分析！

Original 王志山科研后花园 2023-09-08

收录于合集 #R可视化 69个

点击上方

“科研后花园”

关注我们

效果图：

1、设置工作环境并加载R包：

rm(list = ls())setwd("D:/桌面/物种差异分析")#加载R包

library(MicrobiotaProcess) # A comprehensive R package for managing and analyzing microbiome and other ecological data within the tidy framework

library(dplyr) # A Grammar of Data Manipulationlibrary(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphicslibrary(phyloseq) # Handling and analysis of high-throughput microbiome census datalibrary(ggtree) # an R package for visualization of tree and annotation datalibrary(ggtreeExtra) # An R Package To Add Geometric Layers On Circular Or Other Layout Tree Of "ggtree"library(ggstar) # Multiple Geometric Shape Point Layer for 'ggplot2'library(forcats) # Tools for Working with Categorical Variables (Factors)

2、加载数据（与此前进行PCoA分析一致，加载OTU特征表、样品信息表及OTU物种注释表）：

sample <- read.table("sample.txt",check.names = F, row.names = 1, header = 1, sep = "\t")OTU<- read.table("otu.txt",check.names = F, row.names = 1, header = 1, sep = "\t")Tax <- read.table("tax.txt",check.names = F, row.names = 1, header = 1, sep = "\t")

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

ps <- phyloseq(sample_data(sample), otu_table(as.matrix(OTU), taxa_are_rows=TRUE), tax_table(as.matrix(Tax)))

#转换数据格式df <- ps %>% as.MPSE()df# A MPSE-tibble (MPSE object) abstraction: 720# × 20# OTU=60 | Samples=12 | Assays=Abundance,# RareAbundance, RelRareAbundanceBySample,# RelRareAbundanceBySample |# Taxonomy=Kingdom, Phylum, Class, Order,# Family, Genus, Species OTU Sample Abundance RareAbundance <chr> <chr> <int> <int> 1 OTU1 A_1 421 410 2 OTU2 A_1 7 5 3 OTU3 A_1 116 114 4 OTU4 A_1 112 109 5 OTU5 A_1 22 22 6 OTU6 A_1 143 136 7 OTU7 A_1 11 11 8 OTU8 A_1 6 6 9 OTU9 A_1 0 010 OTU10 A_1 1 1# ℹ 710 more rows# ℹ 16 more variables:# RelRareAbundanceBySample <dbl>,# RelRareAbundanceBySample.y <dbl>,# group <chr>, Kingdom <chr>, Phylum <chr>,# Class <chr>, Order <chr>, Family <chr>,# Genus <chr>, Species <chr>, …# ℹ Use `print(n = ...)` to see more rows

4、计算物种的相对丰度：

#物种相对丰度计算df %<>% mp_cal_abundance( # for each samples .abundance = RareAbundance ) %>% mp_cal_abundance( # for each groups .abundance=RareAbundance, .group=group )

5、物种差异分析：

df %<>% mp_diff_analysis( .abundance = RelRareAbundanceBySample, .group = group, tip.level = "OTU", force = FALSE, relative = TRUE, taxa.class = "all", first.test.method = "kruskal.test", first.test.alpha = 0.05, p.adjust = "fdr", filter.p = "fdr", strict = TRUE, fc.method = "generalizedFC", second.test.method = "wilcox.test", second.test.alpha = 0.05, cl.min = 4, cl.test = TRUE, subcl.min = 3, subcl.test = TRUE, ml.method = "lda",# 'lda' or 'rf' normalization = 1e+06, ldascore = 2,#LDA阈值 bootnums = 30, sample.prop.boot = 0.7, ci = 0.95, seed = 123, type = "species" )

6、提取结果并基于ggtree等R包进行可视化：

# 结果提取taxa.tree <- df %>% mp_extract_tree(type="taxatree")taxa.tree

taxa.tree %>% select(label, nodeClass, LDAupper, LDAmean, LDAlower, Sign_group, pvalue, fdr) %>% dplyr::filter(!is.na(fdr))

##通过ggtree和ggtreeExtra可视化ggtree( taxa.tree, layout="radial", size = 0.3) + geom_point(data = td_filter(!isTip), fill="white", size=1, shape=21)+ geom_hilight( data = td_filter(nodeClass == "Phylum"), mapping = aes(node = node, fill = label))+ ggnewscale::new_scale("fill") + geom_fruit( data = td_unnest(RareAbundanceBySample), geom = geom_star, mapping = aes( x = fct_reorder(Sample, group, .fun=min), size = RelRareAbundanceBySample, fill = group, subset = RelRareAbundanceBySample > 0), starshape = 13, starstroke = 0.25, offset = 0.04, pwidth = 0.8, grid.params = list(linetype=2)) + scale_size_continuous( name="Relative Abundance (%)", range = c(.5, 3) ) + scale_fill_manual(values=c("#1B9E77", "#D95F02","blue"))+ ggnewscale::new_scale("fill") + geom_fruit( geom = geom_col, mapping = aes( x = LDAmean, fill = Sign_group, subset = !is.na(LDAmean)), orientation = "y", offset = 0.05, pwidth = 0.5, axis.params = list(axis = "x", title = "Log10(LDA)", title.height = 0.01, title.size = 2, text.size = 1.8, vjust = 1), grid.params = list(linetype = 2))+ geom_tiplab(size=2, offset=11.2)+ ggnewscale::new_scale("size") + geom_point( data=td_filter(!is.na(Sign_group)), mapping = aes(size = -log10(fdr), fill = Sign_group, ), shape = 21, ) + scale_size_continuous(range=c(1, 3)) + scale_fill_manual(values=c("#1B9E77", "#D95F02","blue"))+ theme( legend.key.height = unit(0.3, "cm"), legend.key.width = unit(0.3, "cm"), legend.spacing.y = unit(0.02, "cm"), legend.text = element_text(size = 7), legend.title = element_text(size = 9))

7、基于MicrobiotaProcess包中已有函数进行可视化（简化代码）：

df %>% mp_plot_diff_res( group.abun = TRUE, pwidth.abun=0.1 ) + scale_fill_manual(values=c("deepskyblue", "orange","red")) + scale_fill_manual( aesthetics = "fill_new", values = c("deepskyblue", "orange","red") ) + scale_fill_manual( aesthetics = "fill_new_new", values = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "#999999" ) )

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

测试数据获取：

夸克网盘（手机端下载夸克APP进行搜索提取）：

链接：https://pan.quark.cn/s/db89e3407fd8

提取码：Kmab

参考：

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

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

更多推荐

R绘图模板——气泡图+注释+配对连线！！！

R绘图模板——世界地图的绘制及拓展！！！

基于MicrobiotaProcess包进行PCoA分析！