Mapping and annotating query datasets

Introduction to single-cell reference mapping

In this vignette, we first build an integrated reference and then demonstrate how to leverage this reference to annotate new query datasets. Generating an integrated reference follows the same workflow described in more detail in the integration introduction vignette. Once generated, this reference can be used to analyze additional query datasets through tasks like cell type label transfer and projecting query cells onto reference UMAPs. Notably, this does not require correction of the underlying raw query data and can therefore be an efficient strategy if a high quality reference is available.

Dataset preprocessing

For the purposes of this example, we’ve chosen human pancreatic islet cell datasets produced across four technologies, CelSeq (GSE81076) CelSeq2 (GSE85241), Fluidigm C1 (GSE86469), and SMART-Seq2 (E-MTAB-5061). For convenience, we distribute this dataset through our SeuratData package. The metadata contains the technology (tech column) and cell type annotations (celltype column) for each cell in the four datasets.

library(Seurat)
library(SeuratData)

InstallData("panc8")

To construct a reference, we will identify ‘anchors’ between the individual datasets. First, we split the combined object into a list, with each dataset as an element (this is only necessary because the data was bundled together for easy distribution).

data("panc8")
pancreas.list <- SplitObject(panc8, split.by = "tech")
pancreas.list <- pancreas.list[c("celseq", "celseq2", "fluidigmc1", "smartseq2")]

Prior to finding anchors, we perform standard preprocessing (log-normalization), and identify variable features individually for each. Note that Seurat implements an improved method for variable feature selection based on a variance stabilizing transformation ("vst")

for (i in 1:length(pancreas.list)) {
    pancreas.list[[i]] <- NormalizeData(pancreas.list[[i]], verbose = FALSE)
    pancreas.list[[i]] <- FindVariableFeatures(pancreas.list[[i]], selection.method = "vst", nfeatures = 2000,
        verbose = FALSE)
}

Integration of 3 pancreatic islet cell datasets

Next, we identify anchors using the FindIntegrationAnchors() function, which takes a list of Seurat objects as input. Here, we integrate three of the objects into a reference (we will use the fourth later in this vignette as a query dataset to demonstrate mapping).

• We use all default parameters here for identifying anchors, including the ‘dimensionality’ of the dataset (30; feel free to try varying this parameter over a broad range, for example between 10 and 50).

reference.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")]
pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

We then pass these anchors to the IntegrateData() function, which returns a Seurat object.

• The returned object will contain a new Assay, which holds an integrated (or ‘batch-corrected’) expression matrix for all cells, enabling them to be jointly analyzed.

pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)

After running IntegrateData(), the Seurat object will contain a new Assay with the integrated expression matrix. Note that the original (uncorrected values) are still stored in the object in the “RNA” assay, so you can switch back and forth.

We can then use this new integrated matrix for downstream analysis and visualization. Here we scale the integrated data, run PCA, and visualize the results with UMAP. The integrated datasets cluster by cell type, instead of by technology.

library(ggplot2)
library(cowplot)
library(patchwork)
# switch to integrated assay. The variable features of this assay are automatically set during
# IntegrateData
DefaultAssay(pancreas.integrated) <- "integrated"
# Run the standard workflow for visualization and clustering
pancreas.integrated <- ScaleData(pancreas.integrated, verbose = FALSE)
pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 30, verbose = FALSE)
pancreas.integrated <- RunUMAP(pancreas.integrated, reduction = "pca", dims = 1:30, verbose = FALSE)
p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "tech")
p2 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype", label = TRUE, repel = TRUE) +
    NoLegend()
p1 + p2

Cell type classification using an integrated reference

Seurat also supports the projection of reference data (or meta data) onto a query object. While many of the methods are conserved (both procedures begin by identifying anchors), there are two important distinctions between data transfer and integration:

1. In data transfer, Seurat does not correct or modify the query expression data.
2. In data transfer, Seurat has an option (set by default) to project the PCA structure of a reference onto the query, instead of learning a joint structure with CCA. We generally suggest using this option when projecting data between scRNA-seq datasets.

After finding anchors, we use the TransferData() function to classify the query cells based on reference data (a vector of reference cell type labels). TransferData() returns a matrix with predicted IDs and prediction scores, which we can add to the query metadata.

pancreas.query <- pancreas.list[["fluidigmc1"]]
pancreas.anchors <- FindTransferAnchors(reference = pancreas.integrated, query = pancreas.query,
    dims = 1:30, reference.reduction = "pca")
predictions <- TransferData(anchorset = pancreas.anchors, refdata = pancreas.integrated$celltype,
    dims = 1:30)
pancreas.query <- AddMetaData(pancreas.query, metadata = predictions)

Because we have the original label annotations from our full integrated analysis, we can evaluate how well our predicted cell type annotations match the full reference. In this example, we find that there is a high agreement in cell type classification, with over 96% of cells being labeled correctly.

pancreas.query$prediction.match <- pancreas.query$predicted.id == pancreas.query$celltype
table(pancreas.query$prediction.match)

## 
## FALSE  TRUE 
##    21   617

To verify this further, we can examine some canonical cell type markers for specific pancreatic islet cell populations. Note that even though some of these cell types are only represented by one or two cells (e.g. epsilon cells), we are still able to classify them correctly.

table(pancreas.query$predicted.id)

## 
##             acinar activated_stellate              alpha               beta 
##                 22                 17                253                256 
##              delta             ductal        endothelial              gamma 
##                 22                 30                 12                 18 
##         macrophage               mast            schwann 
##                  1                  2                  5

VlnPlot(pancreas.query, c("REG1A", "PPY", "SST", "GHRL", "VWF", "SOX10"), group.by = "predicted.id")

Unimodal UMAP Projection

In Seurat v4, we also enable projection of a query onto the reference UMAP structure. This can be achieved by computing the reference UMAP model and then calling MapQuery() instead of TransferData().

pancreas.integrated <- RunUMAP(pancreas.integrated, dims = 1:30, reduction = "pca", return.model = TRUE)
pancreas.query <- MapQuery(anchorset = pancreas.anchors, reference = pancreas.integrated, query = pancreas.query,
    refdata = list(celltype = "celltype"), reference.reduction = "pca", reduction.model = "umap")

pancreas.query <- TransferData(anchorset = pancreas.anchors, reference = pancreas.integrated, query = pancreas.query,
    refdata = list(celltype = "celltype"))
pancreas.query <- IntegrateEmbeddings(anchorset = pancreas.anchors, reference = pancreas.integrated,
    query = pancreas.query, new.reduction.name = "ref.pca")
pancreas.query <- ProjectUMAP(query = pancreas.query, query.reduction = "ref.pca", reference = pancreas.integrated,
    reference.reduction = "pca", reduction.model = "umap")

We can now visualize the query cells alongside our reference.

p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype", label = TRUE, label.size = 3,
    repel = TRUE) + NoLegend() + ggtitle("Reference annotations")
p2 <- DimPlot(pancreas.query, reduction = "ref.umap", group.by = "predicted.celltype", label = TRUE,
    label.size = 3, repel = TRUE) + NoLegend() + ggtitle("Query transferred labels")
p1 + p2

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.1.2               cowplot_1.1.1                
##  [3] ggplot2_3.4.1                 thp1.eccite.SeuratData_3.1.5 
##  [5] stxBrain.SeuratData_0.1.1     ssHippo.SeuratData_3.1.4     
##  [7] pbmcsca.SeuratData_3.0.0      pbmcMultiome.SeuratData_0.1.2
##  [9] pbmc3k.SeuratData_3.1.4       panc8.SeuratData_3.0.2       
## [11] ifnb.SeuratData_3.1.0         hcabm40k.SeuratData_3.0.0    
## [13] bmcite.SeuratData_0.3.0       SeuratData_0.2.2             
## [15] SeuratObject_4.1.3            Seurat_4.3.0                 
## 
## loaded via a namespace (and not attached):
##   [1] systemfonts_1.0.4      plyr_1.8.8             igraph_1.4.1          
##   [4] lazyeval_0.2.2         sp_1.6-0               splines_4.2.0         
##   [7] listenv_0.9.0          scattermore_0.8        digest_0.6.31         
##  [10] htmltools_0.5.4        fansi_1.0.4            magrittr_2.0.3        
##  [13] memoise_2.0.1          tensor_1.5             cluster_2.1.3         
##  [16] ROCR_1.0-11            globals_0.16.2         matrixStats_0.63.0    
##  [19] pkgdown_2.0.7          spatstat.sparse_3.0-0  colorspace_2.1-0      
##  [22] rappdirs_0.3.3         ggrepel_0.9.3          textshaping_0.3.6     
##  [25] xfun_0.37              dplyr_1.1.0            crayon_1.5.2          
##  [28] jsonlite_1.8.4         progressr_0.13.0       spatstat.data_3.0-0   
##  [31] survival_3.3-1         zoo_1.8-11             glue_1.6.2            
##  [34] polyclip_1.10-4        gtable_0.3.1           leiden_0.4.3          
##  [37] future.apply_1.10.0    abind_1.4-5            scales_1.2.1          
##  [40] spatstat.random_3.1-3  miniUI_0.1.1.1         Rcpp_1.0.10           
##  [43] viridisLite_0.4.1      xtable_1.8-4           reticulate_1.28       
##  [46] htmlwidgets_1.6.1      httr_1.4.5             RColorBrewer_1.1-3    
##  [49] ellipsis_0.3.2         ica_1.0-3              farver_2.1.1          
##  [52] pkgconfig_2.0.3        sass_0.4.5             uwot_0.1.14           
##  [55] deldir_1.0-6           utf8_1.2.3             tidyselect_1.2.0      
##  [58] labeling_0.4.2         rlang_1.0.6            reshape2_1.4.4        
##  [61] later_1.3.0            munsell_0.5.0          tools_4.2.0           
##  [64] cachem_1.0.7           cli_3.6.0              generics_0.1.3        
##  [67] ggridges_0.5.4         evaluate_0.20          stringr_1.5.0         
##  [70] fastmap_1.1.1          yaml_2.3.7             ragg_1.2.5            
##  [73] goftest_1.2-3          knitr_1.42             fs_1.6.1              
##  [76] fitdistrplus_1.1-8     purrr_1.0.1            RANN_2.6.1            
##  [79] pbapply_1.7-0          future_1.31.0          nlme_3.1-157          
##  [82] mime_0.12              formatR_1.14           compiler_4.2.0        
##  [85] plotly_4.10.1          png_0.1-8              spatstat.utils_3.0-1  
##  [88] tibble_3.1.8           bslib_0.4.2            stringi_1.7.12        
##  [91] highr_0.10             desc_1.4.2             lattice_0.20-45       
##  [94] Matrix_1.5-3           vctrs_0.5.2            pillar_1.8.1          
##  [97] lifecycle_1.0.3        spatstat.geom_3.0-6    lmtest_0.9-40         
## [100] jquerylib_0.1.4        RcppAnnoy_0.0.20       data.table_1.14.8     
## [103] irlba_2.3.5.1          httpuv_1.6.9           R6_2.5.1              
## [106] promises_1.2.0.1       KernSmooth_2.23-20     gridExtra_2.3         
## [109] parallelly_1.34.0      codetools_0.2-18       MASS_7.3-56           
## [112] rprojroot_2.0.3        withr_2.5.0            sctransform_0.3.5     
## [115] parallel_4.2.0         grid_4.2.0             tidyr_1.3.0           
## [118] rmarkdown_2.20         Rtsne_0.16             spatstat.explore_3.0-6
## [121] shiny_1.7.4