❝本节来介绍一款R包「ggpicrust2」,主要用于对「PICRUSt2」输出的结果进行深度操作,ggpicrust2集成了ko丰度到kegg通路丰度转换、通路注释、差异丰度(DA)分析及结果图的可视化等。下面小编来简单介绍下,更多详细的案例内容请参考作者官方文档。「此款R包安装依赖较多,请耐心安装」 ❞
❝https://github.com/cafferychen777/ggpicrust2 ❞
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
pkgs <- c("phyloseq", "ALDEx2", "SummarizedExperiment", "Biobase", "devtools",
"ComplexHeatmap", "BiocGenerics", "BiocManager", "metagenomeSeq",
"Maaslin2", "edgeR", "lefser", "limma", "KEGGREST", "DESeq2")
for (pkg in pkgs) {
if (!requireNamespace(pkg, quietly = TRUE))
BiocManager::install(pkg)
}
devtools::install_github("cafferychen777/ggpicrust2")
library(ggpicrust2)
library(tidyverse)
library(GGally)
library(ggprism)
library(patchwork)
library(ggh4x)
data("metacyc_abundance")
data("metadata")
metacyc_daa_results_df <- pathway_daa(abundance = metacyc_abundance %>%
column_to_rownames("pathway"), metadata = metadata,
group = "Environment", daa_method = "LinDA",
select = NULL, p.adjust = "BH", reference = NULL)
methods <- c("ALDEx2", "DESeq2", "edgeR")
daa_results_list <- lapply(methods, function(method) {
pathway_daa(abundance = metacyc_abundance %>%
column_to_rownames("pathway"), metadata = metadata, group = "Environment", daa_method = method)
})
method_names <- c("ALDEx2_Welch's t test","ALDEx2_Wilcoxon rank test","DESeq2", "edgeR")
# Compare results across different methods
comparison_results <- compare_daa_results(daa_results_list = daa_results_list, method_names = method_names)
> comparison_results
method num_features num_common_features num_diff_features
1 ALDEx2_Welch's t test 124 10 124
2 ALDEx2_Wilcoxon rank test 41 10 41
3 DESeq2 41 10 41
data("kegg_abundance")
data("metadata")
daa_results_df <- pathway_daa(abundance = kegg_abundance,
metadata = metadata, group = "Environment", daa_method = "LinDA")
daa_annotated_results_df <- pathway_annotation(pathway = "KO",
daa_results_df = daa_results_df, ko_to_kegg = TRUE)
data("ko_abundance")
data("metadata")
daa_results_df <- pathway_daa(abundance = ko_abundance %>%
column_to_rownames("#NAME"), metadata = metadata,
group = "Environment", daa_method = "LinDA")
daa_annotated_results_df <- pathway_annotation(pathway = "KO",
daa_results_df = daa_results_df, ko_to_kegg = FALSE)
data("ko_abundance")
data("metadata")
kegg_abundance <- ko2kegg_abundance(data = ko_abundance)
daa_results_df <- pathway_daa(kegg_abundance, metadata = metadata,
group = "Environment", daa_method = "LinDA")
daa_annotated_results_df <- pathway_annotation(pathway = "KO",
daa_results_df = daa_results_df, ko_to_kegg = TRUE)
p <- pathway_errorbar(abundance = kegg_abundance,
daa_results_df = daa_annotated_results_df,
Group = metadata$Environment,
ko_to_kegg = TRUE,
p_values_threshold = 0.05,
order = "pathway_class",
select = NULL,
p_value_bar = TRUE,
colors = NULL,
x_lab = "pathway_name")
data("metacyc_abundance")
data("metadata")
metacyc_daa_results_df <- pathway_daa(abundance = metacyc_abundance %>%
column_to_rownames("pathway"), metadata = metadata,
group = "Environment", daa_method = "LinDA")
metacyc_daa_annotated_results_df <- pathway_annotation(pathway = "MetaCyc",
daa_results_df = metacyc_daa_results_df, ko_to_kegg = FALSE)
p <- pathway_errorbar(abundance = metacyc_abundance %>%
column_to_rownames("pathway"),
daa_results_df = metacyc_daa_annotated_results_df,
Group = metadata$Environment,
ko_to_kegg = FALSE,
p_values_threshold = 0.05,
order = "group",
select = NULL,
p_value_bar = TRUE,
colors = NULL,
x_lab = "description")
data("metacyc_abundance")
data("metadata")
# Perform differential abundance analysis
metacyc_daa_results_df <- pathway_daa(
abundance = metacyc_abundance %>% column_to_rownames("pathway"),
metadata = metadata,
group = "Environment",
daa_method = "LinDA"
)
# Annotate the results
annotated_metacyc_daa_results_df <- pathway_annotation(
pathway = "MetaCyc",
daa_results_df = metacyc_daa_results_df,
ko_to_kegg = FALSE
)
feature_with_p_0.05 <- metacyc_daa_results_df %>%
filter(p_adjust < 0.05)
pathway_heatmap(
abundance = metacyc_abundance %>%
right_join(
annotated_metacyc_daa_results_df %>% select(all_of(c("feature","description"))),
by = c("pathway" = "feature")
) %>%
filter(pathway %in% feature_with_p_0.05$feature) %>%
select(-"pathway") %>%
column_to_rownames("description"),
metadata = metadata,
group = "Environment"
)
pathway_pca(abundance = metacyc_abundance %>% column_to_rownames("pathway"),
metadata = metadata, group = "Environment")