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 contains hematopoietic and stromal bone marrow populations (Baccin et al.). The version used here is a subset of the original dataset to have more similar population sizes and speed up processing.

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.11’

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.

library(Seurat)
library(ggplot2)
library(cowplot)
library(ggsci)
library(dplyr)
library(stringr)

Load the dataset.

so = readRDS(url("https://osf.io/cvnqb/download"))
so

An object of class Seurat 
16701 features across 2821 samples within 1 assay 
Active assay: RNA (16701 features, 2872 variable features)
 3 dimensional reductions calculated: pca, tsne, umap

Check the experiment labels overlaid onto the UMAP visualization.

DimPlot(so, reduction = "umap", group.by = "experiment", cells = sample(colnames(so))) +
  scale_color_nejm()

Check the experiment labels overlaid onto the tSNE visualization, as shown in the original publication (original figure).

DimPlot(so, reduction = "tsne", group.by = "experiment", cells = sample(colnames(so))) +
  scale_color_nejm()

Check the cell type labels overlaid onto the tSNE visualization.

DimPlot(so, reduction = "tsne", group.by = "celltype", cells = sample(colnames(so))) +
  scale_color_igv()

Annotation using SingleR

Load SingleR (install with BiocManager::install("SingleR")).

library(SingleR)

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

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

[1] 16701  2821

MouseRNAseqData 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 will try the mouse dataset from 358 bulk RNA-seq samples of sorted cell populations as the reference. It provides normalized expression values for samples that have been assigned to one of 18 main cell types and 28 subtypes.

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: 21214 358 
metadata(0):
assays(1): logcounts
rownames(21214): Xkr4 Rp1 ... LOC100039574 LOC100039753
rowData names(0):
colnames(358): ERR525589Aligned ERR525592Aligned ... SRR1044043Aligned
  SRR1044044Aligned
colData names(3): label.main label.fine label.ont

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] 13929

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))

Add the assigned labels to the Seurat object.

so_labeled = AddMetaData(so, as.data.frame(singler_pred))

Check the assigned labels to the original labels.

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

.
       Adipocytes           B cells    Cardiomyocytes   Dendritic cells Endothelial cells 
                6               400                 3                 1               201 
     Erythrocytes       Fibroblasts      Granulocytes       Macrophages         Monocytes 
              372               815               191                 2               534 
         NK cells  Oligodendrocytes           T cells 
               72                90               134 

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.

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

.
       Adipocytes           B cells    Cardiomyocytes   Dendritic cells Endothelial cells 
                6               400                 3                 1               197 
     Erythrocytes       Fibroblasts      Granulocytes       Macrophages         Monocytes 
              372               815               191                 2               534 
         NK cells  Oligodendrocytes           T cells              <NA> 
               66                88               115                31 

Visualize MouseRNAseqData cell type labels.

DimPlot(so_labeled, reduction = "tsne", group.by = "labels") +
  scale_color_igv()

Annotation using clustermole

SingleR is able to label cells, but it requires a reference dataset.

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.

Load 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 each cell type and its genes.

markers_tbl = clustermole_markers()
head(markers_tbl)

Check the available cell types (ignore gene columns).

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

Marker gene overlaps

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

Idents(so) = "celltype"
avg_exp_mat = AverageExpression(so, assays = "RNA", slot = "data")

Finished averaging RNA for cluster Adipo-CAR
Finished averaging RNA for cluster Arteriolar-ECs
Finished averaging RNA for cluster Arteriolar-fibro
Finished averaging RNA for cluster B-cell
Finished averaging RNA for cluster Chondrocytes
Finished averaging RNA for cluster Dendritic-cells
Finished averaging RNA for cluster Endosteal-fibro
Finished averaging RNA for cluster Eo-Baso-prog
Finished averaging RNA for cluster Ery-Mk-prog
Finished averaging RNA for cluster Ery-prog
Finished averaging RNA for cluster Erythroblasts
Finished averaging RNA for cluster Fibro-Chondro-p
Finished averaging RNA for cluster Gran-Mono-prog
Finished averaging RNA for cluster large-pre-B
Finished averaging RNA for cluster LMPPs
Finished averaging RNA for cluster Mk-prog
Finished averaging RNA for cluster Mono-prog
Finished averaging RNA for cluster Monocytes
Finished averaging RNA for cluster Myofibroblasts
Finished averaging RNA for cluster Neutro-prog
Finished averaging RNA for cluster Neutrophils
Finished averaging RNA for cluster Ng2-MSCs
Finished averaging RNA for cluster NK-cells
Finished averaging RNA for cluster Osteo-CAR
Finished averaging RNA for cluster Osteoblasts
Finished averaging RNA for cluster pro-B
Finished averaging RNA for cluster Schwann-cells
Finished averaging RNA for cluster Sinusoidal-ECs
Finished averaging RNA for cluster small-pre-B
Finished averaging RNA for cluster Smooth-muscle
Finished averaging RNA for cluster Stromal-fibro
Finished averaging RNA for cluster T-cells

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

Check the average expression matrix.

avg_exp_mat[1:5, 1:5]

        Adipo-CAR Arteriolar-ECs Arteriolar-fibro    B-cell Chondrocytes
Sox17   0.0000000      2.5800606        0.0000000 0.0000000   0.00000000
Mrpl15  0.4315376      0.5284282        0.3304569 0.8481054   0.07307574
Lypla1  0.1990537      0.4477973        0.2448583 0.6352917   0.09397463
Gm37988 0.0000000      0.0000000        0.0000000 0.0000000   0.00000000
Tcea1   0.5620502      0.7077588        0.6135480 0.7798060   0.52126901

Check the cell type names.

levels(Idents(so))

 [1] "Adipo-CAR"        "Arteriolar-ECs"   "Arteriolar-fibro" "B-cell"          
 [5] "Chondrocytes"     "Dendritic-cells"  "Endosteal-fibro"  "Eo-Baso-prog"    
 [9] "Ery-Mk-prog"      "Ery-prog"         "Erythroblasts"    "Fibro-Chondro-p" 
[13] "Gran-Mono-prog"   "large-pre-B"      "LMPPs"            "Mk-prog"         
[17] "Mono-prog"        "Monocytes"        "Myofibroblasts"   "Neutro-prog"     
[21] "Neutrophils"      "Ng2-MSCs"         "NK-cells"         "Osteo-CAR"       
[25] "Osteoblasts"      "pro-B"            "Schwann-cells"    "Sinusoidal-ECs"  
[29] "small-pre-B"      "Smooth-muscle"    "Stromal-fibro"    "T-cells"         

Find markers for the B-cell cluster.

b_genes = rownames(avg_exp_mat[avg_exp_mat[, "B-cell"] == rowMaxs(avg_exp_mat), ])
length(b_genes)

[1] 517

b_markers_df = FindMarkers(so, ident.1 = "B-cell", features = b_genes, verbose = FALSE)

nrow(b_markers_df)

[1] 165

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] "Ms4a1"         "Fcmr"          "Cd74"          "Ly6d"          "Gm43603"      
 [6] "Bank1"         "2010309G21Rik" "Fcer2a"        "H2-DMb2"       "H2-Eb1"       
[11] "Cd79a"         "H2-Aa"         "Tnfrsf13c"     "Ltb"           "Cd79b"        
[16] "Ccr7"          "Fcrl1"         "Siglecg"       "Cd83"          "Srpk3"        

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

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

Registered S3 method overwritten by 'cli':
  method     from    
  print.boxx spatstat

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

head(overlaps_tbl, 10)

Find markers for the Adipo-CAR cluster.

acar_genes = rownames(avg_exp_mat[avg_exp_mat[, "Adipo-CAR"] == rowMaxs(avg_exp_mat), ])
length(acar_genes)

[1] 888

acar_markers_df = FindMarkers(so, ident.1 = "Adipo-CAR", features = acar_genes, verbose = FALSE)

nrow(acar_markers_df)

[1] 326

Check the markers table.

head(acar_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.

acar_markers = rownames(acar_markers_df)
acar_markers = head(acar_markers, 20)
acar_markers

 [1] "Adipoq"        "Kng1"          "Kng2"          "Esm1"          "Cxcl12"       
 [6] "Lpl"           "Gdpd2"         "Agt"           "Dpep1"         "Lepr"         
[11] "Fst"           "Chrdl1"        "Pdzrn4"        "Kitl"          "Ccl19"        
[16] "Ptx3"          "Ackr4"         "1500009L16Rik" "Gas6"          "Serpina12"    

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

overlaps_tbl = clustermole_overlaps(genes = acar_markers, species = "mm")

Check the top scoring cell types for the Adipo-CAR cluster.

head(overlaps_tbl, 10)

Find markers for the Osteoblasts cluster.

o_genes = rownames(avg_exp_mat[avg_exp_mat[, "Osteoblasts"] == rowMaxs(avg_exp_mat), ])
length(o_genes)

[1] 390

o_markers_df = FindMarkers(so, ident.1 = "Osteoblasts", features = o_genes, verbose = FALSE)

nrow(o_markers_df)

[1] 182

Check the markers table.

head(o_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.

o_markers = rownames(o_markers_df)
o_markers = head(o_markers, 20)
o_markers

 [1] "Cpz"           "Smpd3"         "Col22a1"       "Ifitm5"        "Mlip"         
 [6] "Bglap"         "Lipc"          "Cgref1"        "Col13a1"       "Entpd3"       
[11] "Fabp3"         "Bglap2"        "Cthrc1"        "Bglap3"        "Rerg"         
[16] "Cdo1"          "Car3"          "RP23-457J22.1" "Col24a1"       "Bmp3"         

Check overlap of our Osteoblasts markers with all cell type signatures.

overlaps_tbl = clustermole_overlaps(genes = o_markers, species = "mm")

Check the top scoring cell types for the Osteoblasts cluster.

head(overlaps_tbl, 10)

Enrichment of markers

Run enrichment of all cell type signatures across all clusters.

enrich_tbl = clustermole_enrichment(expr_mat = avg_exp_mat, species = "mm")

Some gene sets have size one. Consider setting 'min.sz > 1'.Selecting by score

Most enriched cell types for the B-cell cluster.

enrich_tbl %>% filter(cluster == "B-cell") %>% select(-cluster) %>% head(10)

Most enriched cell types for the Adipo-CAR cluster.

enrich_tbl %>% filter(cluster == "Adipo-CAR") %>% select(-cluster) %>% head(10)

Most enriched cell types for the Osteoblasts cluster.

enrich_tbl %>% filter(cluster == "Osteoblasts") %>% select(-cluster) %>% head(10)

---

previous tutorials