Skip to contents

Gene Set Enrichment Analysis for PICRUSt2 output

Source: R/pathway_gsea.R
pathway_gsea.Rd

This function performs Gene Set Enrichment Analysis (GSEA) on PICRUSt2 predicted functional data to identify enriched pathways between different conditions.

Usage

pathway_gsea(
  abundance,
  metadata,
  group,
  pathway_type = "KEGG",
  method = "camera",
  covariates = NULL,
  contrast = NULL,
  inter.gene.cor = 0.01,
  rank_method = "signal2noise",
  nperm = 1000,
  min_size = 5,
  max_size = 500,
  p_adjust_method = "BH",
  seed = 42,
  go_category = "all",
  organism = "ko",
  p.adjust = NULL,
  comparison = NULL,
  transformation = c("voom", "logCPM"),
  gene_sets = NULL
)

Arguments

abundance

A data frame containing gene/enzyme abundance data, with feature IDs in row names and samples as columns. Data frames may also provide a leading non-numeric feature ID column (for example #NAME, feature, or pathway); it is converted to row names automatically. Values must be finite, non-missing, and non-negative count-like abundances; negative or non-finite values are rejected rather than being coerced to zero. For KEGG analysis: features should be KO IDs (e.g., K00001). For MetaCyc analysis: features should be EC numbers (e.g., EC:1.1.1.1 or 1.1.1.1), NOT pathway IDs. MetaCyc pathway-like identifiers are rejected because GSEA requires gene/enzyme-level input; use pathway_daa for pathway-level MetaCyc abundance tables. For GO analysis: features should be KO IDs that will be mapped to GO terms. NOTE: This function requires gene-level data, not pathway-level abundances. For pathway abundance analysis, use pathway_daa instead

metadata

A data frame containing sample metadata. After sample alignment, all retained samples must have non-missing, non-empty values in the grouping column.

group

A character string specifying the column name in metadata that contains the grouping variable

pathway_type

A single character string specifying the pathway type: "KEGG", "MetaCyc", or "GO"

method

A single character string specifying the GSEA method:

  • "camera": Competitive gene set test using limma's camera function (default). Uses an inter-gene correlation adjustment under a competitive null.

  • "fry": Fast approximation to rotation gene set testing using limma's fry function. Self-contained test that is computationally efficient.

  • "fgsea": Fast preranked GSEA implementation. Note: preranked methods may produce unreliable p-values due to not accounting for inter-gene correlations (Wu et al., 2012).

  • "GSEA" or "clusterProfiler": clusterProfiler's GSEA implementation.

covariates

A character vector specifying column names in metadata to use as covariates for adjustment. Only supported when method is "camera" or "fry"; supplying covariates with preranked methods is an error because those methods use a precomputed rank vector rather than a design matrix. Default is NULL (no covariates). Covariate values must be complete for all aligned samples; rows with missing model variables are rejected rather than being silently dropped by model.matrix(). The resulting design matrix must also be finite and full rank; constant or fully confounded covariates are rejected. Example: covariates = c("age", "sex", "BMI")

contrast

For "camera" or "fry" methods, specify the coefficient or contrast to test. Supplying contrast with preranked methods is an error; use comparison instead. Default NULL automatically tests the single non-reference group coefficient in two-group designs. Multi-group designs must specify this explicitly. A character value must exactly match a design column name or a non-reference group level; substring matching is not used. A numeric scalar is treated as a design-column index, and a numeric vector must have length equal to the number of design columns. Named numeric vectors are matched and reordered by design column names; unnamed numeric vectors are interpreted in design column order.

inter.gene.cor

Numeric value specifying the inter-gene correlation for camera method. Default is 0.01. Use NA to estimate correlation from data for each gene set.

rank_method

A single character string specifying the ranking statistic for preranked methods (fgsea, GSEA, clusterProfiler): "signal2noise", "t_test", "log2_ratio", or "diff_abundance"

nperm

An integer specifying the number of permutations (for clusterProfiler method only). The fgsea method uses adaptive multilevel splitting and does not require a fixed permutation count.

min_size

An integer specifying the minimum gene set size

max_size

An integer specifying the maximum gene set size

p_adjust_method

A character string specifying the p-value adjustment method

seed

An integer specifying the random seed for reproducibility

go_category

A single character string specifying GO category to use. "all" (default) uses all categories present in the reference data. Valid categories are determined by the reference data (currently MF and CC). See table(ko_to_go_reference$category) for available categories.

organism

Deprecated and has no effect. The KEGG and GO reference data bundled with ggpicrust2 are KO-based (organism-independent), so gene sets are returned in KO space regardless of this argument. Retained only for signature compatibility; passing any value other than the default "ko" emits a deprecation warning. Will be removed in a future release.

p.adjust

Deprecated alias for p_adjust_method. Do not supply both parameters with different values.

comparison

For preranked methods only ("fgsea", "GSEA", or "clusterProfiler"), an optional length-2 character vector c(group1, group2) defining the ranking direction. Ranking statistics are calculated as group1 versus group2: positive signal2noise, t_test, diff_abundance, and log2_ratio values indicate higher abundance in group1 (for log2_ratio, log2(group1 / group2)). If NULL, exactly two aligned group levels must be present and their factor-level order is used. Multi-group preranked analyses must specify comparison explicitly.

transformation

Transformation for camera/fry: "voom" (the historical default) uses count-dependent observation weights; "logCPM" uses log2(1e6 * abundance / column_total + 0.5) and abundance-trend variance moderation without voom weights. The latter is invariant to positive per-sample rescaling and can be used for predicted or relative abundances whose units are not observed read counts. It is a different statistical model, not a fallback after a voom failure or a guarantee of calibration. Requires positive sample totals. This parameter does not apply to preranked methods. Neither transformation accounts for uncertainty from the upstream functional prediction. Results also depend on the measured feature universe and, for camera, the correlation model.

gene_sets

Optional named list of feature identifiers for explicitly defined gene sets. When supplied, replaces the bundled sets for all methods; size filtering and multiple-testing adjustment still apply. Identifiers must match the normalized abundance row names. pathway_type still controls input identifier handling. Do not combine custom sets with organism or go_category; these are bundled-reference selection arguments.

Value

A data frame containing GSEA results with columns:

  • pathway_id: Pathway identifier

  • pathway_name: Pathway name/description

  • size: Number of genes in the pathway

  • direction: Direction of enrichment ("Up" or "Down", for camera/fry)

  • pvalue: Raw p-value

  • p.adjust: Adjusted p-value (FDR)

  • method: The method used for analysis

For fgsea/clusterProfiler methods, additional columns include ES, NES, leading_edge, group1, and group2. Positive ES/NES values are in the group1 versus group2 direction. For camera/fry methods, limma does not return NES; ggpicrust2 retains a legacy NES column containing a signed -log10(pvalue) score for visualization compatibility and labels it with score_type and score_label.

Details

Method Selection:

The camera method (default):

  • Uses an inter-gene correlation adjustment (default 0.01)

  • It supports covariate adjustment through the design matrix

  • It performs a competitive test (genes in set vs. genes not in set)

The fry method is a fast alternative that:

  • Performs a self-contained test (are genes in the set differentially expressed?)

  • Is computationally very efficient for large numbers of gene sets

  • Also supports covariate adjustment

Preranked methods use a ranked-feature null rather than camera's competitive model with a correlation adjustment. The hypotheses and input assumptions differ; neither method is universally calibrated for predicted abundances. For these methods, positive ES/NES values indicate gene-set enrichment near the top of the ranked list. With comparison = c(group1, group2), the top of the list corresponds to features higher in group1; reverse comparison to reverse the biological direction.

Covariate Adjustment:

When using method = "camera" or method = "fry", you can adjust for confounding variables by specifying them in the covariates parameter. This is particularly important in microbiome studies where factors like age, sex, BMI, and batch effects can influence results.

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., Lim, E., Vaillant, F., Asselin-Labat, M. L., Visvader, J. E., & Smyth, G. K. (2010). ROAST: rotation gene set tests for complex microarray experiments. Bioinformatics, 26(17), 2176-2182.

Examples

if (FALSE) { # \dontrun{
# Load example data
data(ko_abundance)
data(metadata)

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

# Method 1: Using camera (recommended) - accounts for inter-gene correlations
gsea_results <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "camera"
)

# Method 2: Using camera with covariate adjustment
gsea_results_adj <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  covariates = "Mouse_Sex",
  pathway_type = "KEGG",
  method = "camera"
)

# Method 3: Using fry for fast self-contained testing
gsea_results_fry <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "fry"
)

# Method 4: A different null, using a prespecified preranked comparison
gsea_results_fgsea <- pathway_gsea(
  abundance = abundance_data,
  metadata = metadata,
  group = "Environment",
  pathway_type = "KEGG",
  method = "fgsea",
  comparison = c("Pro-survival", "Pro-inflammatory"),
  seed = 42
)

# Visualize results
visualize_gsea(gsea_results, plot_type = "enrichment_plot", n_pathways = 10)
} # }