--- title: "Gene Set Enrichment Analysis with ggpicrust2" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Gene Set Enrichment Analysis with ggpicrust2} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) ``` ## 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: ```{r installation, eval=FALSE} install.packages(c("ggpicrust2", "MicrobiomeStat", "ggridges", "ggVennDiagram", "circlize", "igraph", "BiocManager")) BiocManager::install(c("limma", "fgsea", "ComplexHeatmap")) ``` ```{r setup, eval=FALSE} library(ggpicrust2) ``` ## 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](https://bioconductor.org/packages/limma) explains the underlying tests; the [fgsea documentation](https://bioconductor.org/packages/fgsea) 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: ```{r basic-gsea, eval=FALSE} # 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. ```{r inspect-results, eval=FALSE} 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. ```{r covariate-gsea, eval=FALSE} # 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. ```{r fry-gsea, eval=FALSE} # 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. ```{r fgsea, eval=FALSE} # 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: ```{r annotate-gsea, eval=FALSE} # 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: ```{r pathway-labels, eval=FALSE} # 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 ```{r barplot, eval=FALSE} # 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 ```{r dotplot, eval=FALSE} # 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 ```{r enrichment-plot, eval=FALSE} # 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, 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. ```{r ridge-plot, eval=FALSE} # 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: ```{r leading-edge-plots, eval=FALSE} 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. ```{r compare-gsea-daa, eval=FALSE} # 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.