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"))Method Selection Guide
Choose the hypothesis and input scale before comparing p-values.
| Method | Null hypothesis / input | Covariates | Score |
|---|---|---|---|
camera (default) |
Competitive gene-set test on a voom fit | Yes | Signed -log10(raw p-value) |
fry |
Self-contained gene-set test on a voom fit | Yes | Signed -log10(raw p-value) |
fgsea |
Preranked gene-set enrichment | No | Normalized enrichment score (NES) |
clusterProfiler |
Preranked GSEA implementation | No | NES |
By default, camera/fry run limma::voom() on
non-negative, count-like input. Do not pass log-transformed or relative
abundances as if they were counts. The default camera correlation
setting is inter.gene.cor = 0.01; it is a supplied
correlation adjustment, not an estimate for every set. When features in
a set are strongly correlated, fixing their correlation at 0.01 can
underestimate the variance and produce anti-conservative p-values. The
camera examples below therefore explicitly estimate within-set
correlation with inter.gene.cor = NA_real_. This addresses
that modeling assumption; it does not establish calibration for every
predicted functional profile. For predicted or relative abundances, the
explicit option transformation = "logCPM" instead uses
log2(1e6 * abundance / column_total + 0.5) with
abundance-trend variance moderation and no voom observation weights. It
is invariant to positive per-sample rescaling, but retains compositional
effects and does not account for upstream prediction uncertainty. The
default remains "voom". 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",
inter.gene.cor = NA_real_,
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 contains shared KOs assigned to disease and
eukaryotic pathway maps as well as microbial pathways. GSEA retains
these maps, whereas ko2kegg_abundance() applies its
prokaryote category filter by default. A disease name among the top
results is therefore a label for a tested KO set, not evidence that the
microbiome carries out a host disease pathway. Inspect its member KOs
and coverage before giving it a biological interpretation.
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 = "camera",
inter.gene.cor = NA_real_
)
# The results now reflect the group effect after adjusting for confounders
head(gsea_results_adjusted)Fast Analysis with fry
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_labelsBarplot
# 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
barplotDotplot
# 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
dotplotEnrichment-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_plotRidge Plot
A ridge plot shows the distribution of member-KO group-mean abundance ratios. It does not reproduce covariate-adjusted model coefficients, voom-weighted 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_plotThe 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$resultsInterpretation
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.