Using ggpicrust2

Introduction

ggpicrust2 provides a practical workflow for PICRUSt2 downstream analysis:

This vignette focuses on the general package workflow. For a deeper GSEA walkthrough, see the dedicated gsea_analysis vignette.

Installation and example data

Install the package and the optional backends used by this tutorial:

install.packages(c("ggpicrust2", "MicrobiomeStat", "BiocManager"))
BiocManager::install(c("KEGGREST", "limma"))
library(ggpicrust2)
library(tibble)
data("ko_abundance")
data("metadata")
alpha <- 0.05

Both workflows below use the same data, LinDA method, and adjusted p-value threshold. LinDA accepts the continuous abundance estimates produced by ko2kegg_abundance(), but its default count-data winsorization rounds values inside MicrobiomeStat. In a direct pathway_daa() call, set linda_winsor = FALSE to preserve fractional values and linda_adaptive = FALSE for fixed pseudo-count handling. A fixed pseudo-count still depends on the units of the input when zeros are present. Count-based methods such as ALDEx2 require integer input and the package rounds non-integer input with a warning. Choose the method for its assumptions and your study design, not to obtain significance. KEGG pathway annotation requires internet access and KEGGREST.

One-command workflow

results <- ggpicrust2(
  data = ko_abundance,
  metadata = metadata,
  group = "Environment",
  pathway = "KO",
  daa_method = "LinDA",
  ko_to_kegg = TRUE,
  order = "pathway_class",
  p_values_bar = TRUE,
  p_values_threshold = alpha,
  x_lab = "pathway_name"
)

# A method's plot is NULL when no pathways can be plotted.
results[[1]]$plot
head(results[[1]]$results)

Stepwise pathway workflow

Run the installation and example-data setup above first. This workflow uses the same analysis settings as the one-command workflow.

Convert KO abundance to KEGG pathway abundance

kegg_pathway_abundance <- ko2kegg_abundance(data = ko_abundance)
head(kegg_pathway_abundance[, 1:3])

Match group labels to sample identifiers

pathway_daa(), pathway_heatmap(), and pathway_pca() accept metadata and a column name, such as group = "Environment". In contrast, pathway_errorbar() has no metadata argument: its capitalized Group parameter requires one group label per abundance column, not a column name.

The bundled metadata and abundance table have different sample orders. Name the group vector with sample IDs so that plotting aligns labels to the correct samples. For your own data, replace sample_name and Environment with your sample-ID and grouping columns. Do not pass an unnamed metadata column unless you have already verified its order against the abundance columns.

stopifnot(
  !anyNA(metadata$sample_name),
  !anyDuplicated(metadata$sample_name),
  setequal(colnames(kegg_pathway_abundance), metadata$sample_name)
)
sample_groups <- setNames(metadata$Environment, metadata$sample_name)

Run differential abundance analysis

daa_results <- pathway_daa(
  abundance = kegg_pathway_abundance,
  metadata = metadata,
  group = "Environment",
  daa_method = "LinDA"
)

head(daa_results)

Annotate pathway results

annotated_daa <- pathway_annotation(
  pathway = "KO",
  daa_results_df = daa_results,
  ko_to_kegg = TRUE,
  p_adjust_threshold = alpha
)

head(annotated_daa)

Visualize pathway-level results

ko_to_kegg = TRUE is required here too: the rows are KEGG pathways and use pathway_name annotations. This flag does not reconvert the abundance matrix in the plotting function.

sig_pathways <- unique(annotated_daa$feature[
  !is.na(annotated_daa$p_adjust) & annotated_daa$p_adjust < alpha
])

p <- NULL
if (length(sig_pathways) > 0) {
  p <- pathway_errorbar(
    abundance = kegg_pathway_abundance,
    daa_results_df = annotated_daa,
    Group = sample_groups,
    ko_to_kegg = TRUE,
    p_values_threshold = alpha,
    order = "pathway_class",
    x_lab = "pathway_name"
  )
} else {
  message("No pathways pass the adjusted p-value threshold; skipping the error bar plot.")
}
p

No significant pathways is a valid analysis outcome. Keep the results table; do not increase the threshold or change methods just to produce a plot. Missing KEGG annotations can also prevent plotting even when significant results exist; check the annotation warnings separately.

if (length(sig_pathways) > 0) {
  pathway_heatmap(
    abundance = kegg_pathway_abundance[sig_pathways, , drop = FALSE],
    metadata = metadata,
    group = "Environment"
  )
}

pathway_pca(
  abundance = kegg_pathway_abundance,
  metadata = metadata,
  group = "Environment"
)

Using ALDEx2 instead

ALDEx2 is an optional Bioconductor dependency. For two groups it returns both Welch and Wilcoxon results, so select one test before annotation and plotting. Otherwise a pathway can occur twice with different p-values. Specify the test in advance. ALDEx2 uses Monte Carlo sampling, so set a seed for reproducibility; results can still differ across package versions.

# Install once with BiocManager::install("ALDEx2").
set.seed(207)
aldex_results <- pathway_daa(
  abundance = kegg_pathway_abundance,
  metadata = metadata,
  group = "Environment",
  daa_method = "ALDEx2"
)
daa_results <- aldex_results[
  aldex_results$method == "ALDEx2_Welch's t test", , drop = FALSE
]

Then rerun the annotation and visualization steps with this daa_results. For more than two groups or multiple contrasts, inspect method, group1, and group2 and select the supported test/contrast explicitly. An ALDEx2 run can have no adjusted p-values below 0.05 even when LinDA finds some; these are different statistical procedures, not equivalent plotting modes.

Taxa contribution workflow

PICRUSt2 contribution files attribute predicted functional abundance to taxa. ggpicrust2 supports both gene-family-level and pathway-level contribution workflows.

Run a synthetic contribution example

This small example illustrates input schemas and aggregation. It is not a biological result, and its sample IDs are separate from the bundled KO dataset.

contrib_input <- expand.grid(
  sample = paste0("S", 1:4),
  function_id = c("K00001", "K00002"),
  taxon = c("ASV1", "ASV2"),
  stringsAsFactors = FALSE
)
contrib_input$taxon_function_abun <- seq_len(nrow(contrib_input))
contrib_data <- read_contrib_file(data = contrib_input)
contrib_metadata <- data.frame(
  sample_name = paste0("S", 1:4),
  Environment = rep(c("Control", "Treatment"), each = 2)
)
taxonomy <- data.frame(
  ASV = c("ASV1", "ASV2"),
  Genus = c("ExampleGenusA", "ExampleGenusB")
)
taxa_contrib <- aggregate_taxa_contributions(
  contrib_data, taxonomy = taxonomy, tax_level = "Genus", top_n = 2
)
head(taxa_contrib)

The aggregation sums the selected contribution column. By default it prefers norm_taxon_function_contrib when supplied; this example supplies raw taxon_function_abun only. A percentage bar subsequently normalizes within each sample/function, so its heights describe the taxonomic composition of that function, not a between-function abundance comparison.

taxa_contribution_bar(
  contrib_agg = taxa_contrib,
  metadata = contrib_metadata,
  group = "Environment",
  facet_by = "function"
)
taxa_contribution_heatmap(contrib_agg = taxa_contrib, n_functions = 2)

Read your own PICRUSt2 files

For real data, replace the synthetic input with one of the readers below and use metadata and taxonomy for those same samples and taxa:

# KO/gene-family contributions:
# contrib_data <- read_contrib_file("pred_metagenome_contrib.tsv")
# Pathway contributions (often MetaCyc):
# contrib_data <- read_pathway_contrib_file("path_abun_contrib.tsv.gz")
# Wide stratified abundance:
# contrib_data <- read_strat_file("pred_metagenome_strat.tsv")

When optional daa_results_df or pathway_ids filters contain KEGG pathway IDs and the contribution table is KO-level, the function expands those pathway IDs to KO members and retains matching KO rows. It does not turn KO contributions into pathway contributions. The output function_id remains a KO identifier. For pathway-level MetaCyc contributions, matching MetaCyc IDs are filtered directly. Do not interpret a member-KO filter as independent evidence that a taxon drives a reconstructed pathway’s activity.

For pathway-level data, use matching pathway annotations. For example:

path_input <- contrib_input
path_input$function_id <- ifelse(path_input$function_id == "K00001",
                                  "GLYCOLYSIS", "PWY-5484")
path_data <- read_pathway_contrib_file(data = path_input)
path_taxa_contrib <- aggregate_taxa_contributions(
  path_data, taxonomy = taxonomy, tax_level = "Genus", top_n = 2
)
pathway_annotation_df <- pathway_annotation(
  data = data.frame(function_id = unique(path_taxa_contrib$function_id)),
  pathway = "MetaCyc"
)

GSEA workflow

Use GSEA when you want pathway-set level inference from KO or EC abundance rather than testing each pathway independently.

gsea_results <- pathway_gsea(
  abundance = ko_abundance %>% column_to_rownames("#NAME"),
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "camera"
)

annotated_gsea <- gsea_pathway_annotation(
  gsea_results = gsea_results,
  pathway_type = "KEGG"
)

visualize_gsea(
  gsea_results = annotated_gsea,
  plot_type = "barplot",
  n_pathways = 15
)

For a method-by-method GSEA explanation, covariate adjustment, and comparison with DAA, see the gsea_analysis vignette.

Summary

The package is easiest to use when you choose the shortest path that matches your question: