Skip to main content

Functional Enrichment Workshop

Making sense of gene and proteins lists with functional enrichment analysis

11 FGSEA

fgsea (“fast gene set enrichment analysis”) is an R package for running pre-ranked GSEA. 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 in the previous chapter, but as a fast, scriptable R function. See the fgsea preprint for algorithmic detail, and the package vignette for the full reference this chapter is based on.

11.1 Input data

We will re-run pre-ranked GSEA against the MSigDB Hallmark gene sets, using the same pre-ranked gene list you used earlier (Pre_Ranked_List_logFC_SYMBOL), so we can compare the R output directly against what you saw in the previous chapter, including the HALLMARK_CHOLESTEROL_HOMEOSTASIS gene set discussed in the earlier challenge question.
All of the data used in this chapter is provided in the data/R_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.

 

11.2 Activity overview

  1. Load the pre-ranked gene list
  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

11.3 Load R packages 

library(readr)
library(dplyr)
library(fgsea)
library(data.table)
library(ggplot2)
library(msigdbr)
library(BiocParallel)

register(SerialParam())

11.4 Load input data

We load the DE genes table derived from the Pezzini et al. 2017 SH-SY5Y neuronal differentiation dataset. The file has one row per gene, with Gene_Symbol and logFC columns and other related statistics columns, comparing differentiated against undifferentiated cells (Pre_Ranked_List_logFC_SYMBOL).

stats_df <- read.table("data/R_data/Pezzini_DE.txt", header = TRUE)
head(stats_df)

11.5 Prepare the ranked statistics vector

fgsea requires the ranking statistics as a single named numeric vector (gene id -> statistic), so we extract the log2 fold change column as a named vector, and arrange it by decreasing log2 fold change. Here, genes at the start of this vector reflect a strong increase in the differentiated cells compared to the undifferentiated cells, while genes at the bottom of this list are strongly decreased in the differentiated cells.

stats <- setNames(stats_df$Log2FC, stats_df$Gene.ID)

# fgsea expects the vector sorted, largest log2 fold change to smallest
stats <- sort(stats, decreasing = TRUE)

11.6 Get MSigDB Hallmark gene sets

Earlier, we 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. There are 50 unique Hallmark gene sets in this collection. After running the code below, have a look at what they are by running unique(hallmark_df$gs_name) in your console.

hallmark_df <- msigdbr(species = "Homo sapiens", collection = "H")
head(hallmark_df)

fgsea expects pathways as a named list of character vectors (gene set name -> genes), so we reshape the msigdbr table into that format. Here we use the ensembl gene IDs, as some of the gene symbols in the input table are duplicated.

pathways <- split(hallmark_df$ensembl_gene, hallmark_df$gs_name)
str(head(pathways))

11.7 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::fgsea

Let’s run fgsea:

set.seed(42)

fgseaRes <- fgsea(pathways = pathways,
                  stats    = stats, #   Named vector of gene-level stats. Names should be the same as in 'pathways'
                  minSize  = 15, #  Minimal size of a gene set to test. All pathways below the threshold are excluded.
                  maxSize  = 500, # Maximal size of a gene set to test. All pathways above the threshold are excluded.
                  eps = 1e-10)

#`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.

head(fgseaRes[order(pval), ])

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

fgseaRes[pathway == "HALLMARK_CHOLESTEROL_HOMEOSTASIS"]

11.8 Save the results to a file

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

title <- "fgsea_Hallmark_results.txt"
fwrite(fgseaRes, file = title, sep = "\t", sep2 = c("", " ", ""))
cat("Table saved to", title, "\n")

11.9 Visualise a single gene set

plotEnrichment reproduces the classic GSEA “running enrichment score” plot for one gene set, the same style of plot produced in the previous chapter:

plotEnrichment(pathways[["HALLMARK_CHOLESTEROL_HOMEOSTASIS"]], stats) +
  labs(title = "HALLMARK_CHOLESTEROL_HOMEOSTASIS")

11.10 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.

topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n = 10), pathway] # filter table by ES > 0, order by increasing pvalue and select top 10 pathways
topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n = 10), pathway] # filter table by ES < 0, order by increasing pvalue and select top 10 pathways
topPathways <- c(topPathwaysUp, rev(topPathwaysDown)) # arrange from highest to lowest

plotGseaTable(pathways[topPathways], stats, fgseaRes, gseaParam = 0.5)

11.11 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:

collapsedPathways <- collapsePathways(fgseaRes[order(pval)][padj < 0.05], pathways, stats)
mainPathways <- fgseaRes[pathway %in% collapsedPathways$mainPathways][order(-NES), pathway]

plotGseaTable(pathways[mainPathways], stats, fgseaRes, gseaParam = 0.5)

11.12 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).