Introduction

CellAnno is a automatic cell annotation tool with built-in human cell reference. To make one unified reference that can be used in any dataset, we compiled, curated and integrated single cell datasets from multiple sources including Human Protein Atlas, Tabula Sapiens, ArrayExpress and GEO and created a comprehensive human cell gene expression reference database named CellMap. It covers 200 cell types/states from 42 human tissue types, including 137 epithelial cell types, 22 stromal cell types and 19 major immune cell types. CellAnno uses the correlation-based algorithm SingleR to assign the most similar cell type from the reference. To make the process efficient and increase clarity of the results, CellAnno applies annotation at cluster level rather than single cell level. For immune cells, an option to run fine grained annotation (e.g. Treg, CD8+ effector T-cell) is available. With this option, a reference generated from ImmGen datasets (GSE15907,GSE37448) is used.

Example: a metastatic prostate cancer sample

The sample is a prostate cancer sample but collected from a metastatic site in liver. Therefore it is possible that captured cells orignate from both prostate and liver. The data can be downloaded from GEO.

library(Seurat)
library(ggplot2)
library(patchwork)
library(CellAnno)
crpc = read.csv("GSM6428953_1779_HMP05_IGO_10726_3_dense.csv.gz", row.names = 1)
crpc = CreateSeuratObject(counts = t(crpc[,-1]))
crpc = cellAnnotate(query.srt = crpc, ref.tissue = "prostate|liver")
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5591
## Number of edges: 222332
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9110
## Number of communities: 17
## Elapsed time: 0 seconds

The cellAnno function return a seurat object with cell type prediction added to the meta data. Now we can label cells with predicted cell types and color cells by prediction score. We indeed see cells from both prostate and liver.

p1 = DimPlot(crpc, group.by = "labels", label = TRUE, label.size = 3, repel = T)+NoLegend()
p2 = FeaturePlot(crpc, features = "score")
p1+p2

We can identify top marker genes for each annotated cell types to verify the accuracy of prediction.

Idents(crpc) = factor(crpc$labels, levels = unique(crpc$labels[order(as.numeric(crpc$cluster))]))
markers = FindAllMarkers(crpc, features = VariableFeatures(crpc), only.pos = TRUE)
topmarkers = getTopMarkers(mks = markers)
DotPlot(crpc, features = topmarkers)+RotatedAxis()+theme(axis.text.x = element_text(size = 7))

## Example: immune cells

The sample is 10k Bone Marrow Mononuclear Cells from 10x.

bmmc = Read10X_h5("10k_BMMNC_5pv2_nextgem_intron_10k_BMMNC_5pv2_nextgem_intron_count_sample_feature_bc_matrix.h5")
bmmc = CreateSeuratObject(counts = bmmc$`Gene Expression`)
bmmc = bmmc[, bmmc$nFeature_RNA>200]
bmmc = cellAnnotate(query.srt = bmmc, ref.tissue = "immune", level2 = TRUE)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 8275
## Number of edges: 320554
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9095
## Number of communities: 17
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 698
## Number of edges: 33896
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7257
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1350
## Number of edges: 52050
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7979
## Number of communities: 6
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 32
## Number of edges: 496
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.5000
## Number of communities: 1
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5255
## Number of edges: 203046
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8518
## Number of communities: 11
## Elapsed time: 0 seconds
p1 = DimPlot(bmmc, group.by = "level1", label = TRUE)+NoLegend()
p2 = DimPlot(bmmc, group.by = "level2", label = TRUE)+NoLegend()
p1+p2

bmmc$level1_level2 = paste(bmmc$level1, bmmc$level2, sep = ".")
bmmc$level1_level2[is.na(bmmc$level2)] = bmmc$level1[is.na(bmmc$level2)]
DimPlot(bmmc, group.by = "level1_level2", label = T)+NoLegend()

Idents(bmmc) = as.factor(bmmc$level1_level2)
markers = FindAllMarkers(bmmc, features = VariableFeatures(bmmc), only.pos = TRUE)
topmarkers = getTopMarkers(mks = markers)
DotPlot(bmmc, features = topmarkers)+RotatedAxis()+theme(axis.text.x = element_text(size = 7))

Predict cell-of-origin for bulk RNA-seq data

The same approach can be used to predict cell-of-origin (COO) for tumors using bulk RNA-seq data. Our comprehensive reference is extremely useful to predict COO for cancer with unknown primary. The testing data here is TCGA clear cell renal cell carcinoma cohort (KIRC). It is known that some samples in this cohort are chromophobe renal cell carcinoma (RCC). The two RCC subtypes originate from different lineages of nephron epithelial cells.

rcc = read.table("TCGA-KIRC.star_tpm.tsv.gz", sep = "\t", row.names = 1, header = T, check.names = FALSE)
rcc.genes = read.table("gencode.v36.annotation.gtf.gene.probemap", header = T)
dup.genes = unique(rcc.genes$gene[duplicated(rcc.genes$gene)])
rcc.genes = rcc.genes[!(rcc.genes$gene %in% dup.genes),]
rcc = rcc[rcc.genes$id, ]
rownames(rcc) = rcc.genes$gene
rcc.coo = cooPredict(logtpm.matrix = rcc, ref = "kidney")
sample.type = substr(colnames(rcc), 14, 15)
table(rcc.coo$labels, sample.type) # 01 is tumor, 11 is normal, 05 is new primary
##                       sample.type
##                         01  05  11
##   collecting_duct        7   0  35
##   intercalated_a        19   0   0
##   proximal_tubule       78   0  33
##   proximal_tubule_b    382   0   2
##   proximal_tubule_c     45   0   0
##   thick_ascending_limb   0   0   2
##   thin_ascending_limb    6   1   0