Gene set enrichment and pathway analysis

Motivation¶

Single-cell RNA-seq (scRNA-seq) provides unprecedented insights into molecular variation, but it is hard to interpret the resulting high-dimensional data, since hundreds of genes may explain the observed variability in a particular biological context.

Gene set enrichment analysis (GSEA) allows to summarize the high-dimensional molecular measurements into interpretable terms such as biological processes or pathways, which enable a better understanding of the obtained results and the generation of mechanistic hypothesis. Note that GSEA is also known as pathway analysis, enrichment analysis, gene set analysis and functional analysis. Accordingly, we use these terms interchangeably throughout this chapter.

Functional analysis requires three different inputs:

• Molecular readouts: a matrix containing gene expression or gene-level statistics (log fold changes, test statistic).

• Prior knowledge: collection of gene sets or footprints associated to biological processes (e.g. epithelial to mesechymal transition, metabolism etc.) or regulatory pathways (e.g. MAPK signalling).

• Inference method: statistical method that computes the “enrichment” scores of the gene sets from the molecular data.

Background¶

Gene set collections¶

Gene sets are a curated list of genes that are known to be involved in a biological process through previous studies and/or experiments and are stored in prior knowledge databases. The Molecular Signatures Database (MSigDB) Subramanian et al., 2005Liberzon et al., 2011 is the most comprehensive database consisting of 9 collections of gene sets. Some commonly used collections are C5, which is the gene ontology (GO) collection, C2 collection of curated gene signatures from published studies that are typically context (e.g. tissue, condition) specific, but also include KEGG and REACTOME gene signatures. The Hallmark collection were developed to reduce the redundancy between gene sets, and for immunologic studies the C7 collection is a common choice. Note that these signatures are mainly derived from Bulk-seq measurements and measure continuous phenotypes. Recently and with the wide-spread availability of scRNA-seq datasets, databases have evolved that provide curated marker lists derived from published single cell studies that define cell types in various tissues and species. These include CellMarker Zhang et al., 2019 and PanglaoDB Franzén et al., 2019. Curated marker lists are not limited to those made available in databases, and can be curated by oneself using own or public data.

Gene sets contain experimentally validated protein members of a given pathway. Running enrichment on those assumes a strong correlation between gene expression, protein abundance, and protein activity. Despite this, most gene sets still provide reliable scores because their members are usually regulated together Szalai & Saez-Rodriguez, 2020. On the other hand, footprints are sets that instead contain genes that changed downstream of said process when it was active, bringing them closer to transcriptomics data. Footprints usually contain interaction weights for each gene in the set, denoting its strength and mode of regulation (positive: activation, negative: inhibition). This concept has been applied to transcription factors, kinases and enzymes Dugourd & Saez-Rodriguez, 2019. Known databases of footprints include PROGENy Schubert et al., 2018 for signaling pathways, CytoSig Jiang et al., 2021 for ligand signatures, and DoRothEA Garcia-Alonso et al., 2019 for transcription factors. Similar to gene sets, these can be curated by oneself derived in a data-driven manner, for example by generating Gene Regulatory Networks.

Applications in single-cell data¶

In scRNA-seq data analysis, enrichment is generally carried out on clusters of cells or cell types, one-at-a-time. For example, one can test how genes change between untreated and treated cells per cell-type and use the obtained gene statistics to see what biological processes are deregulated Squair et al., 2021Crowell et al., 2020. Enrichment can be assessed in any complex experimental design while accounting for batch effects, between-individual variations, gender, strain differences in mouse models, etc. (check the differential gene expression analysis chapter for best practices).

Alternatively, enrichment methods can also be applied directly to the gene expression data of individual cells with the main objective to reduce highly-dimensional expression data into a collection of interpretable latent variables to explore their variability. This dimensionality reduction driven by prior knowledge captures transcriptional heterogeneity and variability at the single-cell level, increasing the power of population comparisons, while keeping biological significance. Compared to other data-driven dimensionality reduction methods like principal component analysis (PCA) or factor analysis, the estimated latent variables generated from enrichment methods are not necessarily independent.

Most enrichment methods have been previously developed for bulk RNA-seq but can be applied to single-cell data. Holland et al. (2020) found that methods have optimal performance in simulated scRNA-seq data, and even partially outperform tools specifically designed for scRNA-seq analysis despite the drop-out events and low library sizes in single-cell data. Holland et al. also concluded that enrichment inference is more sensitive to the choice of gene sets (the prior knowledge used) rather than the statistical methods. This observation though can be explained by the fact that biological processes might be context-specific; for example TF-target associations in one cell-type may actually differ from another cell-type or tissue. There are different strategies to contextualize gene sets, such as inferring gene regulatory networks, or performing factor analysis influenced by prior, as seen in Muvi Qoku & Buettner, 2022 and Spectra Kunes et al., 2022.

In contrast to Holland et al., Zhang et al. Zhang et al. (2020) found that single-cell-based tools, specifically Pagoda2, outperform some bulk-base methods from three different aspects of accuracy, stability and scalability. It should be noted that functional inference tools inherently do not account for batch effects or biological variations other than the biological variation of interest. Therefore, it is up to the data analyst to ensure that the differential gene expression analysis step has worked properly.

Enrichment methods¶

Most enrichment methods are based on hypergeometric tests, rank statistics or empirical estimations of likelihoods based on permutations. Regardless of the statistical method used, their hypothesis tests can be divided in two main classes: competitive or self-contained as defined by Goeman & Bühlmann (2007). A competitive test tests whether the genes in the set are highly ranked relative to the genes not in the set (usually of the same size). The sampling unit here is genes, so the test can be done with a single sample or contrast vector. In self-contained gene set testing, the sampling unit is the subject, so multiple samples per group are required, but it is not required to have genes that are not present in the set. A self-contained test tests whether genes in the test set are differentially expressed without regard to any other gene measured in the dataset. Competitive hypotheses are the most used in enrichment analysis because they depend solely on a single vector of gene-level statistics compared to self-contained methods which are limited by the number of samples Mathur et al., 2018. In a recent benchmark Geistlinger et al., 2020, it has been shown that self-contained methods might be quite sensitive since they can identify gene sets as enriched while containing only a single differentially expressed gene. On the other hand, competitive methods might be more restrictive by testing for excess of differential expression in the gene set when compared to the background level, coming closer to the intuition of an enrichment and tend to rank relevant gene sets systematically higher than self-contained methods. However, it is important to keep in mind that each approach has its pros an cons and that these distinctions between the two null hypotheses mainly affect the interpretation of gene set enrichment results. In practice, competitive methods dominate for their stability, since removing even a single sample can swing self-contained results. For this reason, we will only work with competitive methods in this chapter.

Common competitive enrichment methods include simple hypergeomtric tests or one-tailed Fisher’s exact test (as in Enrichr Chen et al., 2013), the GSEA algorithm Subramanian et al., 2005 and its fast implementation (fgsea) Korotkevich et al., 2021, GSVA Hänzelmann et al., 2013, VISION DeTomaso et al., 2019, AUCell Aibar et al., 2017, Pagoda2 Fan et al., 2016Lake et al., 2018, camera (as in limma Ritchie et al., 2015), scDECAF and simple combined z-score Lee et al., 2008. As mentioned before, while less computational tools that implement these methods use self-contained hypothesis tests, notable examples are roast and fry, both available in the differential testing framework limma Ritchie et al., 2015. These methods leverage linear models and Empirical Bayes moderation of test statistics Smyth, 2005. Note that while most methods work downstream of differential gene testing, methods included in the limma framework (camera, roast and fry) start from the sample level gene expression data. Detailed explanations on these methods can be found in the limma user manual.

There are also enrichment methods that leverage footprints when inferring scores, i.e. they accommodate the usage of gene weights. Common competitive enrichment methods that leverage footprints include VIPER Alvarez et al., 2016, ULM and MLM Badia-i-Mompel et al., 2022. For self-contained enrichment methods, fry and roast do it too.

Enrichment frameworks collect several enrichment methods in a single tool. These facilitate the usage of various enrichment methods on the data since they have a common format and provide utility functions to access prior knowledge resources. Some of them are Piano Väremo et al., 2013 and EGSEA Alhamdoosh et al., 2016 (available in R), and decoupler Badia-i-Mompel et al., 2022 (available both in R and Python).

Technical considerations¶

Filtering out the gene sets with low number of genes¶

A common practice is to exclude any gene sets with a few genes overlapping the data or Highly Variable Genes (HVG) in the pre-processing step. Zhang et al. (2020) found that the performance of both single-cell-based and bulk-based methods drops as gene coverage, that is the number of genes in pathways/gene sets, decreases. Holland et al. (2020) also found that gene sets of smaller size adversely impacts the performance of Bulk-seq enrichment methods on single cell data. These report collectively support that filtering gene sets with few gene members, say less than 10 or 15 genes in the set, is beneficial in functional analysis. Damian & Gorfine (2004) attributed this to the fact that gene variances in gene sets with a smaller number of genes are more likely to be large, whereas gene variances in larger gene sets tend to be smaller. This impacts the accuracy of the test statistics computed to test for enrichment. Zhang et al. additionally found that pathway analysis was susceptible to normalization procedures applied to gene expression measurements.

Data normalization¶

Read counts in single cell experiments are typically normalised early on in the pre-processing pipeline to ensure that measuerments are comparable across cells of various library sizes. Zhang et al. (2020) found that normalisation by SCTransform Hafemeister & Satija, 2019 and scran Lun et al., 2016 generally improves the performance of both single-cell- and bulk-based enrichment scoring tools. They found that the performance of AUCell (a rank-base method) and z-score (transformation to zero mean, unit standard deviation) is particularly affected by normalization with distinct methods.

Functional analysis between conditions¶

Prepare and explore the data¶

We first download the 25K PBMC data and follow the standard scanpy workflow for normalization of read counts and subsetting on the highly variable genes. The dataset contains untreated and IFN-$\beta$ stimulated human PBMC cells Kang et al., 2018. We explore patterns of variation in the data with UMAP representation of 4000 highly variable genes.

While the current object comes with UMAP and PCA embeddings, these have been corrected for stimulation condition, which we don’t want for this analysis. Instead we will recompute these.

Retrieving gene sets¶

Gene sets can be downloaded from their respective databases, usually they are stored in the GMT file format. Gene sets can also be retrieved from databases using wrappers, here is an example on how to download MSigDB using decoupler (which leverages the meta-database OmniPath Türei et al., 2021):

# Retrieving via python
msigdb = dc.op.resource("MSigDB", organism='human')

# Get reactome pathways
reactome = msigdb.query("collection == 'reactome_pathways'")
# Filter duplicates
reactome = reactome[~reactome.duplicated(("geneset", "genesymbol"))]

However, for stability of this tutorial we are using a fixed version of the gene set collection. Here we download and read the GMT file for the REACTOME pathways annotated in the C2 collection of MSigDB:

Here we remove potential noisy gene sets by removing sets with too low (<15) or too high (>500) number of members. Note that these values are completely arbitrary and should be adjusted to the question at hand. As a rule of thumb, we recommend exploring the distribution of gene set lengths to get an idea if there are any outliers.

Cell-level gene set scoring¶

As previously described, gene set scores can be inferred in each individual cell, that is based on absolute gene expression in the cell, regardless of expression of genes in the other cells.

Inference using gene sets with AUCell¶

To demonstrate this concept, we will start by using the method AUCell Aibar et al., 2017 through decoupler.

We now add the scores for the interferon-related REACTOME pathways to the obs field of the AnnData object and annotate the score level of these pathways in each of the cells on the UMAP:

We can also visualize this as distributions in violinplots:

AUCell scores the pathways well-known to be implicated in interferon signalling high in IFN-stimulated cells (except for Megakaryocytes which is expected), while cells in the control condition generally have lower scores for these pathways, demonstrating that gene set scoring with AUCell has been successful.

Inference using footprints with ULM¶

Additionally, we can infer scores using footprints instead of gene sets. Here we will leverage the database PROGENy Schubert et al., 2018. Here is how to access it through decoupler:

# Retrieving via python
progeny = dc.op.progeny(organism='human', top=100)

However, again for stability of this tutorial we are using a fixed version.

Note that in this resource, each gene has a weight associated to the pathway it belongs, denoting its mode of regulation (positive or negative) and its importance.

Now lets estimate pathway scores using ULM from decoupler, a method that accounts for weighted gene sets:

Next we add the scores for the immune-related PROGENy pathways to the obs field of the AnnData object and annotate the score level of these pathways in each of the cells on the UMAP:

As before, we can also visualize this in violinplots:

In this case, the pathway JAK-STAT, triggered by the binding of interferons such as INF-$\beta$, is the one with the highest score after stimulation. The other immune pathways, NFkB and TNFa, remain overall the same. The concordance between pathway scores by AUCell and ULM is promising, given that we know a priori that IFN-related pathways should be the top-ranked terms.

Population contrast enrichment scoring¶

Previously, we enriched functional terms for individual cells. Now, we instead compare cell populations. For this, we pseudobulk our samples and cell type combinations and perform DGE on these pseudobulks to mitigate inflated p-values Squair et al., 2021. We use pydeseq2 for the DGE analysis.

This analysis therefore builds directly on a DGE workflow. If you have not yet performed DGE on your data, see the DGE analysis chapter for details. In this tutorial, we directly load the previously computed results from our DGE analysis of CD14+ monocytes between conditions.

While any statistic can be used, we recommend using t-values rather than log2FCs, as t-values account for the significance of the change. We will transform the obtained t-values, stored in the column "stat", into a wide matrix format so that it can be used with decoupler.

Inference using gene sets with GSEA¶

Once we prepared our gene statistics, we will use those to infer pathway scores using GSEA Subramanian et al., 2005 from decoupler. Running dc.mt.gsea will return enrichment scores and adjusted p-values using Benjamini-Hochberg.

We can then visualize the most active and inactive pathways.

Pathways are shown on the y-axis, and the x-axis shows the GSEA enrichment score. Enrichment scores are signed: positive scores indicate pathway upregulation, while negative scores indicate downregulation. Larger absolute values indicate stronger enrichment. Only the top 15 pathways are displayed, ranked by direction and magnitude. The majority of interferon-related pathways are indeed positively ranked and in the right direction.

Inference using footprints with ULM¶

Let us leverage DGE result again, but this time with the PROGENy footprints and the ULM method.

Like in GSEA, using ULM returns enrichment scores and adjusted p-values using Benjamini-Hochberg. Scores are signed and show the magnitute of up- or downregulation. While only JAK-STAT and NFkB show significant results, we want to plot all 13 pathways available in PROGENy:

As seen in the previous section, JAK-STAT, a INF-$\beta$ induced pathway, is the one with the highest score in stimulated CD14+ monocytes compared to the control ones.

Pseudotime Enrichment Scoring¶

Analyzing omics data is not always limited to comparing two discrete groups, such as treatment versus control. Many biological processes, including cell differentiation, disease progression, and development, are continuous. In these cases, it is important to evaluate how features and their enrichment scores change along the continuous trajectory.

We apply this approach to the adult human bone marrow dataset Setty et al., 2019 from a previous chapter, where we concluded that the Palantir pseudotime should be used for subsequent analyses.

Here, instead of using gene sets representing specific pathways, we score transcription factors (TFs) based on gene regulatory networks (GRNs). TFs are genes whose protein products bind DNA and regulate the expression of target genes, either promoting or inhibiting transcription. Because TF transcript levels often do not reflect their actual activity, we assess TF activity through enrichment of their target genes Badia-i-Mompel et al., 2023. GRNs capture these TF-gene interactions and therefore provide a more accurate representation of regulatory activity. Again, the scores are signed and reflect the magnitude of the effect, indicating the extent to which TFs are activated or deactivated. For more details, please see our GRN chapter.

To do so, we use the CollecTRI network, a comprehensive resource containing a curated collection of TFs and their target genes Müller-Dott et al., 2023. The network can be accessed through decoupler as follows:

collectri = dc.op.collectri(organism="human")

As before, to ensure the stability and reproducibility of this tutorial, we use a fixed version of the network.

TF scores can be easily computed by running the ulm method.

Scores can then be extracted as a new AnnData object.

Next, TFs associated with the inferred pseudotime can be identified.

The top 5 TF markers associated to the pseudotime can then be extracted.

And visualized.

For example, LEF-1 plays a pivotal role in lymphoid differentiation Petropoulos et al., 2008, making its activation in common lymphoid progenitors (CLPs) biologically plausible. Among other factors, FOSB, a subunit of AP-1, plays a role in the functional development of hematopoietic precursor cells into mature blood cells, supporting increased activation at later pseudotime points Liebermann et al., 1998.

Prior to generating more complex visualizations, cells are grouped into bins based on their ordering to facilitate plotting.

Which can be plotted as lines.

Or as a matrix.

Additionally, to better understand the obtained enrichment scores, the targets of a specified TF can be plotted along the trajectory.

The change in the enrichment score results from positive targets of FOSB increasing in expression, while its repressed targets decrease in expression along the trajectory. Consequently, the TF is inactive at the beginning but becomes active by the end.

Spatial Enrichment Analysis¶

Finally, we can also perform GSEA on spatial transcriptomics data. The dataset we use comprises approximately 4,000 mini-bulks (spots) for ~15k genes from a chronic active multiple sclerosis lesion, profiled using Visium Lerma-Martin et al., 2024. Here, we again aim to characterize TF activity using CollecTRI, but this time in its spatial context.

The following step normalizes the data and stores it into a layer for later use.

The spot metadata stored in obs can be inspected.

The dataset also includes tissue annotations from the original publication:

• GM: Grey Matter

• LC: Lesion Core

• LR: Lesion Rim

• PPWM: PeriPlaque White Matter

• VI: Vascular Infiltrating

The H&E image and the different annotated regions (niches) can be visualized.

Spatial Weighting¶

Spatial data are often sparse. To address this, spatial information can be leveraged to impute gene expression based on the expression levels of nearby spatial neighbors Dimitrov et al., 2024. This is achieved by generating spatial weights for each spot and applying them to the gene expression values, effectively smoothing the data.

The spatial weights assigned to each spot can then be visualized.

It can be observed that for a given spot, the spatial weight is maximal at the spot itself and decreases with increasing distance to other spots. Adjust the bandwith (bw) to increase or decrease the radius to consider. The gene expression values can now be spatially weighted.

The original log-normalized counts can be plotted alongside the newly spatially transformed values.

After smoothing, spots surrounded by other spots expressing the same gene will have higher values compared to spots lacking neighboring spots with high values, even if those spots themselves had high values of expression.

Gene expression for these genes appears to exhibit spatial compartmentalization.

Scoring¶

We can now again use the ulm method to compute TF scores using collectri, which we loaded previously.

And extract the scores.

The obtained scores can be visualized.

Here, the TF RFX5, a key regulator for antigen-presenting cells, exhibits a strong spatial pattern, showing greater activity in the lesion rim (LR) compared to other areas of the tissue. This aligns with the known fact that these cells are present in this region to engulf cellular debris coming from the lesion core.

This can be compared with the actual expression values of the RFX5 gene.

When we look at the expression of the TF RFX5 there is no spatial pattern, as it is randomly expressed across the tissue. TF expression tends to be sporadic and insufficiently informative to determine whether a TF is active in a cell, underscoring the importance of using enrichment scores for interpreting biological data Badia-i-Mompel et al., 2023.

Next, marker TFs for each niche can be identified.

The top 3 TF markers per niche can then be extracted.

The obtained markers can be plotted

Individual TFs can be examined by plotting their score distributions.

Here, the TF activities for OLIG1 can be observed. OLIG1 is a known marker TF for oligodendrocytes, which are present in the PPWM.

A note on the redundancy between gene sets¶

Generally, there can be a large overlap between closely related gene sets. The overlap could be due to gene sets belonging to the same ontology, or describing biological processes that are involved in similar functions. For example, in this dataset we have seen the gene sets REACTOME_INTERFERON_ALPHA_BETA_SIGNALING and REACTOME_INTERFERON_SIGNALING both being enriched, the former being the specific immune response and the later is the general one. This overlap impacts the rank of the gene sets in the enrichment results. It is important to interpret what is the general state of the biological system being studied, in this case it is clear that immune signaling might be happening. Moreover, it might be beneficial to select tailored gene sets and footprint resources, to reduce the number of possible false positives. Enrichment frameworks improve biological interpretation by providing testing procedures with different null hypotheses, and including prior knowledge in the computational methods.

Questions¶

Flipcards¶

Multiple choice questions¶

Contributors¶

We gratefully acknowledge the contributions of:

Authors¶

• Soroor Hediyeh-Zadeh

• Pau Badia-i-Mompel

• Isaac Virshup

• Luis Heinzlmeier

Reviewers¶

• Lukas Heumos

• Anastasia Litinetskaya