How to Run Gene Set Enrichment Analysis (GSEA) with fgsea in R: A Complete Workflow for RNA-seq

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-08-21

Key Takeaways

• The fgsea package in R implements a fast, preranked Gene Set Enrichment Analysis (GSEA) algorithm suitable for RNA-seq data, enabling the identification of statistically significant, concordant changes in predefined gene sets between biological states.
• A critical prerequisite is the preparation of a ranked gene list, typically derived from differential expression analysis, where genes are ordered by a robust metric such as the signed test statistic (balancing effect size and statistical reliability) and use consistent gene identifiers (e.g., gene symbols) matching the chosen gene set collection.
• Gene set collections, such as MSigDB Hallmark, Reactome, or KEGG, define the biological hypotheses; Hallmark sets are useful for exploratory analysis, while pathway databases offer granular mechanistic detail for hypothesis testing.
• fgsea execution requires careful parameter selection, including a sufficient number of permutations (e.g., 10,000 for exploratory, 100,000 for definitive results) to ensure accurate p-value estimation and the application of multiple testing correction (e.g., Benjamini-Hochberg) to control the false discovery rate.
• Interpretation relies on normalized enrichment scores (NES) and adjusted p-values, with leading edge genes providing insight into the specific gene set members driving the observed enrichment signal, which can be visualized through enrichment plots, bar charts, and heatmaps.
• Common analytical pitfalls include gene identifier mismatches between the ranked list and gene sets, duplicate gene names, and insufficient permutations, all of which can lead to erroneous or uninterpretable enrichment results.

Gene set enrichment analysis (GSEA) is a computational method that determines whether predefined sets of genes show statistically significant, concordant differences between two biological states. The fgsea package in R implements a fast preranked GSEA algorithm that is widely used for RNA-seq data because it runs quickly on genome-wide expression matrices and produces interpretable enrichment scores. This workflow covers the complete process from preparing a ranked gene list to running fgsea, interpreting results, and visualizing enriched pathways, with attention to multiple testing correction and common analytical pitfalls.

The intended reader is a biology student, researcher, or laboratory professional who has already generated RNA-seq count data and performed differential expression analysis. This article assumes basic familiarity with R syntax and the structure of gene expression matrices. The practical outcome is a reproducible fgsea pipeline that produces defensible biological conclusions from transcriptomic data.

At a Glance

RNA-seq Analysis Context and the Role of GSEA

RNA sequencing has become a standard approach for quantifying cell state across experimental conditions. The volume of public gene expression data continues to grow, and researchers routinely analyze newly generated datasets alongside existing repositories. The NCBI maintains major sequence databases and search systems that support the deposition and retrieval of RNA-seq data, while EMBL-EBI Training provides learning pathways for bioinformatics data resources and practical analysis education. These infrastructure resources underpin the reproducibility expectations for transcriptomic research.

Differential expression analysis identifies individual genes that change between conditions, but biological interpretation often requires understanding coordinated changes across pathways. A typical RNA-seq workflow includes quality control, read alignment, quantification, normalization, and differential expression testing. The Galaxy Training Network offers accessible workflow training and analysis tutorials that cover these steps, and nf-core documentation describes community pipeline standards for reproducible analysis. After differential expression produces a list of genes with test statistics, GSEA provides the next analytical layer by asking whether known biological pathways are enriched among the most changed genes.

The core biological rationale for GSEA is that subtle but coordinated expression changes across many genes in a pathway may be biologically meaningful even when individual genes do not pass significance thresholds. This is particularly relevant for complex phenotypes where no single gene drives the response. The method ranks all genes by a metric such as fold change or signal-to-noise ratio, then walks down the ranked list to calculate an enrichment score that reflects whether members of a given gene set cluster at the top or bottom of the ranking.

Core Principles of fgsea

The fgsea package implements an efficient algorithm for preranked GSEA. The input is a ranked list of genes, and the output is an enrichment score for each gene set tested. The algorithm uses an adaptive permutation scheme that reduces the computational burden compared to the original GSEA implementation while maintaining statistical accuracy.

The enrichment score reflects the degree to which a gene set is overrepresented at the extremes of the ranked list. A positive enrichment score indicates that genes in the set are concentrated at the top of the ranking, while a negative score indicates concentration at the bottom. The normalized enrichment score accounts for gene set size differences, allowing comparison across pathways.

The fgsea package is distributed through Bioconductor, which provides official package documentation, installation instructions, and reproducible genomic analysis workflows. Installing fgsea requires R version 3.5 or higher and uses the standard Bioconductor installation command. The package depends on several other Bioconductor packages for data structures and parallel computing.

Multiple testing correction is essential when testing thousands of gene sets simultaneously. The fgsea output includes raw p-values, and researchers should apply the Benjamini-Hochberg procedure to control the false discovery rate. Pathways with adjusted p-values below a threshold such as 0.05 or 0.1 are typically reported as significantly enriched.

Preparing the Ranked Gene List

The ranked gene list is the most important input for fgsea. The quality of the ranking directly determines the quality of the enrichment results. A poorly constructed ranking can produce misleading enrichment scores even when the underlying differential expression analysis is sound.

Choosing a Ranking Metric

The ranking metric should reflect the biological signal of interest. Common choices include log2 fold change, the signed test statistic from differential expression analysis, or the product of the sign of fold change and the negative logarithm of the p-value. The log2 fold change alone is simple and interpretable, but it does not account for the statistical confidence of each measurement. The test statistic from tools like DESeq2 or limma incorporates both effect size and variance, making it a more robust ranking metric.

For most RNA-seq applications, ranking by the signed test statistic is recommended because it balances effect size with statistical reliability. Genes with large fold changes but high variance will rank lower than genes with moderate fold changes and low variance, which is appropriate for identifying robust biological signals.

Handling Gene Identifiers

The ranked vector must use gene identifiers that match the gene set collection. If the differential expression results use Ensembl gene IDs but the gene sets use gene symbols, the identifiers must be converted before running fgsea. The Bioconductor project provides annotation packages for identifier conversion, and the The Carpentries lessons offer foundational programming training that covers data manipulation in R.

Duplicate gene identifiers in the ranked list will cause errors or produce ambiguous results. Each gene should appear exactly once in the ranked vector. When multiple probes or transcripts map to the same gene, the researcher must decide whether to keep the maximum value, the mean value, or the most significant entry. The choice depends on the data structure and should be documented in the analysis methods.

Sorting and Formatting

The ranked vector must be sorted in decreasing order, with the most upregulated gene first and the most downregulated gene last. The vector names must be the gene identifiers, and the values must be the ranking metric. A typical ranked vector is created from a differential expression results table using the following logic:

ranks <- results$stat
names(ranks) <- results$gene_symbol
ranks <- sort(ranks, decreasing = TRUE)

The resulting vector should be inspected for missing values, duplicate names, and unexpected sorting behavior before proceeding to the fgsea step.

Gene Set Collections and Databases

Gene set collections define the biological hypotheses being tested. The choice of collection shapes the interpretation of results and should be guided by the biological question.

Molecular Signatures Database

The Molecular Signatures Database (MSigDB) is a widely used collection of annotated gene sets maintained by the Broad Institute. It includes hallmark gene sets that summarize specific well-defined biological states, curated canonical pathways from databases such as Reactome and KEGG, and regulatory target gene sets. Hallmark gene sets are often the first choice for exploratory analysis because they reduce noise and provide broad biological coverage.

Access to MSigDB requires registration, and gene sets are downloaded as GMT files. The GMT format is a tab-delimited text format where each line represents one gene set with a name, description, and member genes. The fgsea package includes functions to read GMT files directly.

Gene Ontology and Pathway Databases

Gene Ontology (GO) terms provide a controlled vocabulary for biological processes, molecular functions, and cellular components. GO gene sets can be downloaded from the Gene Ontology Consortium or accessed through Bioconductor annotation packages. Pathway databases such as Reactome and KEGG provide curated pathway definitions that are useful for interpreting mechanism-level biology.

The EMBL-EBI Training resources describe how to access and use these biological databases for practical analysis education. The choice between GO and pathway databases depends on the research question. GO terms are broader and more numerous, while pathway databases provide more mechanistic detail.

Species-Specific Considerations

Gene sets are typically defined using human gene symbols, and applying them to other species requires ortholog mapping. Many gene set collections provide mouse versions, and some provide other model organisms. When working with non-model organisms, the researcher must map gene identifiers to orthologous human genes or use species-specific gene set collections if available.

The NCBI provides ortholog resources and gene information that support cross-species analysis. The mapping step should be documented carefully because incomplete ortholog mapping can reduce the effective size of gene sets and affect enrichment statistics.

Running fgsea

The fgsea function takes the ranked gene list and the gene set collection as primary inputs. Additional parameters control the minimum and maximum gene set size, the number of permutations, and the random seed for reproducibility.

Basic Function Call

The basic fgsea call uses the following structure:

fgsea_results <- fgsea(pathways = gene_sets,
                       stats = ranks,
                       minSize = 15,
                       maxSize = 500,
                       nperm = 10000)

The minSize parameter removes gene sets with fewer than the specified number of genes after the intersection with the ranked list. Small gene sets produce unstable enrichment scores because the statistic is based on few observations. The maxSize parameter removes very large gene sets that may represent broad, nonspecific categories. Common choices are a minimum of 15 genes and a maximum of 500 genes.

Permutation Count and Computational Cost

The number of permutations determines the precision of the p-value estimates. More permutations produce more stable p-values but require more computation time. The fgsea algorithm uses an adaptive approach that performs fewer permutations for gene sets with extreme enrichment scores, which reduces the overall computational burden.

For exploratory analysis, 10,000 permutations provide a reasonable balance between speed and accuracy. For final results that will be reported in publications, 100,000 permutations may be appropriate. The computational cost depends on the number of gene sets tested and the size of the ranked list.

Parallel Processing

The fgsea package supports parallel processing through the BiocParallel framework. Setting up multiple cores can substantially reduce runtime for large analyses. The Bioconductor documentation describes how to configure parallel backends for genomic analysis workflows.

library(BiocParallel)
register(MulticoreParam(4))

The parallel configuration should be documented in the analysis methods for reproducibility.

Interpreting fgsea Results

The fgsea output is a data frame with one row per gene set. Key columns include the pathway name, the enrichment score, the normalized enrichment score, the raw p-value, the adjusted p-value, and the leading edge genes.

Enrichment Score and Normalized Enrichment Score

The enrichment score (ES) reflects the maximum deviation of the running sum statistic from zero as the algorithm walks down the ranked list. The normalized enrichment score (NES) adjusts for gene set size and allows comparison between gene sets. Positive NES values indicate enrichment at the top of the ranked list, and negative values indicate enrichment at the bottom.

The magnitude of the NES provides a measure of effect size. Larger absolute values indicate stronger enrichment. However, the NES should always be interpreted alongside the p-value because a large NES based on few genes may not be statistically significant.

P-values and Multiple Testing Correction

The raw p-value from fgsea estimates the probability of observing an enrichment score at least as extreme as the observed score under the null hypothesis. Because thousands of gene sets are tested simultaneously, the raw p-values must be adjusted for multiple testing.

The Benjamini-Hochberg procedure controls the false discovery rate and is the standard approach for GSEA results. The adjusted p-value is calculated as follows:

fgsea_results$padj <- p.adjust(fgsea_results$pval, method = "BH")

Pathways with adjusted p-values below the chosen threshold are considered significantly enriched. A common threshold is 0.05, but 0.1 may be appropriate for exploratory analyses where the goal is hypothesis generation.

Leading Edge Genes

The leading edge comprises the genes in a gene set that contribute most to the enrichment score. These genes are the core members of the pathway that drive the enrichment signal. Examining the leading edge genes can reveal which specific biological processes are most affected in the comparison.

The leading edge information is particularly useful for interpreting results in the context of the experimental system. For example, in a study of peripheral artery disease muscle microenvironments, the leading edge genes from enriched pathways identified specific endothelial cell and macrophage programs that were altered in disease. Similarly, single-cell RNA dynamics during Salmonella infection revealed that RNA synthesis and degradation dynamics shape the interpretation of total RNA levels, which has implications for interpreting enrichment results from bulk RNA-seq data.

Visualizing Enrichment Results

Visualization is essential for communicating GSEA results to collaborators and readers. The fgsea package includes plotting functions, and additional visualization options are available through base R graphics and the ggplot2 package.

Enrichment Plots

The enrichment plot shows the running enrichment score across the ranked gene list for a single gene set. The plot includes a barcode-like representation of where the gene set members appear in the ranking and a curve showing the running enrichment score. The peak of the curve corresponds to the enrichment score.

plotEnrichment(gene_sets[["HALLMARK_WNT_BETA_CATENIN_SIGNALING"]], ranks)

This plot is the standard figure for reporting individual pathway enrichment. It should be generated for the top enriched pathways and included in supplementary materials.

Bar Charts and Dot Plots

Bar charts of normalized enrichment scores for the top pathways provide a quick overview of the results. The pathways can be ordered by adjusted p-value or by NES, and colored by the direction of enrichment. Dot plots that combine the NES, adjusted p-value, and gene set size on a single figure are also effective for summarizing many pathways at once.

The Phantasus web application provides an interactive environment for gene expression analysis that includes heatmap visualization and downstream analysis capabilities. While fgsea is typically run from the R command line, interactive tools can complement the workflow for exploratory visualization.

Heatmaps of Leading Edge Genes

Heatmaps showing the expression of leading edge genes across samples can validate the enrichment results and provide a visual check on the direction of change. The Methods in Molecular Biology chapter on RNA-seq heat maps describes a protocol for generating heat maps in R that is accessible to users without prior R experience.

The heatmap should use the normalized expression values from the RNA-seq analysis and should be annotated with sample groups. The leading edge genes should show coordinated expression differences between the conditions being compared.

Handling Multiple Testing and Statistical Reporting

The statistical reporting for GSEA results must be transparent and complete. The methods section of a paper or report should describe the ranking metric, the gene set collection, the parameters used for fgsea, and the multiple testing correction approach.

Reporting Standards

The following information should be reported for every GSEA analysis:

1. The version of R and fgsea used
2. The differential expression tool and version used to generate the ranking
3. The ranking metric and how it was calculated
4. The gene set collection and version or download date
5. The minimum and maximum gene set sizes
6. The number of permutations
7. The multiple testing correction method
8. The significance threshold

This level of detail supports reproducibility and allows other researchers to evaluate the analysis. The nf-core documentation emphasizes reproducibility standards for bioinformatics pipelines, and the Galaxy Training Network provides tutorials that model transparent reporting practices.

Common Statistical Pitfalls

A common error is reporting raw p-values without multiple testing correction. With thousands of gene sets tested, raw p-values below 0.05 are expected by chance alone. The adjusted p-value is the appropriate statistic for claiming significance.

Another pitfall is interpreting the NES without considering the p-value. A pathway with a large NES but a high adjusted p-value should not be reported as significantly enriched. The NES and adjusted p-value should be reported together.

The permutation count also affects the precision of p-values. With too few permutations, the minimum possible p-value is limited by the number of permutations. For example, with 1,000 permutations, the minimum p-value is approximately 0.001, which may not be sufficiently precise for highly significant pathways.

Common Failure Patterns and Troubleshooting

Several recurring problems arise when running fgsea on RNA-seq data. Recognizing these patterns and knowing how to address them saves time and prevents incorrect biological conclusions.

Identifier Mismatches

The most common failure is a mismatch between the gene identifiers in the ranked list and the gene set collection. If the ranked list uses Ensembl IDs and the gene sets use gene symbols, the intersection will be small or empty, and fgsea will report that most gene sets have zero genes after filtering.

The solution is to convert identifiers before running fgsea. The Bioconductor annotation packages provide mapping tables, and the conversion should be verified by checking that a reasonable fraction of genes in the ranked list appear in the gene sets.

Duplicate Gene Names

Duplicate gene names in the ranked vector cause fgsea to fail or produce incorrect results. The duplicated() function in R can identify duplicates, and the researcher must decide how to resolve them before running the analysis.

A common approach is to keep the entry with the largest absolute ranking metric for each gene. This preserves the most extreme observation and avoids arbitrary selection.

Gene Set Size Filtering

Gene sets that are too small or too large after intersection with the ranked list are excluded by the minSize and maxSize parameters. If many gene sets are excluded, the researcher should check whether the identifier conversion worked correctly and whether the ranked list has sufficient genome-wide coverage.

Memory Limitations

Running fgsea on very large gene set collections can require substantial memory. The nperm parameter and the number of gene sets both affect memory usage. Reducing the number of permutations or testing a subset of gene sets can address memory limitations.

Quality Control and Reproducibility

Quality control for GSEA extends beyond the fgsea run itself. The input data must be checked for quality before the analysis, and the results should be validated after the analysis.

Input Data Quality Checks

The differential expression results should be inspected for the distribution of p-values and fold changes. An excess of very small p-values may indicate problems with the statistical model or with sample labeling. The single-cell RNA-seq best practices tutorial describes quality control steps that are relevant to RNA-seq analysis generally, including normalization and data correction.

The ranked list should be checked for the number of genes and the distribution of ranking metrics. A ranked list with very few genes will produce unstable enrichment scores because the gene set size filters will exclude most pathways.

Batch Effects and Confounding

Batch effects can distort differential expression results and consequently affect GSEA results. The Caris-ComBat-seq method addresses the challenge of harmonizing large-scale RNA-seq datasets from different platforms, which is relevant when combining data from multiple sources. If the RNA-seq data come from different batches or platforms, batch correction should be performed before differential expression analysis.

Confounding variables such as sex, age, or tissue source should be included in the differential expression model. The variant calling workflow for RNA-seq data demonstrates how RNA-seq data can be used for genomic analysis, but it also highlights the importance of understanding the data generation process.

Reproducibility Practices

The complete analysis should be captured in a script or notebook that can be rerun from the raw data. The The Carpentries lessons provide training on reproducible computing practices, including version control and project organization.

Setting a random seed before running fgsea ensures that the permutation-based p-values are reproducible. The seed should be recorded in the analysis script.

set.seed(42)

The R session information should be saved using sessionInfo() and included in the analysis documentation.

Limitations of fgsea and GSEA

GSEA is a powerful method, but it has limitations that researchers must understand when interpreting results.

Gene Set Composition and Redundancy

Gene sets overlap substantially, particularly within collections like GO terms. A single biological process may be represented by dozens of overlapping gene sets, and the enrichment results will be redundant. This redundancy complicates interpretation because multiple significant gene sets may reflect the same underlying biological signal.

Clustering or summarizing redundant gene sets can help identify the core biological themes. The leading edge genes can be used to identify the shared genes driving enrichment across multiple gene sets.

Sensitivity to Ranking Metric

The choice of ranking metric affects the results. Different metrics may produce different rankings and therefore different enrichment scores. The ranking metric should be chosen based on the biological question and the statistical properties of the differential expression results.

The NF-κB signaling study in pancreatic cancer demonstrates how transcriptome-wide profiling combined with motif analysis can reveal distinct regulatory programs. The choice of ranking metric in such analyses should reflect whether the research question emphasizes effect size, statistical significance, or a combination of both.

Tissue and Cell Type Heterogeneity

Bulk RNA-seq measures the average expression across all cells in a sample. If the tissue contains multiple cell types with different expression programs, the enrichment results reflect the dominant signal instead of cell-type-specific changes. The single-cell compendium of muscle microenvironment in peripheral artery disease illustrates how single-cell analysis can reveal cell-type-specific programs that are masked in bulk data.

When interpreting GSEA results from bulk RNA-seq, the researcher should consider whether the enriched pathways are consistent with the known cell type composition of the tissue. If cell-type-specific effects are of interest, single-cell RNA-seq may be a more appropriate approach.

Pathway Database Coverage

The gene set collection determines which biological processes can be detected. If a relevant pathway is not represented in the collection, the analysis will not detect its enrichment. The choice of gene set collection should be guided by the biological question and should be reported transparently.

The hepatoblastoma multiomic analysis identified WNT signaling pathway enrichment in Beckwith-Wiedemann syndrome tumors, demonstrating how pathway-level analysis can reveal key biological programs. The choice of gene set collection in that study was informed by the known biology of the disease.

Practical Workflow Summary

The following steps summarize the complete fgsea workflow from differential expression results to interpreted biological conclusions.

Step 1: Prepare the Ranked Gene List

Extract the ranking metric from the differential expression results table. Convert gene identifiers to match the gene set collection. Remove duplicate genes and missing values. Sort the vector in decreasing order.

Step 2: Load and Filter Gene Sets

Read the gene set collection from a GMT file or load it from a Bioconductor annotation package. Filter gene sets by size if necessary. Verify that the gene identifiers match the ranked list.

Step 3: Run fgsea

Set the random seed for reproducibility. Configure parallel processing if available. Run the fgsea function with appropriate parameters for the analysis.

Step 4: Apply Multiple Testing Correction

Calculate adjusted p-values using the Benjamini-Hochberg method. Filter the results to retain gene sets that meet the significance threshold.

Step 5: Visualize and Interpret Results

Generate enrichment plots for the top pathways. Create summary visualizations such as bar charts or dot plots. Examine the leading edge genes to understand the biological drivers of enrichment.

Step 6: Document and Report

Record all parameters, versions, and settings. Save the results table and visualization outputs. Write a methods section that describes the analysis completely.

Records and Measurements

Maintaining detailed records of the GSEA analysis supports reproducibility and enables troubleshooting when results are unexpected.

Analysis Log

The analysis log should record the date, the R version, the fgsea version, the input file names, and all parameter settings. This information should be saved alongside the results.

Parameter Documentation

The following parameters should be documented for every fgsea run:

Output File Management

The fgsea results table should be saved as a CSV or TSV file for downstream analysis and sharing. The enrichment plots should be saved as high-resolution image files for publication.

Professional Escalation Criteria

Some analytical situations require consultation with a bioinformatics specialist or statistician. Recognizing these situations prevents incorrect biological conclusions.

When to Seek Specialized Help

Consult a bioinformatics specialist when the differential expression results show unusual patterns, such as an excessive number of genes with very small p-values or a strong inflation of test statistics. These patterns may indicate problems with the statistical model, sample contamination, or batch effects that require specialized expertise.

Seek statistical consultation when the gene set collection is large and the multiple testing burden is substantial. The interpretation of adjusted p-values in the context of highly correlated gene sets requires statistical judgment.

Data Quality Red Flags

The following observations warrant professional review:

1. Very few genes pass the gene set size filters after identifier conversion
2. The enrichment scores are consistently extreme across all gene sets
3. The results change dramatically when the permutation count is increased
4. The leading edge genes for most significant pathways are the same genes

These patterns suggest technical problems with the input data or the analysis configuration.

Building a Decision Framework for Gene Set Selection and Result Prioritization

The choice of gene set collection and the interpretation of enrichment results are the two points in the fgsea workflow where biological judgment matters most. Researchers often default to the MSigDB Hallmark collection because it is convenient, but this convenience can obscure biologically relevant signals that reside in more granular pathway definitions. A structured decision framework helps match the gene set collection to the specific research question and prevents the common failure of reporting redundant or biologically uninformative results.

Matching Gene Set Collections to Research Questions

The Hallmark collection contains 50 gene sets that summarize well-defined biological states by merging multiple source databases. These sets are designed to reduce noise and provide broad coverage, making them appropriate for initial exploratory analysis when the biological hypothesis is not yet refined. However, the Hallmark collection has a critical limitation: it collapses distinct biological processes into single broad categories. For example, the Hallmark WNT beta catenin signaling set combines genes from multiple WNT pathway branches, which may obscure the specific WNT ligand, receptor, or downstream effector that drives the observed expression changes.

For mechanistic hypothesis testing, curated pathway databases such as Reactome or KEGG provide more granular gene sets that separate pathway branches and individual signaling cascades. The EMBL-EBI Training resources describe how to access these pathway databases and explain the differences in gene set granularity across collections. A practical approach is to run fgsea with the Hallmark collection first to identify the broad biological themes, then follow up with a more granular collection such as Reactome to dissect the specific pathway components driving the enrichment signal.

The hepatoblastoma multiomic analysis illustrates this principle in practice. The study identified WNT signaling pathway enrichment in Beckwith-Wiedemann syndrome tumors using pathway-level analysis, but the authors did not stop at the broad pathway annotation. They used pseudotime analysis and single-cell profiling to identify a specific population of transition cells with unique molecular profiles that likely drive the precancer to cancer neoplastic transition. This layered approach, moving from broad pathway enrichment to specific cellular programs, is the model for translating GSEA results into mechanistic insight.

A Three-Tier Gene Set Selection Framework

A practical decision framework for gene set selection uses three tiers based on the research objective and the stage of the analysis.

Tier 1: Exploratory screening. Use the Hallmark collection or a similarly curated broad collection when the goal is to identify which biological processes are altered in the comparison. This tier is appropriate for initial analysis of a new dataset, for hypothesis generation, and for confirming that the experimental system behaves as expected. The output is a short list of broad biological themes that guide subsequent analysis.

Tier 2: Mechanistic dissection. Use curated pathway databases such as Reactome, KEGG, or WikiPathways when the goal is to understand which specific pathway components drive the enrichment signal. This tier is appropriate when the Tier 1 analysis has identified a broad theme and the researcher needs to know which ligands, receptors, or downstream effectors are involved. The output is a more detailed picture of the pathway architecture and the specific nodes that are altered.

Tier 3: Regulatory and functional annotation. Use Gene Ontology biological process terms, transcription factor target sets, or microRNA target sets when the goal is to understand the regulatory logic underlying the expression changes. This tier is appropriate when the research question concerns upstream regulators, downstream functional consequences, or the coordination of multiple biological processes. The output is a set of regulatory hypotheses that can be tested experimentally.

The NF-kB signaling study in pancreatic cancer demonstrates the value of this tiered approach. The authors used transcriptome-wide gene expression profiling to identify differential activation of canonical and noncanonical NF-kB signaling, then used motif enrichment and chromatin accessibility analysis to determine transcription factor binding dynamics. The Tier 3 analysis revealed that RELB exclusively occupies pre-accessible chromatin regions co-enriched for AP1 motifs, providing a mechanistic explanation for the distinct regulatory roles of RELA and RELB that would not have emerged from a Tier 1 analysis alone.

Prioritizing Results Across Multiple Gene Set Collections

When running fgsea with multiple gene set collections, the results must be integrated into a coherent biological narrative. A common failure pattern is reporting the same biological theme multiple times across different collections without recognizing the redundancy. The following prioritization framework addresses this problem.

First, identify the core biological themes by clustering the significant gene sets based on their leading edge gene overlap. Gene sets that share a high proportion of leading edge genes are likely reflecting the same underlying biological signal. The single-cell RNA-seq best practices tutorial describes clustering approaches that can be adapted for gene set results, and the The Carpentries lessons provide foundational programming training for implementing these analyses in R.

Second, for each cluster of redundant gene sets, select the representative gene set that has the lowest adjusted p-value, the largest absolute normalized enrichment score, and the most specific biological annotation. A specific annotation is one that names a particular process, pathway, or molecular function instead of a broad category. For example, a gene set named "WNT ligand biogenesis and trafficking" is more specific than one named "WNT signaling pathway" and should be preferred when both are significant.

Third, order the prioritized gene sets by biological theme instead of by statistical significance alone. This ordering helps the reader understand the major biological programs that are altered in the comparison and prevents the common failure of presenting a list of statistically significant but biologically redundant pathways.

A Record System for Gene Set Selection Decisions

The gene set selection decisions should be documented in the analysis record alongside the fgsea parameters. The following table format captures the relevant information for each gene set collection used in the analysis.

The rationale column should state the biological question that motivated the selection. For example, "Hallmark collection selected to identify broad biological themes in the comparison of treated versus untreated samples" or "Reactome selected to dissect the specific components of the inflammatory response identified in the Tier 1 analysis."

Common Failure Patterns in Gene Set Selection

Several recurring problems arise from poor gene set selection decisions. Recognizing these patterns helps researchers avoid them.

Collection mismatch with biological question. Using the Hallmark collection when the research question concerns a specific pathway branch will produce results that are too broad to be mechanistically informative. The solution is to match the collection granularity to the question, using the three-tier framework as a guide.

Uncritical acceptance of default collections. The default gene set collection in many analysis pipelines is the Hallmark collection, but this default may not be appropriate for all research questions. The Galaxy Training Network provides tutorials that demonstrate how to select and configure gene set collections for different analysis goals, and the nf-core documentation describes how community pipelines handle gene set selection in reproducible workflows.

Ignoring species-specific considerations. Gene set collections are typically defined using human gene symbols, and applying them to other species requires ortholog mapping. The NCBI provides ortholog resources that support cross-species analysis, but the mapping step should be documented carefully because incomplete ortholog mapping can reduce the effective size of gene sets and affect enrichment statistics.

Overlooking gene set redundancy. Reporting all significant gene sets without addressing redundancy produces a results section that is difficult to interpret and may overstate the number of distinct biological findings. The prioritization framework described above addresses this problem by clustering redundant gene sets and selecting representative members.

Validation of Prioritized Results

The prioritized gene sets should be validated before they are reported as biological findings. Validation serves two purposes: confirming that the enrichment signal is robust to analysis parameters and confirming that the biological interpretation is consistent with the experimental system.

Parameter sensitivity analysis. The fgsea results should be checked for sensitivity to the ranking metric and the permutation count. If the top enriched pathways change dramatically when the ranking metric is changed from log2 fold change to the test statistic, the enrichment signal is not robust and should be interpreted with caution. Similarly, if the adjusted p-values change substantially when the permutation count is increased from 10,000 to 100,000, the p-value estimates are not stable.

Biological consistency check. The prioritized gene sets should be consistent with the known biology of the experimental system. For example, in a study of peripheral artery disease muscle microenvironments, the enriched pathways should include angiogenesis, immune cell activation, and extracellular matrix remodeling, which are known to be altered in ischemic muscle. If the top enriched pathways are unrelated to the known biology of the system, the analysis should be reviewed for technical problems.

Cross-validation with complementary data. When possible, the enrichment results should be cross-validated with complementary data types. The single-cell compendium of muscle microenvironment in peripheral artery disease used single-cell RNA sequencing and spatial transcriptomics to validate the bulk RNA-seq findings, revealing cell-type-specific programs that were masked in the bulk data. The scIVNL-seq study of Salmonella infection demonstrated that RNA synthesis and degradation dynamics can shape the interpretation of total RNA levels, which has implications for validating enrichment results from bulk RNA-seq data.

Escalation Criteria for Gene Set Selection Problems

Some gene set selection problems require consultation with a bioinformatics specialist or a domain expert in the relevant biological pathway. The following situations warrant escalation.

Persistent lack of significant enrichment. If no pathways are significantly enriched across multiple gene set collections and parameter settings, the problem may lie in the differential expression analysis instead of the gene set selection. A bioinformatics specialist can review the differential expression model and the quality of the ranked list.

Conflicting results across collections. If different gene set collections produce conflicting biological interpretations, the conflict may indicate that the enrichment signal is driven by a small number of genes that appear in multiple collections. A domain expert can help determine which interpretation is biologically plausible.

Unexpected enrichment of unrelated pathways. If the top enriched pathways are consistently unrelated to the experimental system, the analysis should be reviewed for technical problems such as sample mislabeling, batch effects, or identifier mapping errors. The Caris-ComBat-seq method addresses the challenge of harmonizing large-scale RNA-seq datasets from different platforms, which is relevant when combining data from multiple sources.

Practical Implementation Steps

The following steps implement the decision framework in a typical fgsea analysis.

Step 1: Define the biological question. Write a one-sentence statement of the biological question that the GSEA analysis will address. This statement guides the gene set selection and the interpretation of results.

Step 2: Select the gene set collection tier. Use the three-tier framework to select the appropriate collection or collections. Document the selection rationale in the analysis record.

Step 3: Run fgsea with the selected collections. Use the same ranked gene list and fgsea parameters for all collections to ensure comparability of results.

Step 4: Cluster and prioritize significant gene sets. Identify redundant gene sets by leading edge gene overlap and select representative members for each biological theme.

Step 5: Validate the prioritized results. Check parameter sensitivity and biological consistency. Cross-validate with complementary data when available.

Step 6: Document the selection and prioritization decisions. Record the collection versions, the tier assignments, the rationale for each selection, and the prioritization results in the analysis record.

The Phantasus web application provides an interactive environment for gene expression analysis that can complement this workflow. While fgsea is typically run from the R command line, interactive tools can help researchers explore the results and make informed decisions about gene set selection and prioritization. The application provides streamlined access to public gene expression datasets and allows analysis of user-uploaded datasets, which can be useful for comparing the enrichment results against public data from similar experiments.

Frequently Asked Questions

What is the difference between GSEA and over-representation analysis?

GSEA uses all measured genes ranked by a statistic and tests whether gene set members are concentrated at the extremes of the ranking. Over-representation analysis tests whether genes that pass a significance threshold are enriched in a gene set. GSEA does not require a significance threshold and can detect coordinated changes in genes that individually do not reach significance.

How many permutations should I use for fgsea?

The number of permutations determines the precision of the p-value estimates. For exploratory analysis, 10,000 permutations provide a reasonable balance between speed and accuracy. For final results, 100,000 permutations may be appropriate. The minimum possible p-value is approximately one divided by the number of permutations.

What is the leading edge in GSEA results?

The leading edge comprises the genes in a gene set that contribute most to the enrichment score. These genes appear before the peak of the running enrichment score and represent the core members of the pathway driving the enrichment signal. Examining the leading edge genes helps interpret which specific biological processes are most affected.

How do I convert gene identifiers for fgsea?

Gene identifier conversion can be performed using Bioconductor annotation packages such as org.Hs.eg.db for human data. The conversion maps between identifiers such as Ensembl gene IDs, Entrez IDs, and gene symbols. After conversion, verify that a reasonable fraction of genes in the ranked list appear in the gene set collection.

Can I run fgsea on data from non-human species?

Yes, but the gene set collection must match the species. Many collections provide mouse versions, and some provide other model organisms. For non-model organisms, gene identifiers must be mapped to orthologous human genes or species-specific gene sets must be used. The mapping step should be documented carefully.

What is the difference between enrichment score and normalized enrichment score?

The enrichment score reflects the maximum deviation of the running sum statistic from zero. The normalized enrichment score adjusts for gene set size, allowing comparison between gene sets of different sizes. The normalized enrichment score is the appropriate statistic for comparing enrichment across pathways.

How should I report fgsea results in a publication?

Report the ranking metric, the gene set collection and version, the fgsea parameters including minimum and maximum gene set size and permutation count, the multiple testing correction method, and the significance threshold. Include the R and fgsea versions and the random seed for reproducibility.

What should I do if no pathways are significantly enriched?

First check the quality of the ranked list and the gene set collection. Verify that the gene identifiers match and that the gene set sizes are appropriate. Consider whether the biological contrast is strong enough to produce detectable enrichment. If the differential expression signal is weak, the enrichment analysis may not have sufficient power.

Related Bioinformatics Guides

• Gene Set Enrichment Analysis in R: A Practical Tutorial for Interpreting Omics Data
• RNA-Seq Data Analysis Workflow: From Raw Reads to Insights
• How to Interpret Gene Set Enrichment Analysis Results
• Gene Set Enrichment Analysis Tools: Choosing the Right One
• RNA-Seq Data Analysis in Galaxy: A User-Friendly Platform

Related Clinical & Scientific Guides

• A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data
• Computational Immunology: Modeling the Immune System
• How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices

References and Further Reading

• NCBI Data Resources. National Center for Biotechnology Information.
• EMBL-EBI Training. European Bioinformatics Institute.
• Bioconductor. Bioconductor Project.
• Galaxy Training Network. Galaxy Project.
• nf-core Documentation. nf-core.
• The Carpentries Lessons. The Carpentries.
• A DNA Methylation-based Computational Framework for Tumour-microenvironment State Inference and Molecular Stratification. 2026.
• Beckwith-Wiedemann syndrome multiomic analysis of hepatoblastoma uncovers unique tumour heterogeneity and cellular landscapes, including transition cells leading to tumour formation.. 2026.
• Epigenetic context defines the transcriptional activity of canonical and noncanonical NF-κB signaling in pancreatic cancer.. 2026.
• Phantasus, a web application for visual and interactive gene expression analysis.. 2024.
• Single-cell compendium of muscle microenvironment in peripheral artery disease reveals altered endothelial diversity and LYVE1⁺ macrophage activation.. 2025.
• scIVNL-seq resolves in vivo single-cell RNA dynamics of immune cells during Salmonella infection.. 2025.
• Current best practices in single-cell RNA-seq analysis: a tutorial. Molecular Systems Biology, 2019.
• Tutorial: RNA-seq differential expression & pathway analysis with Sailfish, DESeq2, GAGE, and Pathview. 2015.
• Galaxy Training Material for single-cell RNA-seq tutorial with Plant Datasets. 2021.
• Become Competent in Generating RNA-Seq Heat Maps in One Day for Novices Without Prior R Experience.. Methods in molecular biology, 2021.
• Variant Calling from RNA-seq Data Using the GATK Joint Genotyping Workflow.. Methods in molecular biology, 2022.
• Caris-ComBat-seq: Directionally Harmonizing Large-Scale RNA-seq Datasets. bioRxiv, 2025.

This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.