Example data set

First, load the SingleCellExperiment object of 10X Genomics 2,700 peripheral blood mononuclear cells (PBMC) preprocessed in a previous vignette.

library(SingleCellExperiment)
sce <- readRDS("pbmc3k_tutorial.sce.rds")
sce

## class: SingleCellExperiment 
## dim: 32738 2638 
## metadata(0):
## assays(1): counts
## rownames(32738): ENSG00000243485 ENSG00000237613 ...
##   ENSG00000215616 ENSG00000215611
## rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
## colnames: NULL
## colData names(12): Sample Barcode ... Date_published seurat.ident
## reducedDimNames(2): PCA TSNE
## spikeNames(0):

Run paired t-tests to identify markers

First, we normalize and log-transform the count data.

library(scran)
library(scater)
sce <- computeSumFactors(sce)
sce <- normalize(sce)
sce

## class: SingleCellExperiment 
## dim: 32738 2638 
## metadata(1): log.exprs.offset
## assays(2): counts logcounts
## rownames(32738): ENSG00000243485 ENSG00000237613 ...
##   ENSG00000215616 ENSG00000215611
## rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
## colnames: NULL
## colData names(12): Sample Barcode ... Date_published seurat.ident
## reducedDimNames(2): PCA TSNE
## spikeNames(0):

In this example, we look for positive markers. Namely, up-regulated genes in each cluster relative to all other clusters.

library(scran)
out <- findMarkers(sce, clusters=sce$seurat.ident, direction="up", lfc=1)
out

## DataFrameList of length 8
## names(8): 0 1 2 3 4 5 6 7

The result out is a list of DataFrame. For each cluster, the log fold-change against every other cluster is reported, along with a combined p-value and FDR.

colnames(out[[1]])

##  [1] "Top"     "p.value" "FDR"     "logFC.1" "logFC.2" "logFC.3" "logFC.4"
##  [8] "logFC.5" "logFC.6" "logFC.7"

Obtain gene signatures for relevant cell types

library(xCellData)
library(unisets)
xcell <- xCellData()
xcell

## Sets with 20803 relations between 5079 elements and 489 sets
##             element          set
##         <character>  <character>
##     [1]        C1QA   aDC_HPCA_1
##     [2]        C1QB   aDC_HPCA_1
##     [3]       CCL13   aDC_HPCA_1
##     [4]       CCL17   aDC_HPCA_1
##     [5]       CCL19   aDC_HPCA_1
##     ...         ...          ...
## [20799]       IL2RA Tregs_HPCA_3
## [20800]       KCNA2 Tregs_HPCA_3
## [20801]       LAIR2 Tregs_HPCA_3
## [20802]      MCF2L2 Tregs_HPCA_3
## [20803]        RGS1 Tregs_HPCA_3
## -----------
## elementInfo: IdVector with 0 metadata
##     setInfo: IdVector with 1 metadata (source)

library(org.Hs.eg.db)
TERM2GENE <- as.data.frame(xcell)[, c("set", "element")]
TERM2GENE$element <- mapIds(org.Hs.eg.db, TERM2GENE$element, "ENSEMBL", "SYMBOL")

TERM2NAME <- data.frame(
    id = ids(setInfo(xcell)),
    name = rep("", nSets(xcell))
)

Run the analysis

In this example, we apply Gene Set Enrichment Analysis (GSEA) on the lists of genes ranked by decreasing detection average log-fold change between each cluster and all other clusters. Below, we display the 2 most signicantly enriched cell type gene signature for each cluster. Note that for some clusters, none of the signatures reached signifcance.

library(clusterProfiler)
runGSEA <- function(clusterName, threshold=0, top=10) {
    x <- out[[clusterName]]
    x_lfc_colnames <- grep("^logFC\\.", colnames(x))
    x_rowmeans <- rowMeans(as.matrix(x[, x_lfc_colnames]))
    x_rowmeans <- sort(x_rowmeans, decreasing = TRUE)
    x_rowmeans <- x_rowmeans[x_rowmeans > threshold]
    y = GSEA(x_rowmeans, TERM2GENE=TERM2GENE, TERM2NAME=TERM2NAME, pvalueCutoff=1, verbose = FALSE)
    y <- as.data.frame(y)
    y <- head(y, top)
    y <- cbind(cluster = clusterName, y)
    y
}
resList <- lapply(levels(sce$seurat.ident), runGSEA, threshold=0, top=2)
res <- do.call(rbind, resList)
res <- as.data.frame(res)
show_columns <- setdiff(colnames(res), "core_enrichment")
kable(res[, show_columns])

See also

• Applying Gene Signatures to the Seurat PBMC 3k Tutorial
• Learning Gene Signatures from the Seurat PBMC 3k Tutorial