Introduction

GSVA provides now specific support for single-cell data in the algorithm that runs through the gsvaParam() parameter constructor, and originally described in the publication by Hänzelmann et al. (2013). At the moment, this specific support consists of the following features:

• The input expression data can be stored in different types of data containers prepared to store sparse single-cell data. These types of sparse data containers can be broadly categorized in those that only store the expression values, and those that may store additional row and column metadata. The currently available value-only containers for input are dgCMatrix, SVT_SparseArray, and DelayedMatrix. The currently available container for single-cell data that allows one to input additional row and column metadata is a SingleCellExperiment object.
• While the input single-cell data is always sparse, the output of enrichment scores will be always dense, and therefore, the container storing those scores will be different from the input data, typically a matrix or a dense DelayedMatrix object. The latter will be particularly used when the total number of values exceeds 2^31, which is the largest 32-bit standard integer value in R.
• By default, when the input expression data is stored in a sparse data container, as it typically happens with single-cell data, then a slightly a slightly modified GSVA algorithm will run, if GSVA is the choice of algorithm, by which nonzero values are treated differently from zero values, leading to slightly different results than those obtained by applying the classical GSVA algorithm. If we set the parameter sparse=FALSE in the call to gsvaParam(), the classical GSVA algorithm will be used, which for a typical single-cell data set will result in longer running times and larger memory consumption than running it in the default sparse regime for this type of data.
• The GSVA algorithm can be run either at once through a call to gsva() with a parameter object or in two steps: (1) row normalization and column rank transformation with gsvaRanks(); and (2) column enrichment scores calculation with gsvaScores(). Splitting the GSVA algorithm into these two steps allows one to reuse the output of the first step, which is independent of the gene sets, to calculate enrichment scores for different collections of gene sets, without having to repeat the first step.

In what follows, we will illustrate the use of GSVA on a publicly available single-cell transcriptomics data set of peripheral blood mononuclear cells (PBMCs) published by Zheng et al. (2017).

Import data

We import the PBMC data using the TENxPBMCData package, as a SingleCellExperiment object, defined in the SingleCellExperiment package.

library(SingleCellExperiment)
library(TENxPBMCData)

sce <- TENxPBMCData(dataset="pbmc4k")
sce
class: SingleCellExperiment 
dim: 33694 4340 
metadata(0):
assays(1): counts
rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
  ENSG00000268674
rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
colnames: NULL
colData names(11): Sample Barcode ... Individual Date_published
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

Quality assessment and pre-processing

Here, we perform a quality assessment and pre-processing steps using the package scuttle (McCarthy et al. 2017). We start identifying mitochondrial genes.

library(scrapper)

is_mito <- grepl("^MT-", rowData(sce)$Symbol_TENx)
table(is_mito)
is_mito
FALSE  TRUE 
33681    13 

Calculate quality control (QC) metrics and filter out low-quality cells.

sce <- quickRnaQc.se(sce, subsets=list(mito=is_mito))
sce <- sce[, sce$keep]
dim(sce)
[1] 33694  4147

We filter out genes that are expressed in less than 1% of the cells.

cellsxgene <- rowSums(counts(sce) > 0)
sce <- sce[cellsxgene > floor(ncol(sce)*0.01), ]
dim(sce)
[1] 10799  4147

Calculate library size factors and normalized units of expression in logarithmic scale.

sce <- normalizeRnaCounts.se(sce, size.factors=sce$sum)
assayNames(sce)
[1] "counts"    "logcounts"

Annotate cell types using GSVA

Here, we illustrate how to annotate cell types in the PBMC data using GSVA.

library(GSEABase)
library(GSVA)

fname <- file.path(system.file("extdata", package="GSVAdata"),
                   "pbmc_cell_type_gene_set_signatures.gmt.gz")
gsets <- readGMT(fname)
gsets
GeneSetCollection
  names: B_CELLS_MEMORY, B_CELLS_NAIVE, ..., T_CELLS_REGULATORY_TREGS (22 total)
  unique identifiers: AIM2, BANK1, ..., SKAP1 (248 total)
  types in collection:
    geneIdType: SymbolIdentifier (1 total)
    collectionType: NullCollection (1 total)

gsvaAnnotation(sce) <- ENSEMBLIdentifier("org.Hs.eg.db")
gsvaAnnotation(sce)
geneIdType: ENSEMBL (org.Hs.eg.db)

gsvapar <- gsvaParam(sce, gsets)
gsvapar
A GSVA::gsvaParam object
expression data:
  class: SingleCellExperiment 
  dim: 10799 4147 
  metadata(2): qc annotation
  assays(2): counts logcounts
  rownames(10799): ENSG00000279457 ENSG00000228463 ... ENSG00000273748
    ENSG00000278817
  rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
  colnames: NULL
  colData names(16): Sample Barcode ... keep sizeFactor
  reducedDimNames(0):
  mainExpName: NULL
  altExpNames(0):
using assay: logcounts
using annotation:
  geneIdType: ENSEMBL (org.Hs.eg.db)
gene sets:
  GeneSetCollection
    names: B_CELLS_MEMORY, B_CELLS_NAIVE, ..., T_CELLS_REGULATORY_TREGS (22 total)
    unique identifiers: AIM2, BANK1, ..., SKAP1 (248 total)
    types in collection:
      geneIdType: SymbolIdentifier (1 total)
      collectionType: NullCollection (1 total)
gene set size: [1, Inf]
nonzero values: less than 2^31 (INT_MAX)
ondisk:  auto 
kcdf: auto
kcdfNoneMinSampleSize: 200
tau: 1
maxDiff: TRUE
absRanking: FALSE
sparse:  TRUE 
checkNA: auto
missing data: didn't check
filterRows:  TRUE 

gsvaranks <- gsvaRanks(gsvapar)
gsvaranks
A GSVA::gsvaRanksParam object
expression data:
  class: SingleCellExperiment 
  dim: 10799 4147 
  metadata(2): qc annotation
  assays(3): counts logcounts gsvaranks
  rownames(10799): ENSG00000279457 ENSG00000228463 ... ENSG00000273748
    ENSG00000278817
  rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
  colnames: NULL
  colData names(16): Sample Barcode ... keep sizeFactor
  reducedDimNames(0):
  mainExpName: NULL
  altExpNames(0):
using assay: gsvaranks
using annotation:
  geneIdType: ENSEMBL (org.Hs.eg.db)
gene sets:
  GeneSetCollection
    names: B_CELLS_MEMORY, B_CELLS_NAIVE, ..., T_CELLS_REGULATORY_TREGS (22 total)
    unique identifiers: AIM2, BANK1, ..., SKAP1 (248 total)
    types in collection:
      geneIdType: SymbolIdentifier (1 total)
      collectionType: NullCollection (1 total)
gene set size: [1, Inf]
nonzero values: less than 2^31 (INT_MAX)
ondisk:  auto 
kcdf: auto
kcdfNoneMinSampleSize: 200
tau: 1
maxDiff: TRUE
absRanking: FALSE
sparse:  TRUE 
checkNA: auto
missing data: didn't check
filterRows:  TRUE 

es <- gsvaScores(gsvaranks)
es
class: SingleCellExperiment 
dim: 22 4147 
metadata(1): qc
assays(1): es
rownames(22): B_CELLS_MEMORY B_CELLS_NAIVE ... T_CELLS_GAMMA_DELTA
  T_CELLS_REGULATORY_TREGS
rowData names(1): gs
colnames: NULL
colData names(16): Sample Barcode ... keep sizeFactor
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

geneSets(gsvaranks) <- geneSets(gsvapar)[1:2]
es2 <- gsvaScores(gsvaranks)

Using GSVA scores to assign cell types

Following Amezquita et al. (2020), and some of the steps described in “Chapter 5 Clustering” of the first version of the OSCA book, we use GSVA scores to build a nearest-neighbor graph of the cells using the function buildSNNGraph() from the scran package (Lun et al. 2016). The parameter k in the call to buildSNNGraph() specifies the number of nearest neighbors to consider during graph construction, and here we set k=20 because it leads to a number of clusters close to the expected number of cell types.

library(bluster)

g <- makeSNNGraph(t(assay(es)), k=20)

Second, we use the function cluster_walktrap() from the igraph package (Csardi and Nepusz 2006), to cluster cells by finding densely connected subgraphs. We store the resulting vector of cluster indicator values into the sce object using the function colLabels().

library(igraph)

colLabels(es) <- factor(cluster_walktrap(g)$membership)
table(colLabels(es))

   1    2    3    4    5    6    7    8 
 485  345  714 1047  600  220  197  539 

Similarly to Diaz-Mejia et al. (2019), we apply a simple cell type assignment algorithm, which consists of selecting at each cell the gene set with highest GSVA score, tallying the selected gene sets per cluster, and assigning to the cluster the most frequent gene set, storing that assignment into the sce object with the function colLabels().

whmax <- apply(assay(es), 2, which.max)
gsxlab <- split(rownames(es)[whmax], colLabels(es))
gsxlab <- names(sapply(sapply(gsxlab, table), which.max))
colLabels(es) <- factor(gsub("[0-9]\\.", "", gsxlab))[colLabels(es)]
table(colLabels(es))

   B_CELLS_MEMORY         MONOCYTES  NK_CELLS_RESTING T_CELLS_CD4_NAIVE 
              600              1059               220              1783 
      T_CELLS_CD8 
              485 

We can visualize the cell type assignments by projecting cells dissimilarity in two dimensions with a principal components analysis (PCA) on the GSVA scores, and coloring cells using the previously assigned clusters.

library(RColorBrewer)

res <- prcomp(assay(es))
varexp <- res$sdev^2 / sum(res$sdev^2)
nclusters <- nlevels(colLabels(es))
hmcol <- colorRampPalette(brewer.pal(nclusters, "Set1"))(nclusters)
par(mar=c(4, 5, 1, 1))
plot(res$rotation[, 1], res$rotation[, 2], col=hmcol[colLabels(es)], pch=19,
     xlab=sprintf("PCA 1 (%.0f%%)", varexp[1]*100),
     ylab=sprintf("PCA 2 (%.0f%%)", varexp[2]*100),
     las=1, cex.axis=1.2, cex.lab=1.5)
legend("topright", gsub("_", " ", levels(colLabels(es))), fill=hmcol, inset=0.01)

Cell type assignments of PBMC scRNA-seq data, based on GSVA scores.

Finally, if we want to better understand why a specific cell type is annotated to a given cell, we can use the gsvaEnrichment() function, which will show a GSEA enrichment plot. This function takes as input the output of gsvaRanks(), a given column (cell) in the input singl-cell data, and a given gene set. In Figure @ref(fig:gsvaenrichment) below, we show such a plot for the first cell annotated to the monocytes cell type.

firstmonocytecell <- which(colLabels(es) == "MONOCYTES")[1]
par(mar=c(4, 5, 1, 1))
gsvaEnrichment(gsvaranks, column=firstmonocytecell, geneSet="MONOCYTES",
               cex.axis=1.2, cex.lab=1.5, plot="ggplot")

GSVA enrichment plot of the MONOCYTES gene set in the expression profile of the first cell annotated to that cell type.

In the previous call to gsvaEnrichment() we used the argument plot="ggplot" to produce a plot with the ggplot2 package. By default, if we call gsvaEnrichment() interactively, it will produce a plot using “base R”, but either when we do it non-interactively, or when we set plot="no" it will return a data.frame object with the enrichment data.

Stability of GSVA scores with sparse data

This version of GSVA (>= 2.6.4) fixes a known problem with the stability of GSVA scores with sparse data. More concretely, for a given gene set, its genes with a zero count expression value in a given column/cell, would occupy arbitrary positions in the gene ranking of that column/cell. This, and newer versions of GSVA, place those genes evenly distributed after the nonzerovalues along the ranking. We can see this happening with the MONOCYTES gene set in Figure @ref(fig:gsvaenrichment), where the genes of this gene set with zero expression values are evenly distributed after the left vertical dashed line that in this plot indicates the largest positive leading edge of the enrichment score. We can also verify that those genes have zero expression values in the first cell annotated to the MONOCYTES cell type, as follows.

edat <- gsvaEnrichment(gsvaranks, column=firstmonocytecell,
                       geneSet="MONOCYTES", plot="no")
unname(counts(sce)[edat$gsetidx, firstmonocytecell][order(edat$gsetrnk)])
 [1] 4 5 6 6 6 1 1 5 2 1 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0

Note that the nonzero count values are not perfectly sorted in decreasing order, because rankings are calculated on a expression statistic after column and row normalization. Using the diff() function, which calculates differences between consecutive elements of a vector, we can also verify how exactly these zeros of the gene set are evenly spaced.

diff(sort(edat$gsetrnk))
 [1]  83  27 148  62 197  48   8 109  36  66  59   3 133 887 888 887 887 887 888
[20] 887 887 887 888 887

Finally, we can verify that this strategy results in stable GSVA scores, by running the GSVA algorithm on a random permutation of the rows of the sce object, and checking that the resulting GSVA scores are nearly identical to those obtained with the original sce object.

gsvapar <- gsvaParam(sce[sample(1:nrow(sce)), ], gsets, verbose=FALSE)
es2 <- gsva(gsvapar, verbose=FALSE)
all.equal(assay(es), assay(es2))
[1] "Mean relative difference: 0.0006704369"

The tiny differences in the enrichment scores (less than 1e-3 in these data) are due to the nonzero values in gene sets that are tied in the ranking of the internally calculated expression statistic.

Benchmarking

We are still benchmarking and testing this version of GSVA for single-cell data. If you encounter problems or have suggestions, do not hesitate to contact us by opening an issue in the GSVA GitHub repo.