clustifyr requires query and reference data in the form of normalized expression matrices, corresponding metadata tables, and a list of variable genes ( Figure 1).

clustifyr adopts correlation-based methods to find reference transcriptomes with the highest similarity to query cluster expression profiles, defaulting to Spearman ranked correlation, with options to use Pearson, Kendall, or Cosine correlation instead if desired. clustify() will return a matrix of correlation coefficients for each cell type and cluster, with the row names corresponding to the query cluster number and column names as the reference cell types.

Query clusters are assigned cell types to the highest correlated reference cell type, with an automatic or manual cutoff threshold. Query clusters dissimilar to all available reference cell types are labeled as “unassigned”.

To better integrate with standard workflows that involve S3/S4 R objects, methods for clustifyr are written to directly recognize Seurat ¹⁴ (v2 and v3) and SingleCellExperiment ¹⁵ objects, retrieve the required information, and reinsert classification results back into an output object. A more general wrapper is also included for compatibility with other common data structures and can be easily extended to new object types. This approach also has the added benefit of forgoing certain calculations such as variable gene selection or clustering, which may already be stored within input objects.

res <- clustify( input = sce_small, # an SCE object ref_mat = cbmc_ref , # matrix of expression for each cell type cluster_col = "cell_type1" , # column in meta.data with clusters obj_out = TRUE # output SCE object with cell type ) SingleCellExperiment :: colData(res) [ 1 : 10 , c( "type" , "r" )] #> DataFrame with 10 rows and 2 columns #> type r #> <character> <numeric> #> AZ_A1 pDCs 0.814336567702192 #> AZ_A10 Eryth 0.665800619720566 #> AZ_A11 pDCs 0.682088309107356 #> AZ_A12 Eryth 0.665800619720566 #> AZ_A2 B 0.634114583333333 #> AZ_A3 pDCs 0.814336567702192 #> AZ_A4 pDCs 0.814336567702192 #> AZ_A5 NK 0.655407634437123 #> AZ_A6 pDCs 0.682088309107356 #> AZ_A7 pDCs 0.71424223704931

In the absence of suitable reference data (i.e. RNA-seq or microarray expression matrices), clustifyr can build scRNA-seq reference data by averaging per-cell expression data for each cluster, to generate a transcriptomic snapshot. Direct reference-building from SingleCellExperiment or Seurat objects is supported as well.

new_ref_matrix <- average_clusters( mat = pbmc_matrix_small, metadata = pbmc_meta $ classified, # or use metadata = pbmc_meta, cluster_col = "classified" if_log = TRUE # whether the expression matrix is already log transformed ) new_ref_matrix_sce <- object_ref( input = sce_small, # SCE object cluster_col = "cell_type1" # column in colData with cell identities ) new_ref_matrix_v3 < - seurat_ref( seurat_object = s_small3, # SeuratV3 object cluster_col = "RNA_snn_res.1" # column in meta.data with cell identities )

Data exploration plotting functions, for dimensional reduction scatter plots and heatmaps, are extended from ggplot2 and ComplexHeatmap packages, featuring colorblind-friendly default colors. Gene list-based methods (clustify_lists()) are also implemented via hypergeometric tests, GSEA, jaccard index, or percentage gene detection by cluster, which provide easy to interpret methods to verify the presence of known positive and negative marker genes.

Reference datasets. Multiple scRNA-seq and other cell type references datasets are provided in an ExperimentHub Bioconductor package (clustifyrdatahub). A description of these datasets and others used for benchmarking and optimizing parameters for clustifyr are provided in Table 1.

Correlation method. We benchmarked clustifyr against a suite of comparable datasets, PBMC-bench ^{13, 16} , generated using multiple scRNA-seq methods on aliquots of peripheral blood mononuclear cells (PBMCs) from two individuals. Additional details about each query and reference dataset are provided in Supplemental Table 1. For each single cell technology, average gene expression profiles were generated from annotated cell types and compared across each platform. Notably, for each reference dataset cross-referenced against all other samples, clustifyr achieved a median F1-score (see Benchmarking Methods) of above 0.94 using Spearman ranked correlation ( Figure 2A). Other correlation methods are on par or slightly worse at cross-platform classifications, which is expected based on the nature of ranked vs unranked methods. We therefore selected Spearman as the default method in clustifyr, with other methods also available, as well as a wrapper function to find consensus identities across available correlation methods (call_consensus()).

Correlation minimum cutoff. Recognition of missing reference cell types, so as to avoid misclassification, is another point of great interest in the field. From general usage of clustifyr, we find using a minimum correlation cutoff of 0.5 or 0.4 is generally satisfactory. Alternatively, the cutoff threshold can be determined heuristically using 0.8 * highest correlation coefficient among the clusters. One example is shown in Figure 2B, using PBMC rejection benchmark data modified by the SciBet package ¹⁷ . Megakaryocytes were removed from the reference melanoma immune cells data, but retained in the test data to mimic the situation when the reference data does not contain a rare cell type. clustifyr analysis successfully found the megakaryocytes to be dissimilar to all available reference cell types, and hence left as “unassigned” under the default minimum threshold cutoff.

Variable gene selection and normalization. As the core function of clustifyr is ranked correlation, feature selection to focus on highly variable genes is critical. To illustrate the importance of feature selection we used clustifyr to classify pancreatic cell types generated using the inDrops platform using a reference built from a dataset generated on the Smart-Seq2 platform ^{18, 19} . In Figure 2C, we compare correlation coefficients using all detected genes (>10,000) vs feature selection by the package M3Drop. A basic level of feature selection, e.g. using M3Drop, Seurat VST (default uses top 2,000 variable genes), or simply 1,000 genes with the highest variance in the reference data, is sufficient to classify the pancreatic cells. In the case of other cell type mixtures, especially ones without complete knowledge of the expected cell types, further optimization of clustering and feature selection may be of greater importance. clustifyr does not provide novel clustering, feature selection, or normalization methods on its own, but instead is built to maintain flexibility to incorporate methods from other, and future, packages. We recommend that users use normalized reference and query data and match normalization methods between datasets when possible. We view these questions as fast-moving fields ^{20, 21} , and hope to benefit from new advances, while keeping the general clustifyr framework intact.

Minimum cells per cluster. We next applied clustifyr to a larger general reference set built from the Mouse Cell Atlas ²² and examined cell type classification of another mouse cell atlas, the Tabula Muris dataset ⁵ . clustifyr assigned cell types with a median accuracy of 1. Using these test datasets we sought to determine the minimum number of query cells necessary in a cluster to obtain accurate cell type annotation. We subsampled the query data ( Figure 2D) and as expected, with further downsampling of the number of cells in each query cluster, we observe decreased accuracy. Yet, even at 15 cells per tested cluster, clustifyr still performed well, with a further increase in speed. Based on these results, we set the default parameters in clustifyr to exclude or warn users of classification on clusters containing less than 10 cells. These results also suggest that clustering the query dataset to obtain more refined clusters (e.g fewer cells per cluster) could be employed to aid in the identification of rarer or less well-defined cell subsets. clustifyr can also be used to classify individual cells, although we do not recommend per cell classification because of the reduced accuracy observed with decreasing numbers of cells per cluster.

Subclustering. clustifyr also provides functionality to assess the quality of the cell type annotations. An intentional overclustering and classification function based on k-means clustering (overcluster_test()) is implemented in clustifyr for exploration of cell type annotation at increasing numbers of clusters ( Figure 2E). This approach provides a rapid visualization to determine if cell type annotations are stable with varying numbers of clusters. For example, scRNA-seq data from the Seurat PBMC 3k tutorial was reclassified at multiple clustering levels using Cord Blood Mononuclear Cells (CBMCs) as reference, which demonstrated largely stable cell type assignments in the presence of overclustered query data ( Figure 2E) ²³ . When using scRNA-seq data as the reference data, matrices are built by averaging per-cell expression data for each cluster (average_clusters()), to generate a transcriptomic snapshot similar to bulk RNA-seq or microarray data. An additional argument to subcluster the reference single cell clusters is also available, to generate more than one expression profile per reference cell type, in a manner analogous to overcluster_test(), but applied to the reference scRNA-seq dataset. The number of subclusters for each reference cell type is dependent on the number of cells in the cluster (n), and the sub-clustering power argument (x), following the formula n ^{x 9} . This approach does not improve classification in the PBMC-bench data ( Figure 2F), whose reference and query clustering are already consistent. However, we envision its utility would greatly depend on the granularity of the clustering in the reference dataset.

Using clustifyr, PBMC clusters from the Seurat PBMC 3k tutorial are correctly labeled using either bulk-RNA seq references generated from processed microarray data of purified cell types ²⁴ , the ImmGen database of bulk-RNA-seq ^{9, 25} , or previously annotated scRNA-seq results from the Seurat CBMC CITE-seq tutorial ^{14, 23} ( Figure 3).

To assess the performance of clustifyr, we used the Tabula Muris dataset ⁵ , which contains data generated from 12 matching tissues using both 10x Genomics 3’ end seq (“drop”) and Smart-Seq2 (“facs”) platforms. We attempted to assign cell type identities to clusters in “drop” Seurat objects using references built from “facs” Seurat objects, which contain pre-computed variable genes generated by the Seurat mean.var.plot (dispersion z-scores based on expression bins) approach. For each method we used the recommended variable gene selection approach. clustifyr uses variable genes supplied by the user and for benchmarking we used the variable genes stored in the Seurat object. scmap calculates variable genes using a modified approach based on M3drop. SingleR selects variable genes by identifying marker genes between clusters. scPred, in contrast, selects informative principal components as a feature selection procedure whereas ACTINN does not perform feature selection for classification.

In benchmarking results, clustifyr is comparably accurate versus other automated classification packages ( Figure 4A). Cross-platform comparisons are inherently more difficult, and the approach used by clustifyr is aimed at being platform- and normalization-agnostic. Mean runtime, including both reference building and test data classification, in Tabular Muris classifications was ~ 1 second if the required variable gene list is extracted from the query Seurat object. Alternatively, variable genes can be recalculated by other methods such as M3Drop ²⁶ , to reach similar results (clustifyr (m3drop)).

Signature marker gene lists are an additional reference data type that is commonly used to guide cluster cell type classification. We therefore sought to determine if a gene list enrichment approach could provide comparable classification power to using correlation. clustifyr provides a function clustify_lists() which compares marker genes between query clusters to a list of marker genes per reference cell type. clustify_lists will calculate enrichment with a hypergeometric test, marker overlap with the jaccard index, or use the percent of cells expressing marker genes to annotate cell types. Alternatively, if ranked gene lists are available, Gene Set Enrichment Analysis (GSEA) using the fGSEA package ²⁷ or Spearman ranked correlation can be employed. We find that using gene expression for clustifyr classification had higher accuracy than gene list enrichment using a hypergeometric test or the jaccard index, however this approach could be very useful for datasets without scRNA or bulk RNA-seq data for use as a reference. ( Figure 4A).

For scalability benchmarking, we adapted an existing benchmark dataset, scRNAseq_Benchmark subsampling, which contains query and reference data with downsampled numbers of cells from the Smart-seq2 Tabula Muris dataset ^{5, 13} . Once again, clustifyr is accurate and efficient, compared to other developed methods ( Figure 4B). As a further comparison, we also examined classification of cell types in murine brain datasets generated by the Allen Institute Brain Atlas, and provided by the scRNAseq_Benchmark pipeline ¹³ . The two murine brain regions contained 34 shared cell types and clustifyr was also able to reach similarly satisfactory cell annotation compared to other annotation methods. ( Figure 4C).

Lastly, we applied clustifyr to a series of increasingly challenging datasets from the scRNAseq_Benchmark ¹³ unseen population rejection test ( Figure 4D). This test assess how frequently cells will be mis-assigned when the corresponded cell types are not present in the dataset. The PBMC dataset contains different T-cell subsets, which do not often cluster into discrete well-defined cell types solely based on gene expression. Without the corresponding cell type references, 57.5% of T cells were rejected and unassigned. When only CD4+ references were removed, 28.2% of test CD4+ T cells were rejected and unassigned. clustifyr was unable to reject CD4+/CD45RO+ memory T cells, mislabeling them as CD4+/CD25 T Reg instead when the exact reference was unavailable. However, these misclassifications are also observed with other classification tools benchmarked in the scRNAseq_Benchmark study ( Figure 4D) ¹³ .

clustifyr was tested against scmap v1.8.0 ⁸ , SingleR v1.0.1 ⁹ , Seurat v3.1.1 ¹⁴ , latest GitHub versions of ACTINN ¹¹ and scPred ¹² , and SVM as implemented in python3 scikit-learn v0.19.1 ²⁸ . scRNA-seq Tabula Muris data was downloaded as seuratV2 objects. Human pancreas data was downloaded as SCE objects. In all instances, to mimic the usage case of clustifyr, clustering and dimension reduction projections are acquired from available metadata, in lieu of new analysis.

An R script was modified to benchmark clustifyr following the approach and datasets of scRNAseq_Benchmark ¹³ , using M3Drop ²⁶ to generate variable genes for clustifyr. R code used for benchmarking, and preprocessing of other datasets, in the form of matrices and tables, are documented in R scripts available in the clustifyr and clustifyrdatahub GitHub repositories.

Classification accuracy was measured using two approaches depending on the datasets compared. For datasets where the query and reference data contain identical cell types, an F1-score, the harmonic mean of the precision and recall, was calculated for each cell type (PBMC-bench, Allen Brain Institute Atlas, and Smart-Seq2 Tabula Muris subsampling). When summarizing classification accuracy across an entire dataset the median F1-score is reported. Datasets with varying cell types in the query and reference data cannot be characterized with an F1-score and instead accuracy, defined as the ratio between the number of correctly classified clusters and the overall number of clusters, is reported (Mouse cell atlas vs. Tabula Muris and Tabula Muris Smart-Seq2 vs. 10x Genomics).

clustifyr is distributed as part of the Bioconductor R package repository and is compatible with Mac OS X, Windows, and major Linux operating systems. Package dependencies and system requirements are documented in the clustifyr Bioconductor repository.