Skip to main content

Functional Enrichment Workshop

Making sense of gene and proteins lists with functional enrichment analysis

13 FGSEA

fgsea (“fast gene set enrichment analysis”) is an R/Bioconductor package for running pre-ranked GSEA. Rather than testing a fixed list of “interesting” genes as ORA does, GSEA uses the ranks (or ranking statistic, e.g. log2 fold change or t-statistic) of all measured genes to ask whether a gene set is skewed towards the top or bottom of that ranking.

fgsea implements the same underlying pre-ranked GSEA method as the GSEA Desktop/MSigDB web tool you used earlier in this workshop, but as a fast, scriptable R function, using an adaptive multi-level Monte-Carlo scheme that lets it accurately estimate very small p-values without needing a large number of permutations. See the fgsea preprint for algorithmic detail, and the package vignette for the full reference this chapter is based on.

We will re-run pre-ranked GSEA against the MSigDB Hallmark gene sets, using the same pre-ranked gene list you used earlier in the GSEA Desktop tool (Pre_Ranked_List_logFC_SYMBOL), so we can compare the R output directly against what you saw on the web tool, including the HALLMARK_CHOLESTEROL_HOMEOSTASIS gene set discussed in the earlier challenge question.

 

In this chapter we will:

  1. Load the pre-ranked gene list used in the earlier GSEA Desktop activity
  2. Prepare the ranked statistics vector required by fgsea
  3. Retrieve the MSigDB Hallmark gene sets with msigdbr
  4. Run pre-ranked GSEA with the fgsea function
  5. Save the tabular results to a file
  6. Visualise an individual gene set with plotEnrichment
  7. Visualise many gene sets at once with plotGseaTable
  8. Collapse redundant/overlapping gene sets with collapsePathways
  9. Run an over-representation test with fora for comparison against ORA

All of the data used in this chapter is provided in the data/ folder of this workshop’s materials, so there’s nothing extra to download - just read on and run the code chunks below in order (in RStudio, open this same .Rmd file directly to run and edit the code yourself).

 

13.1 Load R packages

13.2 1. Load input data

Earlier, we used a pre-ranked gene list derived from the Pezzini et al. 2017 SH-SY5Y neuronal differentiation dataset: genes ranked by their log2 fold change between differentiated and undifferentiated cells (Pre_Ranked_List_logFC_SYMBOL).

We load that same file here so that our fgsea results are directly comparable to what you saw on the web tools.

The file has one row per gene, with Gene_Symbol and logFC columns, already sorted from most strongly up-regulated to most strongly down-regulated in the differentiated condition.

13.3 2. Prepare the ranked statistics vector

fgsea requires the ranking statistics as a single named numeric vector (gene symbol -> statistic), rather than a data frame. A small number of gene symbols in the table are duplicated (for example where multiple probes/transcripts map to the same gene symbol); for these we keep the value with the largest absolute log fold change.

13.4 3. Get MSigDB Hallmark gene sets

Earlier, the GSEA Desktop and STRING GSEA activities queried the MSigDB Hallmark gene set collection. The msigdbr package provides convenient, versioned, tidy access to MSigDB collections directly in R, without needing to download .gmt files by hand.

fgsea expects pathways as a named list of character vectors (gene set name -> genes), so we reshape the msigdbr table into that format:

13.5 4. Run fgsea

Before running the below code chunk, review the parameters for the fgsea function using the Help pane of RStudio (?fgsea::fgsea). Note the similarities to the parameters you saw in the GSEA Desktop tool, e.g. minSize/maxSize (gene set size filters).

fgsea has a default lower bound (eps = 1e-10) on the p-values it estimates, for speed. If a gene set could plausibly have an even smaller p-value than this, you will see a warning. To estimate p-values more precisely (at the cost of a slower run), set eps = 0:

The result table contains, per gene set: the enrichment score (ES), normalised enrichment score (NES), raw and adjusted p-values, and the leadingEdge genes driving the enrichment.

Let’s check whether HALLMARK_CHOLESTEROL_HOMEOSTASIS (discussed in the earlier challenge question) shows up as significantly enriched here too:

13.6 5. Save the results to a file

data.table::fwrite handles the list-column (leadingEdge) automatically, unlike write.csv/write.table.

13.7 6. Visualise a single gene set

plotEnrichment reproduces the classic GSEA “running enrichment score” plot for one gene set, the same style of plot produced by the GSEA Desktop tool:

13.8 7. Visualise the top gene sets as a table plot

plotGseaTable draws a compact summary of many gene sets at once, showing each gene set’s ranked-list position, running score trace, and NES/p-value.

13.9 8. Collapse redundant gene sets

Hallmark gene sets are curated to minimise overlap, but many gene set collections (e.g. GO, Reactome) contain highly overlapping terms. collapsePathways identifies a smaller, largely independent subset of the significant gene sets whose leading-edge genes are not already explained by another, more significant gene set:

13.10 9. Over-representation test with fora

fgsea also includes fora, a fast implementation of a classic hypergeometric over-representation test (like the ORA you ran earlier with g:Profiler and STRING). It requires a foreground gene list, a background/universe gene list, and a set of pathways.

We reuse the same DE_198set_Gene_Symbols.txt gene list you used for ORA in g:Profiler and STRING earlier as our foreground, and all genes measured in the experiment as our background, so we can compare the GSEA and ORA results for this dataset directly.

Question: Do the top gene sets from fora (ORA) match the top gene sets from fgsea (GSEA)? Why might they differ, given that ORA only sees the 198 significant genes while GSEA uses the ranking of all ~14,000 genes?

13.11 Performance considerations

fgsea is parallelised via the BiocParallel package. By default, the first registered backend returned by bpparam() is used (above we explicitly registered SerialParam() for reproducibility). To use multiple cores, set the nproc argument of fgsea(), or supply your own BPPARAM, e.g. BPPARAM = MulticoreParam(workers = 4).