Introduction
This vignette demonstrates the powerful and intuitive workflow for single-cell RNA-seq annotation provided by the easybio package. The process is designed to combine the speed of automated database matching with the reliability of interactive verification and manual curation.
The core workflow follows three logical steps:
- Automated Annotation: Use
match_ref()to quickly get a list of potential cell types for each cluster based on its marker genes. - Verification & Exploration: Interactively investigate the automated results using
check_marker()andplot_seurat_dot()to build confidence in the annotations. This step helps answer two critical questions:- “Why was this annotation made?” (Which of my genes matched the database?)
- “Is this annotation correct?” (Are the canonical markers for this cell type expressed in my cluster?)
- Final Curation: Based on the evidence gathered, use
finsert()to assign the final, high-confidence cell type labels.
You can also view the R script for this workflow by running:
file.show(system.file(package = 'easybio', 'example-single-cell.R'))
Setup
First, let’s load the necessary libraries and the example marker data included with easybio. This data is derived from the 10x Genomics 3k PBMC dataset.
litedown::reactor(warning = FALSE) # vignette setting
library(easybio)
library(Seurat)
library(data.table)
# The pbmc.markers dataset is included in easybio
head(pbmc.markers)
| p_val | avg_log2FC | pct.1 | pct.2 | p_val_adj | cluster | gene | |
|---|---|---|---|---|---|---|---|
| RPS12 | 0 | 0.739 | 1.000 | 0.991 | 0 | 0 | RPS12 |
| RPS6 | 0 | 0.693 | 1.000 | 0.995 | 0 | 0 | RPS6 |
| RPS27 | 0 | 0.737 | 0.999 | 0.992 | 0 | 0 | RPS27 |
| RPL32 | 0 | 0.627 | 0.999 | 0.995 | 0 | 0 | RPL32 |
| RPS14 | 0 | 0.634 | 1.000 | 0.994 | 0 | 0 | RPS14 |
| RPS25 | 0 | 0.769 | 0.997 | 0.975 | 0 | 0 | RPS25 |
Step 1: Automated Annotation with match_ref
We begin by feeding the cluster markers (from Seurat::FindAllMarkers) into match_ref(). This function compares our markers against the CellMarker 3.0 database and returns a ranked list of potential cell types for each cluster.
Bounding the search to the right tissue
By default match_ref() searches the whole database, which mixes every organ: for PBMC data that puts hepatocytes and alveolar cells in the same candidate list as monocytes. Each database entry carries two labels that can narrow it:
tissue_class, the organ or system the entry came from (Blood,Bone marrow,Lung, …), andtissue_type, the sample the original study described (Peripheral blood,Peripheral blood mononuclear cell,Whole blood, …).
They are two labels on an entry rather than a hierarchy — one tissue_type can be listed under several classes — and match_ref() combines them with AND. Bounding tissue_class is the useful move here; bounding tissue_type as well removes evidence instead of correcting for anything: restricted to "Peripheral blood" alone, cluster 8’s platelet markers keep only a fraction of their support and it comes out as a dendritic cell instead. Our sample is peripheral blood mononuclear cells, so blood and bone marrow, where the monocytes and dendritic cells circulating in it are made, is the natural bound:
# the organ-level classes, and what one of them actually holds
length(available_tissue_class("Human"))
#> [1] 87
available_tissue_type("Human", tissue_class = "Blood")
#> [1] "Peripheral blood" "Peripheral blood mononuclear cell"
#> [3] "Blood" "Umbilical cord blood"
#> [5] "Lymph" "Serum"
#> [7] "Artery" "Whole blood"
#> [9] "Plasma" "Thymus"
marker_matched <- match_ref(
marker = pbmc.markers, n = 50, spc = "Human",
tissue_class = c("Blood", "Bone marrow")
)
# Let's look at the top 2 potential cell types for each cluster
marker_matched[, head(.SD, 2), by = cluster]
| cluster | cell_name | uniqueN | N | ordered_symbol | orderN | pct_with |
|---|---|---|---|---|---|---|
| 0 | T cell | 30 | 80 | CCR7,LEF1,CD5,MAL,LRRN3,NELL2,… | 9,8,6,5,4,4,… | 0.447,0.342,0.127,0.268,0.038,0.158,… |
| 0 | CD4+ T cell | 21 | 81 | MAL,TSHZ2,LEF1,CCR7,TRABD2A,FHIT,… | 8,8,7,6,6,5,… | 0.268,0.094,0.342,0.447,0.199,0.200,… |
| 1 | Monocyte | 31 | 178 | CD14,S100A9,S100A8,CREB5,S100A12,NRG1,… | 51,13,12, 9, 9, 8,… | 0.667,0.996,0.975,0.031,0.277,0.146,… |
| 1 | Myeloid cell | 16 | 43 | CD14,S100A8,S100A9,CD93,S100A12,CLEC4E,… | 6,6,6,4,4,3,… | 0.667,0.975,0.996,0.077,0.277,0.150,… |
| 2 | T cell | 21 | 50 | GPR171,CD2,TNFRSF25,CD40LG,AQP3,CCR6,… | 9,6,6,5,2,2,… | 0.132,0.651,0.162,0.263,0.420,0.063,… |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 6 | T cell | 9 | 27 | PRF1,FGFBP2,SPON2,GNLY,GZMB,KLRF1,… | 6,4,4,3,3,3,… | 0.948,0.877,0.729,0.961,0.961,0.271,… |
| 7 | Dendritic cell | 24 | 102 | FCER1A,CLEC4C,CLEC10A,ENHO,FLT3,LILRA4,… | 15, 9, 8, 7, 7, 7,… | 0.812,0.188,0.688,0.438,0.125,0.156,… |
| 7 | Conventional dendritic cell | 21 | 101 | FCER1A,CD1E,ENHO,PKIB,CLEC10A,CLIC2,… | 11,10, 9, 9, 8, 8,… | 0.812,0.062,0.438,0.188,0.688,0.281,… |
| 8 | Platelet | 34 | 111 | PPBP,PF4,GP9,TUBB1,ITGA2B,SPARC,… | 13, 9, 6, 5, 4, 4,… | 1.000,1.000,0.923,0.846,0.846,1.000,… |
| 8 | Conventional dendritic cell | 15 | 43 | GP1BA,TMEM40,AP001189.4,CLU,GNG11,GP9,… | 4,4,3,3,3,3,… | 0.615,0.769,0.846,0.846,1.000,0.923,… |
The output table gives us uniqueN (the number of unique matching markers) and N (the total number of matches), which helps rank the potential annotations.
We can create a quick preliminary annotation by taking the top hit for each cluster.
cl2cell_auto <- marker_matched[, head(.SD, 1), by = .(cluster)]
cl2cell_auto <- setNames(cl2cell_auto[["cell_name"]], cl2cell_auto[["cluster"]])
print("Initial automated annotation:")
#> [1] "Initial automated annotation:"
cl2cell_auto
#> 0 1 2
#> "T cell" "Monocyte" "T cell"
#> 3 4 5
#> "B cell" "T cell" "Monocyte"
#> 6 7 8
#> "Natural killer cell" "Dendritic cell" "Platelet"
We can also get a global view of all possible annotations using plot_possible_cell.
plot_possible_cell(marker_matched[, head(.SD), by = .(cluster)], min_unique_n = 2)
The same view can be coloured by how well the evidence is actually detected. With value = "pct" the fill shows the share of each candidate’s matched markers whose detection rate reaches min_pct, which is the quickest way to spot annotations that rest on barely expressed genes.
plot_possible_cell(marker_matched[, head(.SD), by = .(cluster)], min_unique_n = 2, value = "pct")
Step 2: Verification and Exploration
This is the most critical step. Instead of blindly trusting the automated result, we use easybio’s tools to verify it.
Answering “Why was this annotation made?”
To see the evidence behind an annotation, we use check_marker() with cis = TRUE. This shows us which of our own marker genes from our data matched the database for a given annotation.
# Let's investigate clusters 1, 5, and 7
local_evidence <- check_marker(marker_matched, cl = c(1, 5, 7), top_cell_n = 2, cis = TRUE)
print(local_evidence)
#> $Monocyte
#> [1] "CD14" "S100A9" "S100A8" "CREB5" "S100A12" "NRG1"
#> [7] "ASGR1" "FCAR" "MS4A6A" "ALDH1A1" "CLEC4E" "LINC00937"
#> [13] "MITF" "CCDC149" "CD93" "FCGR1A" "FOLR3" "HGF"
#> [19] "MGST1" "CYP27A1" "DYSF" "LGALS2" "LIN7A" "MPO"
#> [25] "QPCT" "RNASE2" "TMEM150B" "CCL2" "ELANE" "PLA2G7"
#> [31] "SIRPA"
#>
#> $`Myeloid cell`
#> [1] "CD14" "S100A8" "S100A9" "CD93" "S100A12" "CLEC4E" "MS4A6A"
#> [8] "LGALS2" "MPO" "ALDH1A1" "CCDC149" "CYP27A1" "ELANE" "FCGR1A"
#> [15] "MGST1" "PLA2G7"
#>
#> $Monocyte
#> [1] "MS4A7" "FCGR3B" "CDKN1C" "LYPD2" "C1QA" "CASP5"
#> [7] "SFTPD" "C1QB" "GPR20" "HES4" "SMPDL3A" "ARHGAP29"
#> [13] "CKB" "CYP4F22" "PPP1R17" "VMO1"
#>
#> $Macrophage
#> [1] "C1QA" "C1QB" "C3" "FCGR3B" "KCNMA1" "MS4A7" "MYO10"
#> [8] "SMPDL3A"
#>
#> $`Dendritic cell`
#> [1] "FCER1A" "CLEC4C" "CLEC10A" "ENHO" "FLT3" "LILRA4"
#> [7] "TNFRSF21" "CLIC2" "IL1R2" "P2RY6" "CD1E" "GAS6"
#> [13] "SERPINF1" "AXL" "DNASE1L3" "LAMP5" "LRRC26" "PROC"
#> [19] "SCT" "TIFAB" "CRYM" "PMP22" "SEZ6L" "SMPD3"
#>
#> $`Conventional dendritic cell`
#> [1] "FCER1A" "CD1E" "ENHO" "PKIB" "CLEC10A" "CLIC2"
#> [7] "IL1R2" "FBLN2" "FLT3" "IDO1" "P2RY6" "AXL"
#> [13] "DNASE1L3" "SERPINF2" "GAS6" "GEM" "GFRA2" "PLS3"
#> [19] "PMP22" "SDS" "SEZ6L"
#>
Answering “Is this annotation correct?”
To validate an annotation, we use check_marker() with cis = FALSE (the default). This fetches the canonical markers for the suggested cell type from the database. We can then check if these well-known markers are expressed in our cluster.
canonical_markers <- check_marker(marker_matched, cl = c(1, 5, 7), top_cell_n = 2, cis = FALSE)
print(canonical_markers)
#> $Monocyte
#> [1] "CD14" "FCN1" "LYZ" "SLC11A1" "VCAN" "CD300E"
#> [7] "IRAK3" "CLEC12A" "SERPINA1" "FCGR3A"
#>
#> $`Myeloid cell`
#> [1] "LYZ" "CD68" "FCN1" "CD14" "S100A9" "S100A8" "SULF2" "CD86"
#> [9] "KLF4" "FPR1"
#>
#> $Macrophage
#> [1] "CD68" "CD163" "C1QA" "C1QB" "CD14" "C1QC" "MRC1"
#> [8] "SLCO2B1" "CXCL2" "PHLDA1"
#>
#> $`Dendritic cell`
#> [1] "FCER1A" "CD1C" "CLEC4C" "CLEC10A" "PLD4" "CST3" "NEGR1"
#> [8] "ENHO" "FLT3" "SLC41A2"
#>
#> $`Conventional dendritic cell`
#> [1] "CD1C" "FCER1A" "NDRG2" "CD1E" "PKIB" "ENHO" "PDLIM1" "CYP2S1"
#> [9] "CCSER1" "COL9A2"
#>
Visual Confirmation with plot_seurat_dot
The best way to check marker expression is visually. plot_seurat_dot is designed to work seamlessly with check_marker.
The entire pipeline from annotation to visualization can be done in a single, elegant pipe:
# For this example to be runnable, we need a Seurat object.
# We'll create a minimal one. In your real workflow, you would use your own srt object.
marker_genes <- unique(pbmc.markers$gene)
counts <- matrix(
abs(rnorm(length(marker_genes) * 50, mean = 1, sd = 2)),
nrow = length(marker_genes),
ncol = 50
)
rownames(counts) <- marker_genes
colnames(counts) <- paste0("cell_", 1:50)
srt <- CreateSeuratObject(counts = counts)
# Assign clusters that match the pbmc.markers data
srt$seurat_clusters <- sample(0:8, 50, replace = TRUE)
Idents(srt) <- "seurat_clusters"
# Now, let's plot the evidence for clusters 1, 5, and 7
match_ref(marker = pbmc.markers, n = 50, spc = "Human") |>
check_marker(cl = c(1, 5, 7), top_cell_n = 2, cis = TRUE) |>
plot_seurat_dot(srt = srt)
This dot plot clearly shows the expression of the genes that led to the annotations for clusters 1, 5, and 7, allowing us to confidently assess the results.
Step 3: Final Manual Curation
After reviewing the evidence from the dot plots, we can make our final, informed decision. The finsert function provides a convenient way to create the final annotation vector.
# Based on our exploration, we finalize the annotations
cl2cell_final <- finsert(
list(
c(3) ~ "B cell",
c(8) ~ "Megakaryocyte",
c(7) ~ "DC",
c(1, 5) ~ "Monocyte",
c(0, 2, 4) ~ "Naive CD8+ T cell",
c(6) ~ "Natural killer cell"
),
len = 9 # Ensure vector length covers all clusters (0-8)
)
print("Final curated annotation:")
#> [1] "Final curated annotation:"
cl2cell_final
#> 0 1 2
#> "Naive CD8+ T cell" "Monocyte" "Naive CD8+ T cell"
#> 3 4 5
#> "B cell" "Naive CD8+ T cell" "Monocyte"
#> 6 7 8
#> "Natural killer cell" "DC" "Megakaryocyte"
This cl2cell_final vector can now be added to your Seurat object’s metadata for downstream analysis and plotting.
Using a Custom Marker Database
For specialized analyses, such as focusing on a specific tissue, working with a non-model organism, or using a proprietary list of markers, you can provide your own custom reference to match_ref.
The reference must be a data.frame (or data.table) with at least two columns: cell_name and marker. The easiest way to create this is from a named list.
Step 1: Create a named list of your custom markers.
custom_ref_list <- list(
"T-cell" = c("CD3D", "CD3E", "CD3G"),
"B-cell" = c("CD79A", "MS4A1"),
"Myeloid" = c("LYZ", "CST3", "AIF1")
)
print(custom_ref_list)
#> $`T-cell`
#> [1] "CD3D" "CD3E" "CD3G"
#>
#> $`B-cell`
#> [1] "CD79A" "MS4A1"
#>
#> $Myeloid
#> [1] "LYZ" "CST3" "AIF1"
#>
Step 2: Convert the list to the required data.frame format.
easybio provides the list2dt helper function for this.
custom_ref_df <- list2dt(custom_ref_list, col_names = c("cell_name", "marker"))
head(custom_ref_df)
| cell_name | marker |
|---|---|
| T-cell | CD3D |
| T-cell | CD3E |
| T-cell | CD3G |
| B-cell | CD79A |
| B-cell | MS4A1 |
| Myeloid | LYZ |
Step 3: Run match_ref with the ref parameter.
When ref is provided, the function ignores the spc, tissue_class, and tissue_type parameters for matching.
marker_custom <- match_ref(
marker = pbmc.markers,
n = 50,
ref = custom_ref_df
)
# Note that the cell_name column now contains our custom cell types
marker_custom[, head(.SD, 2), by = cluster]
| cluster | cell_name | uniqueN | N | ordered_symbol | orderN | pct_with |
|---|---|---|---|---|---|---|
| 3 | B-cell | 2 | 2 | CD79A,MS4A1 | 1,1 | 0.936,0.855 |
Additional Utilities
easybio also provides functions for direct queries.
get_marker()
Directly retrieve markers for any cell type of interest.
get_marker(spc = "Human", cell = c("Monocyte", "Neutrophil"), number = 5, min_count = 1)
#> $Monocyte
#> [1] "FCN1" "CD14" "S100A8" "S100A9" "LILRA5"
#>
#> $Neutrophil
#> [1] "S100A9" "S100A8" "FCGR3B" "CSF3R" "CXCL8"
#>
plot_marker_distribution()
Check the distribution of a specific marker across all cell types and tissues in the database.
plot_marker_distribution(mkr = "CD68")
