--- 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 | |--------|-------------------------|------------|-------| | `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](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", 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 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. ```{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 = "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. ```{r depth-gsea, eval=FALSE} # 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. ```{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, 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. ```{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.