Cell Type Annotation

Introduction

This is a brief tutorial on automatic cell type annotation of single-cell RNA sequencing (scRNA-seq) data. The primary dataset used here is a set of 2,700 Peripheral Blood Mononuclear Cells (PBMCs) processed using the standard Seurat workflow as demonstrated in Seurat’s clustering tutorial. It is a reasonably small dataset with well-established cell types that is commonly used in scRNA-seq benchmarking studies.

Load data

This tutorial includes some Bioconductor dependencies. Before proceeding, confirm that Bioconductor is installed and its version is at least 3.10.

if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::version()

[1] ‘3.10’

If this is not the case, remove all versions of BiocVersion with remove.packages("BiocVersion"). Then update Bioconductor packages using BiocManager::install().

Since we are using a Seurat object, load Seurat and related packages.

if (!require("Seurat")) {
   BiocManager::install("Seurat", update = FALSE)
   library(Seurat)
}
library(ggplot2)
library(cowplot)

Load other relevant packages.

if (!require("dplyr")) {
   BiocManager::install("dplyr", update = FALSE)
   library(dplyr)
}
if (!require("stringr")) {
   BiocManager::install("stringr", update = FALSE)
   library(stringr)
}

Load the PBMC dataset using the SeuratData package.

if (!require("SeuratData")) {
   BiocManager::install("satijalab/seurat-data", update = FALSE)
   library(SeuratData)
}
InstallData("pbmc3k")
data("pbmc3k")
pbmc3k

An object of class Seurat 
13714 features across 2700 samples within 1 assay 
Active assay: RNA (13714 features)

The dataset includes both the raw and the processed versions.

pbmc3k.final

An object of class Seurat 
13714 features across 2638 samples within 1 assay 
Active assay: RNA (13714 features)
 2 dimensional reductions calculated: pca, umap

Check the original Seurat cell type labels overlaid onto the UMAP visualization. See the original Seurat guided clustering tutorial for details on how the data was processed and the labels were assigned.

DimPlot(pbmc3k.final, reduction = "umap") +
  scale_color_brewer(palette = "Set1")

Annotation (using SingleR)

Load SingleR.

if (!require("SingleR")) {
   BiocManager::install("SingleR", update = FALSE)
   library(SingleR)
}

SingleR expects the input as a matrix or a SummarizedExperiment object.

exp_mat = GetAssayData(pbmc3k.final, assay = "RNA", slot = "data")
exp_mat = as.matrix(exp_mat)
dim(exp_mat)

[1] 13714  2638

HPCA main cell types

The SingleR package provides normalized expression values and cell types labels based on bulk RNA-seq, microarray, and single-cell RNA-seq data from several different datasets. See the SingleR vignette for the description.

We can try the Human Primary Cell Atlas dataset as the reference. It provides normalized expression values for 713 microarray samples that have been assigned to one of 37 main cell types and 157 subtypes.

singler_se = HumanPrimaryCellAtlasData()

It is possible that this fails with a No internet connection using 'localHub=TRUE' error. This may be resolved by running ExperimentHub::setExperimentHubOption("PROXY", "http://127.0.0.1:10801").

singler_se

class: SummarizedExperiment 
dim: 19363 713 
metadata(0):
assays(1): logcounts
rownames(19363): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
rowData names(0):
colnames(713): GSM112490 GSM112491 ... GSM92233 GSM92234
colData names(2): label.main label.fine

Restrict to common genes between the test and reference datasets.

common_genes = intersect(rownames(exp_mat), rownames(singler_se))
common_genes = sort(common_genes)
exp_common_mat = exp_mat[common_genes, ]
singler_se = singler_se[common_genes, ]
length(common_genes)

[1] 11375

Perform SingleR annotation.

singler_pred = SingleR(
  test = exp_common_mat,
  ref = singler_se,
  labels = singler_se$label.main
)

Each row of the output data frame contains prediction results for a single cell.

head(as.data.frame(singler_pred))

SingleR provides a method to display the scores for all cells across all reference labels to inspect the confidence of the predicted labels across the dataset.

plotScoreHeatmap(singler_pred, show.labels = TRUE, fontsize = 8)

As expected, the non-immune populations have very low scores in these cells.

Add the assigned labels to the Seurat object.

pbmc_labeled_obj = AddMetaData(pbmc3k.final, as.data.frame(singler_pred))

Compare the assigned labels to the original labels.

pbmc_labeled_obj@meta.data %>% select(labels, seurat_annotations) %>% table(useNA = "ifany")

                  seurat_annotations
labels             Naive CD4 T Memory CD4 T CD14+ Mono   B CD8 T FCGR3A+ Mono  NK  DC
  B_cell                     0            0          0 333     0            0   0   0
  CMP                        4            0          0   0     0            0   0   0
  Monocyte                   0            0        434   0     0          155   0  27
  NK_cell                    5            1          0   0    37            0 150   0
  Platelets                  0            0          0   0     0            0   0   0
  Pre-B_cell_CD34-           3            1         46   1     1            6   1   4
  Pro-B_cell_CD34+           0            0          0   4     0            0   0   0
  T_cells                  685          481          0   6   233            1   4   1
                  seurat_annotations
labels             Platelet
  B_cell                  0
  CMP                     0
  Monocyte                1
  NK_cell                 0
  Platelets              12
  Pre-B_cell_CD34-        0
  Pro-B_cell_CD34+        1
  T_cells                 0

SingleR also attempts to remove low quality or ambiguous assignments. Ambiguous assignments are based on the difference between the score for the assigned label and the median across all labels for each cell. Tuning parameters can be adjusted with pruneScores(). We can check how many cells are considered ambiguous and which populations they potentially belong to.

pbmc_labeled_obj@meta.data %>% select(pruned.labels, labels) %>% table(useNA = "ifany")

                  labels
pruned.labels      B_cell  CMP Monocyte NK_cell Platelets Pre-B_cell_CD34-
  B_cell              331    0        0       0         0                0
  CMP                   0    3        0       0         0                0
  Monocyte              0    0      614       0         0                0
  NK_cell               0    0        0     186         0                0
  Platelets             0    0        0       0        12                0
  Pre-B_cell_CD34-      0    0        0       0         0               63
  Pro-B_cell_CD34+      0    0        0       0         0                0
  T_cells               0    0        0       0         0                0
  <NA>                  2    1        3       7         0                0
                  labels
pruned.labels      Pro-B_cell_CD34+ T_cells
  B_cell                          0       0
  CMP                             0       0
  Monocyte                        0       0
  NK_cell                         0       0
  Platelets                       0       0
  Pre-B_cell_CD34-                0       0
  Pro-B_cell_CD34+                5       0
  T_cells                         0    1394
  <NA>                            0      17

Check the HPCA cell type labels overlaid onto the original UMAP visualization.

DimPlot(pbmc_labeled_obj, reduction = "umap", group.by = "labels") +
  scale_color_brewer(palette = "Set1")

HPCA subtypes

In addition to the 37 main cell types, the HPCA dataset also contains 157 subtypes. You can perform SingleR annotation for the subtypes.

singler_pred = SingleR(
  test = exp_common_mat,
  ref = singler_se,
  labels = singler_se$label.fine
)

Add the assigned labels to the Seurat object.

pbmc_labeled_obj = AddMetaData(pbmc3k.final, as.data.frame(singler_pred))

Check the assigned labels.

pbmc_labeled_obj@meta.data %>%
  select(labels) %>%
  table(useNA = "ifany")

.
                        B_cell                B_cell:immature 
                            26                            212 
                 B_cell:Memory                   B_cell:Naive 
                            17                             83 
            B_cell:Plasma_cell                            CMP 
                             2                              3 
                           MEP                       Monocyte 
                             1                              4 
                Monocyte:CD16-                 Monocyte:CD16+ 
                           381                            264 
       Monocyte:leukotriene_D4                        NK_cell 
                             1                            131 
          NK_cell:CD56hiCD62L+                    NK_cell:IL2 
                            23                             10 
                     Platelets               Pre-B_cell_CD34- 
                            11                             31 
                   T_cell:CD4+     T_cell:CD4+_central_memory 
                            33                            606 
   T_cell:CD4+_effector_memory              T_cell:CD4+_Naive 
                           353                            146 
                   T_cell:CD8+     T_cell:CD8+_Central_memory 
                           244                              5 
   T_cell:CD8+_effector_memory T_cell:CD8+_effector_memory_RA 
                            42                              1 
             T_cell:CD8+_naive             T_cell:gamma-delta 
                             5                              3 

Compare the assigned labels to the original labels. Subset to T-cells to keep the output more compact.

pbmc_labeled_obj@meta.data %>%
  select(labels, seurat_annotations) %>%
  filter(str_detect(seurat_annotations, "T")) %>%
  droplevels() %>%
  table(useNA = "ifany")

                                seurat_annotations
labels                           Naive CD4 T Memory CD4 T CD8 T
  CMP                                      3            0     0
  MEP                                      1            0     0
  NK_cell                                  1            0    22
  NK_cell:CD56hiCD62L+                     3            1     3
  Pre-B_cell_CD34-                         4            1     1
  T_cell:CD4+                             23            9     0
  T_cell:CD4+_central_memory             355          244     5
  T_cell:CD4+_effector_memory             78          192    83
  T_cell:CD4+_Naive                      134           12     0
  T_cell:CD8+                             51           18   151
  T_cell:CD8+_Central_memory               2            0     3
  T_cell:CD8+_effector_memory             38            4     0
  T_cell:CD8+_effector_memory_RA           0            0     1
  T_cell:CD8+_naive                        3            0     2
  T_cell:gamma-delta                       1            2     0

Compare the assigned labels to the original labels. Subset to monocytes to keep the output more compact.

pbmc_labeled_obj@meta.data %>%
  select(labels, seurat_annotations) %>%
  filter(str_detect(seurat_annotations, "Mono")) %>%
  droplevels() %>%
  table(useNA = "ifany")

                         seurat_annotations
labels                    CD14+ Mono FCGR3A+ Mono
  Monocyte                         3            0
  Monocyte:CD16-                 359            0
  Monocyte:CD16+                  95          162
  Monocyte:leukotriene_D4          1            0
  Pre-B_cell_CD34-                21            0
  T_cell:CD4+                      1            0

Annotation (using clustermole)

SingleR is able to label cells, but it requires a reference dataset. This is not a problem for PBMCs, but a perfect reference dataset may not always be available or you may have some unexpected populations.

A more exploratory and unbiased approach is possible with clustermole, an R package that provides a collection of cell type markers for thousands of human and mouse cell populations sourced from a variety of databases as well as methods to query them.

if (!require("clustermole")) {
   install.packages("clustermole")
}
library(clustermole)

Available cell types

We can retrieve a list of all markers in the ClusterMole database to see all the available options. Each row in the returned data frame is a combination of a single gene and its associated cell type.

markers_tbl = clustermole_markers()
head(markers_tbl)

Check the available cell types (ignore gene columns).

markers_tbl %>% distinct(celltype, organ, db)

Marker gene overlaps

We can find markers for the B-cell cluster Seurat’s FindMarkers function.

b_markers_df = FindMarkers(pbmc3k.final, ident.1 = "B", verbose = FALSE)

Check the markers table.

head(b_markers_df)

With the default cutoffs, this gives us a data frame with hundreds of genes. Let’s subset to just the top 20 genes.

b_markers = rownames(b_markers_df)
b_markers = head(b_markers, 20)
b_markers

 [1] "CD79A"     "MS4A1"     "CD79B"     "LINC00926" "TCL1A"     "HLA-DQA1" 
 [7] "VPREB3"    "HLA-DQB1"  "CD74"      "HLA-DRA"   "FCER2"     "HLA-DPB1" 
[13] "BANK1"     "HLA-DRB1"  "HLA-DPA1"  "TSPAN13"   "HLA-DQA2"  "FCRLA"    
[19] "CD37"      "HLA-DRB5" 

Check overlap of our B-cell markers with all cell type signatures.

overlaps_tbl = clustermole_overlaps(genes = b_markers, species = "hs")

Check the top scoring cell types for the B-cell cluster.

head(overlaps_tbl, 10)

Enrichment of markers

Calculate the average expression levels for the clusters for enrichment analysis using Seurat’s AverageExpression function, convert to a matrix, and log-transform.

avg_exp_mat = AverageExpression(pbmc3k.final, assays = "RNA", slot = "data")

Finished averaging RNA for cluster Naive CD4 T
Finished averaging RNA for cluster Memory CD4 T
Finished averaging RNA for cluster CD14+ Mono
Finished averaging RNA for cluster B
Finished averaging RNA for cluster CD8 T
Finished averaging RNA for cluster FCGR3A+ Mono
Finished averaging RNA for cluster NK
Finished averaging RNA for cluster DC
Finished averaging RNA for cluster Platelet

avg_exp_mat = as.matrix(avg_exp_mat$RNA)
avg_exp_mat = log1p(avg_exp_mat)

Check the average exression matrix.

head(avg_exp_mat)

              Naive CD4 T Memory CD4 T CD14+ Mono          B      CD8 T FCGR3A+ Mono
AL627309.1    0.006109960  0.005909767 0.04740195 0.00000000 0.02033764   0.00000000
AP006222.2    0.000000000  0.008172591 0.01082590 0.00000000 0.01184445   0.00000000
RP11-206L10.2 0.007425455  0.000000000 0.00000000 0.02043998 0.00000000   0.00000000
RP11-206L10.9 0.000000000  0.000000000 0.01044641 0.00000000 0.00000000   0.01192865
LINC00115     0.018938463  0.024390599 0.03685000 0.03814842 0.01929541   0.01458683
NOC2L         0.403772470  0.307346059 0.24101294 0.46155233 0.44279234   0.40564359
                     NK        DC Platelet
AL627309.1    0.0000000 0.0000000        0
AP006222.2    0.0000000 0.0000000        0
RP11-206L10.2 0.0000000 0.0812375        0
RP11-206L10.9 0.0000000 0.0000000        0
LINC00115     0.0569015 0.0000000        0
NOC2L         0.5311082 0.2841343        0

Run enrichment of all cell type signatures across all clusters.

enrichment_tbl = clustermole_enrichment(expr_mat = avg_exp_mat, species = "hs")

Selecting by score

Check the top scoring cell types for the B-cell cluster.

enrichment_tbl %>% filter(cluster == "B") %>% head(10)