Read gene sets in GMT format

First, we fetch a collection of 22 leukocyte gene set signatures, containing a total 547 genes, which should help to distinguish among 22 mature human hematopoietic cell type populations isolated from peripheral blood or in vitro culture conditions, including seven T cell types: naïve and memory B cells, plasma cells, NK cell, and myeloid subsets. These gene sets have been used in the benchmarking publication by Diaz-Mejia et al. (2019), and were originally compiled by the CIBERSORT developers, where they called it the LM22 signature (Newman et al. 2015). The LM22 signature is stored in the GSVAdata experiment data package as a compressed text file in GMT format, which can be read into R using the readGMT() function from the GSVA package, which will return the gene sets, by default, into a GeneSetCollection object, defined in the GSEABase package. This default argument can be changed to return the gene sets into a base list object by setting valueType="list" in the call to readGMT().

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)

Add gene identifier type metadata

Note that while gene identifers in the sce object correspond to Ensembl stable identifiers (ENSG...), the gene identifiers in the gene sets are HGNC gene symbols. This, in principle, precludes matching directly what gene in the single-cell data object sce corresponds to what gene set in the GeneSetCollection object gsets. However, the GSVA package can do that matching as long as the appropriate metadata is present in both objects.

In the case of a GeneSetCollection object, its geneIdType metadata slot stores the type of gene identifier. In the case of a SingleCellExperiment object, such as the previous sce object, such metadata is not present. However, using the function gsvaAnnotation() from the GSVA package, and the helper function ENSEMBLIdentifier() from the GSEABase package, we add such metadata to the sce object as follows.

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

Build parameter object

We first build a parameter object using the function gsvaParam(). By default, the expression values in the logcounts assay will be selected for downstream analysis.

gsvapar <- gsvaParam(sce, gsets)
gsvapar
class: GSVA::gsvaParam
expression data dim: 10799 4147
number of gene sets: 22
details: use 'details(object)'

Calculate GSVA scores

While at this point, we could already run the entire GSVA algorithm with a call to the gsva(gsvapar) function. We show here how to do it in three steps. First we calculate row-normalized expression values using the function gsvaRowNorm(), which if, as in this example, the given input is a SingleCellExperiment object, then the output will be the same object with an additional assay called gsvarownr containing the row-normalized expression values.

gsvarownorm <- gsvaRowNorm(gsvapar)
gsvarownorm
class: SingleCellExperiment 
dim: 10799 4147 
metadata(3): qc annotation gsvaParam
assays(3): counts logcounts gsvarownr
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):
assayNames(gsvarownorm)
[1] "counts"    "logcounts" "gsvarownr"

Second, we calculate GSVA column rank values using the function gsvaColRanks(), which takes as input the output of gsvaRowNorm(), and returns the column rank values in a new assay called gsvaranks, if the input is a SingleCellExperiment object.

gsvacolranks <- gsvaColRanks(gsvarownorm)
gsvacolranks
class: SingleCellExperiment 
dim: 10799 4147 
metadata(3): qc annotation gsvaParam
assays(4): counts logcounts gsvarownr 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):
assayNames(gsvacolranks)
[1] "counts"    "logcounts" "gsvarownr" "gsvaranks"

Third, we finally calculate the GSVA scores using the output of gsvaColRanks() as input to the function gsvaColScores(). By default, this function will calculate the scores for all gene sets specified in the input parameter object given in the call to gsvaRowNorm().

es <- gsvaColScores(gsvacolranks)
es
class: SingleCellExperiment 
dim: 22 4147 
metadata(2): qc gsvaParam
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):

However, we could calculate the scores for another collection of gene sets, without having to calculate the column ranks again, by giving this other collection of gene sets as second argument to the call to gsvaColScores().

es2 <- gsvaColScores(gsvacolranks, alternative_gsets)

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 makeSNNGraph() from the bluster package. The parameter k in the call to makeSNNGraph() 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 
 919 1081 1017  589  205  336 

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_NAIVE         MONOCYTES  NK_CELLS_RESTING T_CELLS_CD4_NAIVE 
              589              1017               205              2000 
      T_CELLS_CD8 
              336 

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)
mask <- colLabels(es) == "NK_CELLS_RESTING"
points(res$rotation[mask, 1], res$rotation[mask, 2], ## show the overlap better
       col=hmcol[colLabels(es)[mask]], pch=19)
legend("bottomright", 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 single-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(gsvacolranks, column=firstmonocytecell, geneSet="MONOCYTES",
               cex.axis=1.2, cex.lab=1.5, plot="ggplot")

GSVA enrichment plot of the EOSINOPHILS 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.