Pathway analysis: Gene Set Enrichment Analysis (GSEA)

Objectives

This notebook will demonstrate how to:

• Prepare tabular data of gene-level statistics for use with Gene Set Enrichment Analysis (GSEA), including how to remove duplicate gene identifiers
• Perform GSEA with the clusterProfiler package
• Visualize GSEA results with the enrichplot package

---

In this notebook, we will perform Gene Set Enrichment Analysis (GSEA) on the neuroblastoma cell line differential gene expression (DGE) results we generated during the RNA-seq module.

To refresh our memory, we analyzed data from Harenza et al. (2017) and specifically tested for DGE between with MYCN amplified cell lines and non-amplified cell lines using DESeq2. We have a table of results that contains our log2 fold changes and adjusted p-values.

GSEA is a functional class scoring (FCS) approach to pathway analysis that was first introduced in Subramanian et al. (2005). The rationale behind FCS approaches is that small changes in individual genes that participate in the same biological process or pathway can be significant and of biological interest. FCS methods are better suited for identifying these pathways that show coordinated changes than ORA. In ORA, we pick a cutoff that typically only captures genes with large individual changes.

There are 3 general steps in FCS methods (Khatri et al. 2012):

1. Calculate a gene-level statistic (we’ll use log2 fold change from DESeq2 here)
2. Gene-level statistics are aggregated into a pathway-level statistic
3. Assess the statistical significance of the pathway-level statistic

Other resources

• For another example using clusterProfiler for GSEA, see Intro to DGE: Functional Analysis. from Harvard Chan Bioinformatics Core Training.
• The way we’ll use clusterProfiler here uses fgsea (Fast Gene Set Enrichment Analysis) under the hood. You can read more about fgsea in Korotkevich et al. (2021).
• refine.bio examples Gene set enrichment analysis - RNA-seq from which this material has been adapted.

Set up

Libraries

# Package to run GSEA
library(clusterProfiler)

clusterProfiler v4.12.0  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/

If you use clusterProfiler in published research, please cite:
T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141

Attaching package: 'clusterProfiler'

The following object is masked from 'package:stats':

    filter

# Package that contains the MSigDB gene sets in tidy format
library(msigdbr)

Directories and Files

Directories

# Where the DGE results are stored
input_dir <- file.path("..", "RNA-seq", "results", "NB-cell")

# We will create a directory to specifically hold our GSEA results if it does
# not yet exist
output_dir <- file.path("results", "NB-cell")
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE)
}

Input files

# DGE results
dge_results_file <- file.path(input_dir,
                             "NB-cell_DESeq_amplified_v_nonamplified_results.tsv")

Output files

# GSEA pathway-level scores and statistics
gsea_results_file <- file.path(output_dir,
                               "NB-cell_gsea_results.tsv")

Gene Set Enrichment Analysis

Adapted from refine.bio examples

Figure 1. Subramanian et al. (2005).

GSEA calculates a pathway-level metric, called an enrichment score (sometimes abbreviated as ES), by ranking genes by a gene-level statistic. This score reflects whether or not a gene set or pathway is over-represented at the top or bottom of the gene rankings (Subramanian et al. 2005; Yu)

Specifically, all genes are ranked from most positive to most negative based on their statistic and a running sum is calculated: Starting with the most highly ranked genes, the running sum increases for each gene in the pathway and decreases for each gene not in the pathway. The enrichment score for a pathway is the running sum’s maximum deviation from zero. GSEA also assesses statistical significance of the scores for each pathway through permutation testing. As a result, each input pathway will have a p-value associated with it that is then corrected for multiple hypothesis testing (Subramanian et al. 2005; Yu).

The implementation of GSEA we use in here examples requires a gene list ordered by some statistic and input gene sets. When you use previously computed gene-level statistics with GSEA, it is called GSEA pre-ranked.

Gene sets

In the previous notebook, we used KEGG pathways for over-representation analysis. We identified pathways that were significantly over-represented and found that the significant pathways shared genes. FCS methods analyze each pathway independently and essentially ignore this overlap between gene sets.

MSigDB includes a collection called Hallmark gene sets (Liberzon et al. 2015) Here’s an excerpt of the collection description:

Hallmark gene sets summarize and represent specific well-defined biological states or processes and display coherent expression. These gene sets were generated by a computational methodology based on identifying gene set overlaps and retaining genes that display coordinate expression. The hallmarks reduce noise and redundancy and provide a better delineated biological space for GSEA.

We’ll use the Hallmark collection for GSEA. Notably, there are only 50 gene sets included in this collection. The fewer gene sets we test, the lower our multiple hypothesis testing burden.

We can retrieve only the Hallmark gene sets by specifying category = "H" to the msigdbr() function.

hs_hallmark_df <- msigdbr(species = "Homo sapiens",
                          category = "H")

Differential gene expression results

# Read in the DGE results
dge_results_df <- readr::read_tsv(dge_results_file)

Rows: 24912 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): gene_id, gene_symbol
dbl (5): baseMean, log2FoldChange, lfcSE, pvalue, padj

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Let’s take a peek at the contents of this file.

head(dge_results_df)

Since this data frame of DGE results includes gene symbols, we do not need to perform any kind of gene identifier conversion. We do, however, need to check for duplicate gene symbols. We can accomplish this with duplicated(), which returns a logical vector (e.g., TRUE or FALSE). The function sum() will count TRUE values as 1s and FALSE as 0s, so using it with duplicated() will count the number of duplicate values.

sum(duplicated(dge_results_df$gene_symbol))

[1] 926

This will cause a problem when we go to run GSEA.

Removing duplicates

The GSEA approach requires on discriminating between genes that are in a gene set and those that are not. Practically speaking, gene sets are just collections of gene identifiers! When the function we use for GSEA pre-ranked gets a list with duplicated gene identifiers, it can produce unexpected results.

Compared to the total number of genes that are in our results, there are not a lot of duplicates but we’ll still need to make a decision about how to handle them.

Let’s get a vector of the duplicated gene symbols so we can use it to explore our filtering steps.

duplicated_gene_symbols <- dge_results_df |>
  dplyr::filter(duplicated(gene_symbol)) |>
  dplyr::pull(gene_symbol)

Now we’ll look at the values for the the duplicated gene symbols.

dge_results_df |>
  dplyr::filter(gene_symbol %in% duplicated_gene_symbols) |>
  dplyr::arrange(gene_symbol)

We can see that the associated values vary for each row.

Let’s keep the gene symbols associated with the higher absolute value of the log2 fold change.

Retaining the instance of the gene symbols with the higher absolute value of a gene-level statistic means that we will retain the value that is likely to be more highly- or lowly-ranked or, put another way, the values less likely to be towards the middle of the ranked gene list. We should keep this decision in mind when interpreting our results. For example, if all the duplicate identifiers happened to be in a particular gene set, we may get an overly optimistic view of how perturbed that gene set is because we preferentially selected instances of the identifier that have a higher absolute value of the statistic used for ranking.

In the next chunk, we are going to filter out the duplicated rows using the dplyr::distinct() function after sorting by absolute value of the log2 fold change. This will keep the first row with the duplicated value thus keeping the row with the largest absolute value.

filtered_dge_df <- dge_results_df |>
  # Sort so that the highest absolute values of the log2 fold change are at the
  # top
  dplyr::arrange(dplyr::desc(abs(log2FoldChange))) |>
  # Filter out the duplicated rows using `dplyr::distinct()`
  dplyr::distinct(gene_symbol, .keep_all = TRUE)

Let’s see what happened to our duplicate identifiers.

# Subset to & arrange by gene symbols that were duplicated in the original
# data frame of results
filtered_dge_df |>
  dplyr::filter(gene_symbol %in% duplicated_gene_symbols) |>
  dplyr::arrange(gene_symbol)

Now we’re ready to prep our pre-ranked list for GSEA.

Pre-ranked list

The GSEA() function takes a pre-ranked (sorted) named vector of statistics, where the names in the vector are gene identifiers. This is step 1 – gene-level statistics.

lfc_vector <- filtered_dge_df |>
  # Extract a vector of `log2FoldChange` named by `gene_symbol`
  dplyr::pull(log2FoldChange, name = gene_symbol)
lfc_vector <- sort(lfc_vector, decreasing = TRUE)

Let’s look at the top ranked values.

# Look at first entries of the log2 fold change vector
head(lfc_vector)

  GAGE12D    GABBR1   FABP5P7  KIAA0355      MNX1    GAGE2A 
23.620487 23.130027 22.682631 22.562940  6.813764  6.801245 

And the bottom of the list.

# Look at the last entries of the log2 fold change vector
tail(lfc_vector)

     MUC15      TFPI2      LENG1       MT1M       DAZ2   CSPG4P4Y 
 -8.017004  -8.768669 -22.754783 -25.318887 -25.444428 -26.095612 

Run GSEA

Now for the analysis!

We can use the GSEA() function to perform GSEA with any generic set of gene sets, but there are several functions for using specific, commonly used gene sets (e.g., gseKEGG()).

gsea_results <- GSEA(geneList = lfc_vector,  # ordered ranked gene list
                     minGSSize = 25,  # minimum gene set size
                     maxGSSize = 500,  # maximum gene set set
                     pvalueCutoff = 0.05,
                     pAdjustMethod = "BH",  # correction for multiple hypothesis testing
                     TERM2GENE = dplyr::select(hs_hallmark_df,
                                               gs_name,
                                               gene_symbol))

using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

preparing geneSet collections...

GSEA analysis...

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (0.1% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

leading edge analysis...

done...

The warning about ties means that there are multiple genes that have the same log2 fold change value. This percentage is small and unlikely to impact our results. A large number of ties might tell us there’s something wrong with our DGE results (Ballereau et al. 2018).

Let’s take a look at the GSEA results.

View(gsea_results@result |>
       dplyr::arrange(dplyr::desc(NES))
)

Normalized enrichment scores (NES) are enrichment scores that are scaled to make gene sets that contain different number of genes comparable.

Let’s write these results to file.

gsea_results@result |> readr::write_tsv(gsea_results_file)

Visualizing GSEA results

We can visualize GSEA results for individual pathways or gene sets using enrichplot::gseaplot(). Let’s take a look at 3 different pathways – one with a highly positive NES, one with a highly negative NES, and one that was not a significant result – to get more insight into how ES are calculated.

Highly Positive NES

The gene set HALLMARK_MYC_TARGETS_V1 had high positive log2 fold changes. Recall a positive log2 fold change means a it had a higher expression value in MYCN amplified cell lines.

enrichplot::gseaplot(gsea_results,
                     geneSetID = "HALLMARK_MYC_TARGETS_V1",
                     title = "HALLMARK_MYC_TARGETS_V1",
                     color.line = "#0066FF")

Notice how the genes that are in the gene set, indicated by the black bars, tend to be on the left side of the graph indicating that they have positive gene-level scores.

Highly Negative NES

The gene set HALLMARK_INTERFERON_ALPHA_RESPONSE had a highly negative NES.

enrichplot::gseaplot(gsea_results,
                     geneSetID = "HALLMARK_INTERFERON_ALPHA_RESPONSE",
                     title = "HALLMARK_INTERFERON_ALPHA_RESPONSE",
                     color.line = "#0066FF")

This gene set shows the opposite pattern – genes in the pathway tend to be on the right side of the graph.

A non-significant result

The @results slot will only show gene sets that pass the pvalueCutoff threshold we supplied to GSEA(), but we can plot any gene set so long as we know its name. Let’s look at HALLMARK_P53_PATHWAY, which was not in the results we viewed earlier.

enrichplot::gseaplot(gsea_results,
                     geneSetID = "HALLMARK_P53_PATHWAY",
                     title = "HALLMARK_P53_PATHWAY",
                     color.line = "#0066FF")

Genes in the pathway are distributed more evenly throughout the ranked list, resulting in a more “middling” score.

Note: The plots returned by enrichplot::gseaplot are ggplots, so we could use ggplot2::ggsave() to save them to file if we wanted to.

Session Info

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] msigdbr_7.5.1          clusterProfiler_4.12.0 optparse_1.7.5        

loaded via a namespace (and not attached):
  [1] DBI_1.2.2               gson_0.1.0              shadowtext_0.1.3       
  [4] gridExtra_2.3           rlang_1.1.3             magrittr_2.0.3         
  [7] DOSE_3.30.0             compiler_4.4.1          RSQLite_2.3.6          
 [10] png_0.1-8               vctrs_0.6.5             reshape2_1.4.4         
 [13] stringr_1.5.1           pkgconfig_2.0.3         crayon_1.5.2           
 [16] fastmap_1.1.1           XVector_0.44.0          labeling_0.4.3         
 [19] ggraph_2.2.1            utf8_1.2.4              HDO.db_0.99.1          
 [22] rmarkdown_2.26          tzdb_0.4.0              enrichplot_1.24.0      
 [25] UCSC.utils_1.0.0        purrr_1.0.2             bit_4.0.5              
 [28] xfun_0.43               zlibbioc_1.50.0         cachem_1.0.8           
 [31] aplot_0.2.2             GenomeInfoDb_1.40.0     jsonlite_1.8.8         
 [34] blob_1.2.4              highr_0.10              BiocParallel_1.38.0    
 [37] tweenr_2.0.3            parallel_4.4.1          R6_2.5.1               
 [40] bslib_0.7.0             stringi_1.8.3           RColorBrewer_1.1-3     
 [43] jquerylib_0.1.4         GOSemSim_2.30.0         Rcpp_1.0.12            
 [46] knitr_1.46              readr_2.1.5             IRanges_2.38.0         
 [49] Matrix_1.7-0            splines_4.4.1           igraph_2.0.3           
 [52] tidyselect_1.2.1        qvalue_2.36.0           yaml_2.3.8             
 [55] viridis_0.6.5           codetools_0.2-20        lattice_0.22-6         
 [58] tibble_3.2.1            plyr_1.8.9              treeio_1.28.0          
 [61] Biobase_2.64.0          withr_3.0.0             KEGGREST_1.44.0        
 [64] evaluate_0.23           gridGraphics_0.5-1      scatterpie_0.2.2       
 [67] getopt_1.20.4           polyclip_1.10-6         Biostrings_2.72.0      
 [70] ggtree_3.12.0           pillar_1.9.0            stats4_4.4.1           
 [73] ggfun_0.1.4             generics_0.1.3          vroom_1.6.5            
 [76] hms_1.1.3               S4Vectors_0.42.0        ggplot2_3.5.1          
 [79] tidytree_0.4.6          munsell_0.5.1           scales_1.3.0           
 [82] glue_1.7.0              lazyeval_0.2.2          tools_4.4.1            
 [85] data.table_1.15.4       fgsea_1.30.0            babelgene_22.9         
 [88] fs_1.6.4                graphlayouts_1.1.1      fastmatch_1.1-4        
 [91] tidygraph_1.3.1         cowplot_1.1.3           grid_4.4.1             
 [94] ape_5.8                 tidyr_1.3.1             AnnotationDbi_1.66.0   
 [97] colorspace_2.1-0        nlme_3.1-164            patchwork_1.2.0        
[100] GenomeInfoDbData_1.2.12
 [ reached getOption("max.print") -- omitted 20 entries ]