Option 1: Built-In gene sets

data("escape.gene.sets", package="escape")

Option 2: MSigDB via getGeneSets()

Gene set enrichment analysis begins by identifying the appropriate gene sets for your study. The getGeneSets() function simplifies this process by extracting one or more gene set libraries from the Molecular Signature Database (MSigDB) and returning them as a GSEABase GeneSetCollection object. Note that the first time you run getGeneSets(), it downloads a complete local copy of the gene sets, which may take a little while. Future calls will use the cached version, greatly improving performance.

To retrieve gene set collections from MSigDB, specify the library or libraries of interest using the library parameter. For example, to select multiple libraries, use:

In addition, the function supports further subsetting through these parameters:

• subcategory: Narrow down your selection by specifying subcategories within a library. Examples include “CGN”, “CGP”, “CP:BIOCARTA”, “CP:KEGG”, “GO:BP”, “IMMUNESIGDB”, etc.
• gene.sets: Isolate individual pathways or gene sets by providing their specific names.

GS.hallmark <- getGeneSets(library = "H")

Option 3: Define personal gene sets

gene.sets <- list(Bcells = c("MS4A1","CD79B","CD79A","IGH1","IGH2"),
                        Myeloid = c("SPI1","FCER1G","CSF1R"),
                        Tcells = c("CD3E", "CD3D", "CD3G", "CD7","CD8A"))

Option 4: Using msigdbr

msigdbr is an alternative R package to access the Molecular Signature Database in R. There is expanded support for species in the package as well as a mix of accessible versus downloadable gene sets, so it can be faster than caching a copy locally.

GS.hallmark <- msigdbr(species  = "Homo sapiens",  
                       category = "H")

ssGSEA

This method calculates the enrichment score using a rank-normalized approach and generating an empirical cumulative distribution function for each individual cell. The enrichment score is defined for a gene set (G) using the number of genes in the gene set (NG) and total number of genes (N).

\[ ES(G,S) \sum_{i = 1}^{n} [P_W^G(G,S,i)-P_{NG}(G,S,i)] \]

Please see the following citation for more information.

GSVA

GSVA varies slightly by estimating a Poisson-based kernel cumulative density function. But like ssGSEA, the ultimate enrichment score reported is based on the maximum value of the random walk statistic. GSVA appears to have better overall consistency and runs faster than ssGSEA.

\[ ES_{jk}^{max} = V_{jk} [max_{l=1,...,p}|v_{jk}(l)] \]

Please see the following citation for more information.

AUCell

In contrast to ssGSEA and GSVA, AUCell takes the gene rankings for each cell and step-wise plots the position of each gene in the gene set along the y-axis. The output score is the area under the curve for this approach.

Please see the following citation for more information.

UCell

UCell calculates a Mann-Whitney U statistic based on the gene rank list. Importantly, UCell has a cut-off for ranked genes (\[r_{max}\]) at 1500 - this is per design as drop-out in single-cell can alter enrichment results. This also substantially speeds the calculations up.

The enrichment score output is then calculated using the complement of the U statistic scaled by the gene set size and cut-off.

\[ U_j^` = 1-\frac{U_j}{n \bullet r_{max}} \]

Please see the following citation for more information.

Choosing a backend

Each method above runs on escape’s own engine by default. Setting backend = "plaid" reroutes the calculation through plaid, which reimplements these scores as a single sparse matrix operation and uses substantially less memory. On a 166,000-spot Xenium SpatialExperiment with three gene sets, UCell scoring took 28.6 s natively and 6.8 s through plaid, a 4.2x speedup. plaid also provides three scores escape has no native implementation for, which select the plaid backend on their own.

Those correlations come from a simulated 2000-gene, 120-cell matrix with two planted gene sets. On narrower matrices they drop: pbmc_small (230 genes) gives roughly 0.79 for UCell and 0.50 for ssGSEA.

# native (default)
scores <- escape.matrix(pbmc_small, gene.sets = GS.hallmark, method = "UCell")

# same method, plaid engine, log-normalized input picked automatically
fast <- escape.matrix(pbmc_small, gene.sets = GS.hallmark, method = "UCell",
                      backend = "plaid")

# method-specific tuning goes through backend.args, never through ...
tuned <- escape.matrix(pbmc_small, gene.sets = GS.hallmark, method = "ssGSEA",
                       backend = "plaid", backend.args = list(alpha = 0.25))

Two arguments are easy to confuse. normalize is escape’s post-hoc drop-out scaling (see performNormalization() below). backend.args$normalize is plaid’s median normalization of the raw scores, applied inside the calculation. They are independent, and turning both on stacks them.

Install plaid with BiocManager::install("plaid"). If it is missing, the native path is unaffected and any backend = "plaid" call errors with install instructions. Anyone using a plaid backend should cite Zito et al., Bioinformatics 2025, 41(12):btaf621 in addition to the original method paper and escape.

Choosing the input assay

input.assay controls which expression matrix is scored. The default "auto" reads raw counts for the native backend and log-normalized values for plaid, which is what each expects. "counts" and "logcounts" resolve to the right place for both object types - the counts and data layers in Seurat, the counts and logcounts assays in a SummarizedExperiment - so the same call works across classes. Any other string is taken literally as a layer or assay name.

Note that the four native methods are rank or CDF based, so a monotone transform of the input generally leaves their scores unchanged. input.assay matters most for the plaid paths, which average in log space.

Spatial data

SpatialExperiment objects work anywhere a SingleCellExperiment does, spots being treated as cells. Nothing special is required - pass the object straight to runEscape().

To keep this vignette self-contained, we build a small SpatialExperiment from pbmc_small with made-up coordinates rather than downloading a real slide:

library(SpatialExperiment)

set.seed(42)
spe <- SpatialExperiment(
  assays        = list(counts = counts(sce.pbmc)),
  colData       = colData(sce.pbmc),
  spatialCoords = matrix(runif(2 * ncol(sce.pbmc)), ncol = 2,
                         dimnames = list(colnames(sce.pbmc), c("x", "y")))
)

spe <- runEscape(spe,
                 gene.sets      = escape.gene.sets,
                 method         = "UCell",
                 min.size       = 3,
                 new.assay.name = "escape")

spe

## class: SpatialExperiment 
## dim: 230 80 
## metadata(1): escape_backend
## assays(1): counts
## rownames(230): MS4A1 CD79B ... SPON2 S100B
## rowData names(0):
## colnames(80): ATGCCAGAACGACT CATGGCCTGTGCAT ... GGAACACTTCAGAC
##   CTTGATTGATCTTC
## colData names(9): orig.ident nCount_RNA ... ident sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(1): escape
## spatialCoords names(2) : x y
## imgData names(0):

runEscape() attaches the scores as an altExp named by new.assay.name and leaves the rest of the object alone - spatialCoords(), imgData(), and sample_id all survive unchanged:

altExpNames(spe)

## [1] "escape"

dim(spatialCoords(spe))

## [1] 80  2

The scores themselves carry no spatial information. Coordinates stay on the parent object, so to plot enrichment against position you pull both from spe:

scores <- t(assay(altExp(spe, "escape")))

df <- data.frame(spatialCoords(spe),
                 Proinflammatory = scores[, "Proinflammatory"])

ggplot(df, aes(x = x, y = y, color = Proinflammatory)) +
  geom_point(size = 2) +
  scale_color_gradientn(colors = hcl.colors(7, "inferno")) +
  coord_fixed() +
  theme_classic()

Objects carrying several sample_id values need no special treatment. Note that all downstream escape functions - performNormalization(), performPCA(), and the visualization gallery below - read the altExp the same way they would for a SingleCellExperiment, so the rest of this vignette applies unchanged.

escape.matrix

escape has 2 major functions - the first being escape.matrix(), which serves as the backbone of enrichment calculations. Using count-level data supplied from a single-cell object or matrix, escape.matrix() will produce an enrichment score for the individual cells with the gene sets selected and output the values as a matrix.

method

• AUCell
• GSVA
• ssGSEA
• UCell
• PLAID, singscore, scSE (these require the plaid backend and select it on their own)

groups

• The number of cells to calculate at once.

min.size

• The minimum size of detectable genes in a gene set. Gene sets less than the min.size will be removed before the calculation.

normalize

• Use the number of genes from the gene sets in each cell to normalize the enrichment scores. The default value is FALSE.

make.positive

• During normalization, whether to shift the enrichment values to a positive range (TRUE) or not (FALSE). The default value is FALSE.

Cautionary note: make.positive was added to allow for differential analysis downstream of enrichment as some methods may produce negative values. It preserves log-fold change, but ultimately modifies the enrichment values and should be used with caution.

backend

• Which engine computes the scores, “native” (the default) or “plaid”. See Choosing a backend above.

backend.args

• A named list of method-specific tuning passed to the underlying plaid function, e.g. list(alpha = 0.25). Names are checked against the target function, so a typo errors rather than being ignored. Ordinary ... arguments go to the native engines and are refused on the plaid path.

input.assay

• Which expression matrix to score. “auto” (the default) reads counts for the native backend and logcounts for plaid. See Choosing the input assay.

enrichment.scores <- escape.matrix(pbmc_small, 
                                   gene.sets = escape.gene.sets, 
                                   groups = 1000, 
                                   min.size = 5)

ggplot(data = as.data.frame(enrichment.scores), 
      mapping = aes(enrichment.scores[,1], enrichment.scores[,2])) + 
  geom_point() + 
  theme_classic() + 
  xlab(colnames(enrichment.scores)[1]) + 
  ylab(colnames(enrichment.scores)[2]) 

Multi-core support is for all methods is available through BiocParallel. To add more cores, use the argument BPPARAM to escape.matrix(). Here we will use the SnowParam() for it’s support across platforms and explicitly call 2 workers (or cores).

enrichment.scores <- escape.matrix(pbmc_small, 
                                   gene.sets = escape.gene.sets, 
                                   groups = 1000, 
                                   min.size = 3, 
                                   BPPARAM = BiocParallel::SnowParam(workers = 2))

runEscape

Alternatively, we can use runEscape() to calculate the enrichment score and directly attach the output to a single-cell object. The additional parameter for ``runEscape is new.assay.name, in order to save the enrichment scores as a custom assay in the single-cell object.

pbmc_small <- runEscape(pbmc_small, 
                        method = "ssGSEA",
                        gene.sets = escape.gene.sets, 
                        groups = 1000, 
                        min.size = 3,
                        new.assay.name = "escape.ssGSEA")

sce.pbmc <- runEscape(sce.pbmc, 
                      method = "UCell",
                      gene.sets = escape.gene.sets, 
                      groups = 1000, 
                      min.size = 5,
                      new.assay.name = "escape.UCell")

We can quickly examine the attached enrichment scores using the visualization/workflow we prefer - here we will use just FeaturePlot() from the Seurat R package.

#Define color palette 
colorblind_vector <- hcl.colors(n=7, palette = "inferno", fixup = TRUE)

FeaturePlot(pbmc_small, "Proinflammatory") + 
  scale_color_gradientn(colors = colorblind_vector) + 
  theme(plot.title = element_blank())

performNormalization

Although we glossed over the normalization that can be used in escape.matrix() and runEscape(), it is worth mentioning here as normalization can affect all downstream analyses.

There can be inherent bias in enrichment values due to drop out in single-cell expression data. Cells with larger numbers of features and counts will likely have higher enrichment values. performNormalization() will normalize the enrichment values by calculating the number of genes expressed in each gene set and cell. This is similar to the normalization in classic GSEA and it will be stored in a new assay.

pbmc_small <- performNormalization(input.data = pbmc_small, 
                                   assay = "escape.ssGSEA", 
                                   gene.sets = escape.gene.sets)

An alternative for scaling by expressed gene sets would be to use a scaling factor previously calculated during normal single-cell data processing and quality control. This can be done using the scale.factor argument and providing a vector.

pbmc_small <- performNormalization(input.data = pbmc_small, 
                                   assay = "escape.ssGSEA", 
                                   gene.sets = escape.gene.sets, 
                                   scale.factor = pbmc_small$nFeature_RNA)

performNormalization() has an additional parameter make.positive. Across the individual gene sets, if negative normalized enrichment scores are seen, the minimum value is added to all values. For example if the normalized enrichment scores (after the above accounting for drop out) ranges from -50 to 50, make.positive will adjust the range to 0 to 100 (by adding 50). This allows for compatible log2-fold change downstream, but can alter the enrichment score interpretation.

Three details are worth knowing:

• It always reads raw counts. The scale factor is the number of genes in each set with a non-zero count, so performNormalization() ignores whatever input.assay was used to compute the scores. Detection is identical in count and log space, so this is deliberate.
• Supplied scores win. If you pass enrichment.data explicitly it is used even when the named assay also exists on the object, and you get a warning that both were available.
• Setting assay = NULL returns a matrix rather than attaching a new assay to the object. This is how escape.matrix(normalize = TRUE) calls it internally.

Gene-set names are matched to the enrichment columns literally first, then against the hyphenated form Seurat produces when it builds an assay. That means underscored library names such as HALLMARK_HYPOXIA normalize correctly whether the scores live in a Seurat assay, an altExp, or a bare matrix.

heatmapEnrichment

We can examine the enrichment values across our gene sets by using heatmapEnrichment(). This visualization will return the mean of the group.by variable. As a default - all visualizations of single-cell objects will use the cluster assignment or active identity as a default for visualizations.

heatmapEnrichment(pbmc_small, 
                  group.by = "ident",
                  gene.set.use = "all",
                  assay = "escape.ssGSEA")

Most of the visualizations in escape have a defined set of parameters.

group.by

• The grouping variable for the comparison.

facet.by

• Using a variable to facet the graph into separate visualizations.

scale

• TRUE - z-transform the enrichment values.
• FALSE - leave raw values (DEFAULT).

In addition, heatmapEnrichment() allows for the reclustering of rows and columns using Euclidean distance of the enrichment scores and the Ward2 methods for clustering using cluster.rows and cluster.columns.

heatmapEnrichment(sce.pbmc, 
                  group.by = "ident",
                  assay = "escape.UCell",
                  scale = TRUE,
                  cluster.rows = TRUE,
                  cluster.columns = TRUE)

Each visualization has an additional argument called **palette that supplies the coloring scheme to be used - available color palettes can be viewed with hcl.pals().

hcl.pals()

##   [1] "Pastel 1"      "Dark 2"        "Dark 3"        "Set 2"        
##   [5] "Set 3"         "Warm"          "Cold"          "Harmonic"     
##   [9] "Dynamic"       "Grays"         "Light Grays"   "Blues 2"      
##  [13] "Blues 3"       "Purples 2"     "Purples 3"     "Reds 2"       
##  [17] "Reds 3"        "Greens 2"      "Greens 3"      "Oslo"         
##  [21] "Purple-Blue"   "Red-Purple"    "Red-Blue"      "Purple-Orange"
##  [25] "Purple-Yellow" "Blue-Yellow"   "Green-Yellow"  "Red-Yellow"   
##  [29] "Heat"          "Heat 2"        "Terrain"       "Terrain 2"    
##  [33] "Viridis"       "Plasma"        "Inferno"       "Rocket"       
##  [37] "Mako"          "Dark Mint"     "Mint"          "BluGrn"       
##  [41] "Teal"          "TealGrn"       "Emrld"         "BluYl"        
##  [45] "ag_GrnYl"      "Peach"         "PinkYl"        "Burg"         
##  [49] "BurgYl"        "RedOr"         "OrYel"         "Purp"         
##  [53] "PurpOr"        "Sunset"        "Magenta"       "SunsetDark"   
##  [57] "ag_Sunset"     "BrwnYl"        "YlOrRd"        "YlOrBr"       
##  [61] "OrRd"          "Oranges"       "YlGn"          "YlGnBu"       
##  [65] "Reds"          "RdPu"          "PuRd"          "Purples"      
##  [69] "PuBuGn"        "PuBu"          "Greens"        "BuGn"         
##  [73] "GnBu"          "BuPu"          "Blues"         "Lajolla"      
##  [77] "Turku"         "Hawaii"        "Batlow"        "Blue-Red"     
##  [81] "Blue-Red 2"    "Blue-Red 3"    "Red-Green"     "Purple-Green" 
##  [85] "Purple-Brown"  "Green-Brown"   "Blue-Yellow 2" "Blue-Yellow 3"
##  [89] "Green-Orange"  "Cyan-Magenta"  "Tropic"        "Broc"         
##  [93] "Cork"          "Vik"           "Berlin"        "Lisbon"       
##  [97] "Tofino"        "ArmyRose"      "Earth"         "Fall"         
## [101] "Geyser"        "TealRose"      "Temps"         "PuOr"         
## [105] "RdBu"          "RdGy"          "PiYG"          "PRGn"         
## [109] "BrBG"          "RdYlBu"        "RdYlGn"        "Spectral"     
## [113] "Zissou 1"      "Cividis"       "Roma"

heatmapEnrichment(pbmc_small, 
                  assay = "escape.ssGSEA",
                  palette = "Spectral") 

Alternatively, we can add an additional layer to the ggplot object that is returned by the visualizations using something like scale_fill_gradientn() for continuous values or scale_fill_manual() for the categorical variables.

heatmapEnrichment(sce.pbmc, 
                  group.by = "ident",
                  assay = "escape.UCell") + 
  scale_fill_gradientn(colors = rev(brewer.pal(11, "RdYlBu"))) 

geyserEnrichment

We can also focus on individual gene sets - one approach is to use geyserEnrichment(). Here individual cells are plotted along the Y-axis with graphical summary where the central dot refers to the median enrichment value and the thicker/thinner lines demonstrate the interval summaries referring to the 66% and 95%.

geyserEnrichment(pbmc_small,
                 assay = "escape.ssGSEA",
                 gene.set.use = "T1-Interferon")

To show the additional parameters that appear in visualizations of individual enrichment gene sets - we can reorder the groups by the mean of the gene set using order.by = “mean”.

geyserEnrichment(pbmc_small,
                 assay = "escape.ssGSEA",
                 gene.set.use = "T1-Interferon",
                 order.by = "mean")

What if we had 2 separate samples or groups within the data? Another parameter we can use is facet.by to allow for direct visualization of an additional variable.

geyserEnrichment(pbmc_small,
                 assay = "escape.ssGSEA",
                 gene.set.use = "T1-Interferon",
                 facet.by = "groups")

Lastly, we can select the way the color is applied to the plot using the color.by parameter. Here we can set it to the gene set of interest “HALLMARK-INTERFERON-GAMMA-RESPONSE”.

geyserEnrichment(pbmc_small,
                 assay = "escape.ssGSEA",
                 gene.set.use = "T1-Interferon",
                 color.by  = "T1-Interferon")

ridgeEnrichment

Similar to the geyserEnrichment() the ridgeEnrichment() can display the distribution of enrichment values across the selected gene set. The central line is at the median value for the respective grouping.

ridgeEnrichment(sce.pbmc,
                assay = "escape.UCell",
                gene.set.use = "T2_Interferon")

We can get the relative position of individual cells along the x-axis using the add.rug parameter.

ridgeEnrichment(sce.pbmc,
                assay = "escape.UCell",
                gene.set.use = "T2_Interferon",
                add.rug = TRUE,
                scale = TRUE)

splitEnrichment

Another distribution visualization is a violin plot, which we separate and directly compare using a binary classification. Like ridgeEnrichment(), this allows for greater use of categorical variables. For splitEnrichment(), the output will be two halves of a violin plot based on the split.by parameter with a central boxplot with the relative distribution across all samples.

splitEnrichment(pbmc_small,
                assay = "escape.ssGSEA",
                gene.set.use = "Lipid-mediators",
                split.by = "groups")

If selecting a split.by variable with more than 2 levels, splitEnrichment() will convert the violin plots to dodge.

splitEnrichment(pbmc_small,
                assay = "escape.ssGSEA",
                gene.set.use = "Lipid-mediators",
                split.by = "ident",
                group.by = "groups")

gseaEnrichment

gseaEnrichment() reproduces the two-panel GSEA graphic from Subramanian et al. (2005): * Panel A – the running enrichment score (RES) as you “walk” down the ranked list. * Panel B – a rug showing exact positions of each pathway gene.

It works on escape’s per-cell ranks, but collapses them across cells with a summary statistic (summary.fun = “mean” by default).

How it works:

1. Rank all genes in each group by summary.fun of expression/statistic.
2. Perform the weighted Kolmogorov–Smirnov walk: +w when the next gene is in the set, −1/(N − NG) otherwise.
3. ES = maximum signed deviation; permutation on gene labels (or phenotypes) to derive NES and p.

gseaEnrichment(pbmc_small,
               gene.set.use = "T2_Interferon",
               gene.sets    = escape.gene.sets,
               group.by     = "ident",
               summary.fun  = "mean",
               nperm        = 1000)

densityEnrichment

densityEnrichment() is a method to visualize the mean rank position of the gene set features along the total feature space by group. Instead of the classic GSEA running-score, it overlays kernel-density traces of the gene ranks (1 = most highly expressed/ranked gene) for every group or cluster. High densities at the left-hand side mean the pathway is collectively up-regulated; peaks on the right imply down-regulation.

Anatomy of the plot

1. X-axis – gene rank (1 … N). Left = top-ranked genes.
   
2. Y-axis – density estimate (area under each curve = 1).
   
3. One coloured line per level of group.by – default is Seurat/SCE cluster.

densityEnrichment(pbmc_small, 
                  gene.set.use = "T2_Interferon", 
                  gene.sets = escape.gene.sets)

scatterEnrichment

It may be advantageous to look at the distribution of multiple gene sets - here we can use scatterEnrichment() for a 2 gene set comparison. The color values are based on the density of points determined by the number of neighbors, similar to the Nebulosa R package. We just need to define which gene set to plot on the x.axis and which to plot on the y.axis.

scatterEnrichment(pbmc_small, 
                  assay = "escape.ssGSEA",
                  x.axis = "T2-Interferon",
                  y.axis = "Lipid-mediators")

The scatter plot can also be converted into a hexbin, another method for summarizing the individual cell distributions along the x and y axis, by setting style = “hex”.

scatterEnrichment(sce.pbmc, 
                  assay = "escape.UCell",
                  x.axis = "T2_Interferon",
                  y.axis = "Lipid_mediators",
                  style = "hex")

Principal Component Analysis (PCA)

escape has its own PCA function performPCA() which will work on a single-cell object or a matrix of enrichment values. This is specifically useful for downstream visualizations as it stores the eigenvalues and rotations. If we want to look at the relative contribution to overall variance of each component or a Biplot-like overlay of the individual features, use performPCA().

Alternatively, other PCA-based functions like Seurat’s RunPCA() or scater’s ``runPCA() can be used. These functions are likely faster and would be ideal if we have a larger number of cells and/or gene sets.

pbmc_small <- performPCA(pbmc_small, 
                         assay = "escape.ssGSEA",
                         n.dim = 1:10)

escape has a built in method for plotting PCA pcaEnrichment() that functions similarly to the scatterEnrichment() function where x.axis and y.axis are the components to plot.

pcaEnrichment(pbmc_small, 
              dimRed = "escape.PCA",
              x.axis = "PC1",
              y.axis = "PC2")

pcaEnrichment() can plot additional information on the principal component analysis.

add.percent.contribution will add the relative percent contribution of the x and y.axis to total variability observed in the PCA.

display.factors will overlay the magnitude and direction that the features/gene sets contribute to the selected components. The number of gene sets is determined by number.of.factors. This can assist in understanding the underlying differences in enrichment across different cells.

pcaEnrichment(pbmc_small, 
              dimRed = "escape.PCA",
              x.axis = "PC1",
              y.axis = "PC2",
              add.percent.contribution = TRUE,
              display.factors = TRUE,
              number.of.factors = 10)

Precomputed Rank Lists

Functional enrichment is not limited to per-cell scores. Many workflows start with differential-expression (DE) statistics (e.g. Seurat’s FindMarkers(), DESeq2’s results(), edgeR’s topTags()). Those produce a ranked gene list that can be fed into a classical Gene-Set Enrichment Analysis (GSEA).

DEG.markers <- FindMarkers(pbmc_small, 
                           ident.1 = "0", 
                           ident.2 = "1")

GSEA.results <- enrichIt(input.data = DEG.markers, 
                         gene.sets = escape.gene.sets, 
                         ranking_fun = "signed_log10_p")               

head(GSEA.results)

##            pathway      pval      padj    log2err         ES        NES  size
##             <char>     <num>     <num>      <num>      <num>      <num> <int>
## 1:   T2_Interferon 0.2478541 0.8083678 0.08336341 -0.8228147 -1.1007827    10
## 2:        M1.Macro 0.3233471 0.8083678 0.06767604 -0.7714736 -1.0683506    15
## 3: Proinflammatory 0.7155963 0.9487459 0.15740290  0.4002094  0.6185734     8
## 4: Lipid_mediators 0.8472554 0.9487459 0.02921031 -0.6849303 -0.8766974     6
## 5:  CD8_Activation 0.9487459 0.9487459 0.01770381 -0.5167583 -0.6855810     9
##                                                                                leadingEdge
##                                                                                     <char>
## 1:           HLA-DRB1;HLA-DPB1;HLA-DPA1;HLA-DRA;HLA-DRB5;HLA-DQB1;HLA-DQA1;HLA-DMB;HLA-DMA
## 2: HLA-DRB1;HLA-DPB1;HLA-DPA1;HLA-DRA;HLA-DRB5;CD14;HLA-DQB1;HLA-DQA1;HLA-DMB;HLA-DMA;IL1B
## 3:                                                                               GZMA;PRF1
## 4:                                                                    CTSS;CFD;LGALS3;IL1B
## 5:                                                                                 HLA-DRA
##    geneRatio
##        <num>
## 1: 0.9000000
## 2: 0.7333333
## 3: 0.2500000
## 4: 0.6666667
## 5: 0.1111111

## (1) Bar plot –20 most significant per database
enrichItPlot(GSEA.results)               

## (2) Dot plot – coloured by –log10 padj, sized by core-hits
enrichItPlot(GSEA.results, "dot", top = 10) 

## (3) C-net plot – network of pathways ↔ leading-edge genes
enrichItPlot(GSEA.results, "cnet", top = 5) 

Differential Enrichment

Differential enrichment analysis can be performed similar to differential gene expression analysis. For the purposes of finding the differential enrichment values, we can first normalize the enrichment values for the ssGSEA calculations. Notice here, we are using make.positive = TRUE in order to adjust any negative values. This is a particular issue when it comes to ssGSEA and GSVA enrichment scores.

pbmc_small <- performNormalization(pbmc_small, 
                                   assay = "escape.ssGSEA", 
                                   gene.sets = escape.gene.sets,
                                   make.positive = TRUE)

all.markers <- FindAllMarkers(pbmc_small, 
                              assay = "escape.ssGSEA_normalized", 
                              min.pct = 0,
                              logfc.threshold = 0)

head(all.markers)

##                        p_val avg_log2FC pct.1 pct.2    p_val_adj cluster
## T2-Interferon   2.210898e-09 -0.4499809     1     1 1.989808e-08       0
## CD8-Activation  2.791543e-06  0.9948584     1     1 2.512388e-05       0
## M1.Macro        9.243851e-06 -0.3021382     1     1 8.319466e-05       0
## Cytolytic       1.104433e-04  1.9352033     1     1 9.939895e-04       0
## T1-Interferon   4.895518e-04 -0.4474723     1     1 4.405966e-03       0
## Lipid-mediators 6.294119e-03 -0.2174715     1     1 5.664707e-02       0
##                            gene
## T2-Interferon     T2-Interferon
## CD8-Activation   CD8-Activation
## M1.Macro               M1.Macro
## Cytolytic             Cytolytic
## T1-Interferon     T1-Interferon
## Lipid-mediators Lipid-mediators