其他
R可视化——基于MicrobiotaProcess包进行物种差异分析!
点击上方
“科研后花园”
关注我们
效果图:
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 Manipulation
library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics
library(phyloseq) # Handling and analysis of high-throughput microbiome census data
library(ggtree) # an R package for visualization of tree and annotation data
library(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)
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")
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 0
10 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
#物种相对丰度计算
df %<>%
mp_cal_abundance( # for each samples
.abundance = RareAbundance
) %>%
mp_cal_abundance( # for each groups
.abundance=RareAbundance,
.group=group
)
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"
)
# 结果提取
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))
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/