Gene Set Enrichment Analysis with ggpicrust2

Introduction

This vignette tests sets of KO identifiers in PICRUSt2 predicted functional profiles. It covers method-specific hypotheses and scores, group contrasts, covariate adjustment, and visualizations. The input is predicted abundance, not measured gene expression or pathway activity.

Installation

The examples use bundled KO abundance and metadata. Install the optional backends once, then load the package:

install.packages(c("ggpicrust2", "MicrobiomeStat", "ggridges", "ggVennDiagram",
                   "circlize", "igraph", "BiocManager"))
BiocManager::install(c("limma", "fgsea", "ComplexHeatmap"))

library(ggpicrust2)

Method Selection Guide

Choose the hypothesis and input scale before comparing p-values.

Method	Null hypothesis / input	Covariates	Score
`fry` (default)	Self-contained gene-set test on log-CPM abundance	Yes	Signed -log10(raw p-value)
`camera`	Competitive gene-set test on log-CPM abundance	Yes	Signed -log10(raw p-value)
`fgsea`	Preranked gene-set enrichment	No	Normalized enrichment score (NES)
`clusterProfiler`	Preranked GSEA implementation	No	NES

Since version 2.6.0, fry is the default method. It asks whether the KOs in a set change in a consistent direction between groups. camera asks whether a set changes more than the other KOs. In a predicted metagenome, a shift in a dominant taxon moves many KOs together, so a competitive test rarely singles out a set even when its members change. Choose camera when relative enrichment is the question.

By default (since version 2.6.0), camera/fry analyze log2(1e6 * abundance / column_total + 0.5) with abundance-trend variance moderation, and camera estimates the inter-gene correlation separately for each set (inter.gene.cor = NA). Both choices follow from how PICRUSt2 builds its output. Each predicted KO abundance is a sum of taxon abundances multiplied by reference-genome copy numbers, so the table has no read-count units, and KOs carried by the same genomes move together. The log-CPM model gives the same result when a table is multiplied by a constant, and the estimated correlation prevents camera from treating correlated members of a set as independent evidence. Neither choice removes the uncertainty of the upstream prediction.

The defaults before 2.6.0 remain available explicitly: transformation = "voom" for genuine read counts, and a fixed inter.gene.cor such as 0.01 when a small correlation is justified for the data. A fixed small correlation underestimates the variance of strongly correlated sets and makes camera anti-conservative. An optional named list, gene_sets, replaces the bundled reference sets; its member IDs must match the input feature IDs, and the same size filters and multiple-testing adjustment apply. The limma documentation explains the underlying tests; the fgsea documentation explains preranked enrichment and leading-edge genes.

Preranked statistics are calculated from the supplied abundance, without an implicit voom or library-size normalization step. Account for the input scale in the study design. Tied ranks can affect the leading edge; a fixed seed makes a run reproducible but does not remove that ambiguity.

Competitive analysis with camera

This example uses the competitive camera test:

# Load example data
data(ko_abundance)
data(metadata)
metadata$Environment <- factor(
  metadata$Environment, levels = c("Pro-inflammatory", "Pro-survival")
)

# Prepare abundance data
abundance_data <- as.data.frame(ko_abundance)
rownames(abundance_data) <- abundance_data[, "#NAME"]
abundance_data <- abundance_data[, -1]

# Run the competitive camera test
gsea_results <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "camera",
  min_size = 5,
  max_size = 500,
  p_adjust_method = "BH"
)

# View the top results
head(gsea_results)

For this factor order, camera/fry test Pro-survival minus Pro-inflammatory. direction = "Up" refers to that contrast. The compatibility column NES is signed -log10(pvalue) for these methods; it is neither a true NES nor an effect size. score_type and score_label identify its meaning.

table(gsea_results$direction)
sum(gsea_results$p.adjust < 0.05, na.rm = TRUE)
unique(gsea_results[, c("method", "score_type", "score_label")])

The KEGG reference also assigns bacterial KOs to host disease and other eukaryotic pathway maps. Since version 2.6.0, KEGG gene sets apply the same prokaryote filter as ko2kegg_abundance() by default, so GSEA and pathway DAA test the same maps: those that KEGG draws for prokaryotic genomes, without Organismal Systems and Human Diseases maps other than bacterial infection and antimicrobial resistance. With filter_for_prokaryotes = FALSE, a disease name among the top results is a label for a tested KO set, not evidence that the microbiome carries out a host disease pathway.

Covariate Adjustment

One of the most powerful features of the camera and fry methods is the ability to adjust for confounding variables. This is particularly important in microbiome studies where factors like age, sex, BMI, and batch effects can influence results.

# Mouse_Sex is observed in the bundled metadata and varies within both groups.

gsea_results_adjusted <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  covariates = "Mouse_Sex",
  pathway_type = "KEGG",
  method = "fry"
)

# The results now reflect the group effect after adjusting for confounders
head(gsea_results_adjusted)

Sequencing depth can act as a confounder too. When one group was sequenced more deeply than the other and some samples have few reads, zeros and low counts fall unevenly between the groups, and fry on logCPM can report pathways that do not change. Supply the read depth of the amplicon data through sequencing_depth; its centered log is added to the design. Use the column sums of the ASV table given to PICRUSt2, not the column sums of the predicted table, which depend on community composition.

# asv_counts: the ASV-by-sample count table used as PICRUSt2 input
depth <- colSums(asv_counts)

gsea_results_depth <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  covariates = "Mouse_Sex",
  sequencing_depth = depth
)

Self-contained analysis with fry (default)

Use fry for the self-contained null that genes in the set have no group effect. It does not test camera’s competitive null, so discovery counts can differ without a software error.

# Fast rotation gene set test
gsea_results_fry <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "fry",
  min_size = 5,
  max_size = 500
)

head(gsea_results_fry)

Preranked GSEA (fgsea)

For a prespecified ranked-list analysis, fgsea returns true NES values and leading-edge genes. The explicit comparison below makes a positive score mean higher abundance in Pro-survival, matching the camera/fry contrast above. The input and ranking choices are part of this demonstration, not a universal recommendation for every study.

# Preranked testing uses a different null from camera/fry.
gsea_results_fgsea <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "fgsea",
  rank_method = "signal2noise",
  comparison = c("Pro-survival", "Pro-inflammatory"),
  min_size = 10,
  max_size = 500,
  p_adjust_method = "BH",
  seed = 42
)

# View the top results
head(gsea_results_fgsea)

Annotating GSEA Results

To make the results more interpretable, we can annotate them with pathway names and descriptions:

# Annotate GSEA results
annotated_results <- gsea_pathway_annotation(
  gsea_results = gsea_results,
  pathway_type = "KEGG"
)

# View the annotated results
head(annotated_results)

Visualizing GSEA Results

The ggpicrust2 package provides several visualization options for GSEA results. The visualize_gsea() function uses annotation names when available. It displays the top n_pathways after sorting; it does not automatically filter by significance. Inspect adjusted p-values before calling displayed pathways significant.

Pathway Label Options

The visualize_gsea() function offers flexible pathway labeling:

# Option 1: Use raw GSEA results (shows pathway IDs)
plot_with_ids <- visualize_gsea(
  gsea_results = gsea_results,
  plot_type = "barplot",
  n_pathways = 10
)

# Option 2: Use annotated results (automatically shows pathway names)
plot_with_names <- visualize_gsea(
  gsea_results = annotated_results,
  plot_type = "barplot",
  n_pathways = 10
)

# Option 3: Explicitly specify which column to use for labels
plot_custom_labels <- visualize_gsea(
  gsea_results = annotated_results,
  plot_type = "barplot",
  pathway_label_column = "pathway_name",
  n_pathways = 10
)

# Compare the plots
plot_with_ids
plot_with_names
plot_custom_labels

Barplot

# Create a barplot of the top-ranked pathways
barplot <- visualize_gsea(
  gsea_results = annotated_results,
  plot_type = "barplot",
  n_pathways = 20,
  sort_by = "p.adjust"
)

# Display the plot
barplot

Dotplot

# Create a dotplot of the top-ranked pathways
dotplot <- visualize_gsea(
  gsea_results = annotated_results,
  plot_type = "dotplot",
  n_pathways = 20,
  sort_by = "p.adjust"
)

# Display the plot
dotplot

Enrichment-score summary

# This is a score-summary bar chart, not a running enrichment curve.
enrichment_plot <- visualize_gsea(
  gsea_results = annotated_results,
  plot_type = "enrichment_plot",
  n_pathways = 10,
  sort_by = "p.adjust"
)

# Display the plot
enrichment_plot

Ridge Plot

A ridge plot shows the distribution of member-KO group-mean abundance ratios. It does not reproduce covariate-adjusted model coefficients, model-based effects, or a gene-set significance test. A pathway can contain KO ratios of both signs even when its test has one enrichment direction.

# Create a ridge plot for GSEA results
# Note: Requires ggridges package to be installed
ridge_plot <- pathway_ridgeplot(
  gsea_results = gsea_results,
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  comparison = c("Pro-inflammatory", "Pro-survival"),
  n_pathways = 10,
  sort_by = "p.adjust",
  show_direction = TRUE,
  colors = c("Down" = "#3182bd", "Up" = "#de2d26")
)

# Display the plot
ridge_plot

The ridge plot shows: - Each pathway as a density ridge - Color indicates enrichment direction (Up = red, Down = blue) - The log2 ratio of mean supplied KO abundance, Pro-survival / Pro-inflammatory - A data-derived pseudocount added to both group means; see ?pathway_ridgeplot - A vertical dashed line at 0 for reference

Leading-edge network and heatmap

Only preranked methods provide leading-edge genes. camera/fry return an empty leading_edge field; passing those results to these displays produces no leading-edge information. Do not interpret an empty graph as no biological relationships. Use the preranked result created above for this example:

annotated_fgsea <- gsea_pathway_annotation(gsea_results_fgsea, pathway_type = "KEGG")
leading_results <- annotated_fgsea[
  !is.na(annotated_fgsea$leading_edge) & nzchar(annotated_fgsea$leading_edge), , drop = FALSE
]
if (nrow(leading_results) > 0) {
  network_plot <- visualize_gsea(
    leading_results, plot_type = "network", n_pathways = 10,
    network_params = list(similarity_measure = "jaccard", similarity_cutoff = 0.2)
  )
  print(network_plot)
  leading_heatmap <- visualize_gsea(
    leading_results, plot_type = "heatmap", n_pathways = 10,
    abundance = abundance_data, metadata = metadata, group = "Environment",
    heatmap_params = list(cluster_rows = TRUE, cluster_columns = TRUE,
                          show_rownames = TRUE)
  )
  ComplexHeatmap::draw(leading_heatmap)
}

Edges measure overlap between leading-edge sets, not biochemical interaction or causation. Each heatmap row is the mean supplied abundance over one pathway’s matched leading-edge genes, standardized across samples. It is not a row per gene or measured gene expression. Overlapping sets reuse genes.

Comparing GSEA and DAA Results

Compare results at the same identifier level: GSEA returns KEGG pathway IDs, so aggregate KO abundance to KEGG pathways before DAA. Annotation changes labels, not the unit of analysis. The two procedures test different hypotheses; overlap is descriptive agreement, not independent validation. This example uses LinDA consistently with the general workflow tutorial. Restrict the comparison to the shared tested universe: gene-set size filters and pathway filters can otherwise make an untested pathway look like a method-specific discovery. P-values retain their original analysis-wide adjustment.

# Compare KEGG pathways to KEGG pathways, not individual KO identifiers.
kegg_pathway_abundance <- ko2kegg_abundance(data = ko_abundance)
daa_results <- pathway_daa(
  abundance = kegg_pathway_abundance,
  metadata = metadata,
  group = "Environment",
  daa_method = "LinDA"
)

# Compare only pathways that both procedures actually tested.
# Keep each analysis's original multiple-testing adjustment.
common_pathways <- intersect(annotated_results$pathway_id, daa_results$feature)
gsea_common <- annotated_results[annotated_results$pathway_id %in% common_pathways, , drop = FALSE]
daa_common <- daa_results[daa_results$feature %in% common_pathways, , drop = FALSE]
comparison <- compare_gsea_daa(
  gsea_results = gsea_common,
  daa_results = daa_common,
  plot_type = "venn",
  p_threshold = 0.05
)

# Display the comparison plot
comparison$plot

# View the comparison results
comparison$results

Interpretation

Use a prespecified contrast and a method whose input assumptions fit the study. Camera, fry, and preranked tests ask different questions, so do not select the method with the most significant pathways. The same predicted abundance data underlie DAA and GSEA; overlap is descriptive agreement. Report the tested feature universe, pathway-size filters, method, score type, and adjusted p-value threshold. These outputs do not validate actual pathway activity or propagate PICRUSt2 prediction uncertainty.

References

Wu, D., & Smyth, G. K. (2012). Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Research, 40(17), e133.
Wu, D., et al. (2010). ROAST: rotation gene set tests for complex microarray experiments. Bioinformatics, 26(17), 2176-2182.