A common application of single-cell RNA sequencing (RNA-seq) data is
to identify discrete cell types. To take advantage of the large collection
of well-annotated scRNA-seq datasets, scClassify
package implements
a set of methods to perform accurate cell type classification based on
ensemble learning and sample size calculation.
This vignette demonstrates the usage of scClassify
,
providing a pithy description of each method with workable examples.
First, install scClassify
via BiocManager
.
# installation of scClassify
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("scClassify")
We assume that you have log-transformed (size-factor normalized) matrices where each row is a gene and each column a cell for a reference dataset and a query dataset. For demonstration purposes, we will take a subset of single-cell pancreas datasets from two independent studies (Wang et al., and Xin et al.).
library("scClassify")
data("scClassify_example")
xin_cellTypes <- scClassify_example$xin_cellTypes
exprsMat_xin_subset <- scClassify_example$exprsMat_xin_subset
wang_cellTypes <- scClassify_example$wang_cellTypes
exprsMat_wang_subset <- scClassify_example$exprsMat_wang_subset
exprsMat_xin_subset <- as(exprsMat_xin_subset, "dgCMatrix")
exprsMat_wang_subset <- as(exprsMat_wang_subset, "dgCMatrix")
The original cell type annotations and compositions of the example datasets can be easily accessed as shown below.
table(xin_cellTypes)
#> xin_cellTypes
#> alpha beta delta gamma
#> 285 261 49 79
table(wang_cellTypes)
#> wang_cellTypes
#> acinar alpha beta delta ductal gamma stellate
#> 5 206 118 10 96 21 45
We can see that Xin et al. data only have 4 cell types, while Wang et al. has 7 cell types.
We first perform non-ensemble scClassify
by using Xin et al.
as our reference dataset and Wang et al. data as ur query data.
We use WKNN
as the KNN algorithm, DE
(differential expression genes)
as the gene selection method, and lastly pearson
as
the similarity calculation method.
scClassify_res <- scClassify(exprsMat_train = exprsMat_xin_subset,
cellTypes_train = xin_cellTypes,
exprsMat_test = list(wang = exprsMat_wang_subset),
cellTypes_test = list(wang = wang_cellTypes),
tree = "HOPACH",
algorithm = "WKNN",
selectFeatures = c("limma"),
similarity = c("pearson"),
returnList = FALSE,
verbose = FALSE)
We can check the cell type tree generated by the reference data:
scClassify_res$trainRes
#> Class: scClassifyTrainModel
#> Model name: training
#> Feature selection methods: limma
#> Number of cells in the training data: 674
#> Number of cell types in the training data: 4
plotCellTypeTree(cellTypeTree(scClassify_res$trainRes))
Noted that scClassify_res$trainRes
is a scClassifyTrainModel
class.
Check the prediction results.
table(scClassify_res$testRes$wang$pearson_WKNN_limma$predRes, wang_cellTypes)
#> wang_cellTypes
#> acinar alpha beta delta ductal gamma stellate
#> alpha 0 206 0 0 0 2 0
#> beta 0 0 118 0 1 0 0
#> beta_delta_gamma 0 0 0 0 25 0 0
#> delta 0 0 0 10 0 0 0
#> gamma 0 0 0 0 0 19 0
#> unassigned 5 0 0 0 70 0 45
We next perform ensemble scClassify
by using Xin et al.
as our reference dataset and Wang et al. data as our query data.
We use WKNN
as the KNN algorithm, DE
as the gene selection method,
and pearson
and spearman
as the similarity calculation methods.
Thus, we will generate two combinations of gene selection models and
similarity metrics as training classifiers:
WKNN
+ DE
+ pearson
WKNN
+ DE
+ spearman
Here, we will weight these two classifiers equally by setting
weighted_ensemble = FALSE
. By default this is set as TRUE
,
so each base classifier will be weighted by the accuracy rates trained
in the reference data.
scClassify_res_ensemble <- scClassify(exprsMat_train = exprsMat_xin_subset,
cellTypes_train = xin_cellTypes,
exprsMat_test = list(wang = exprsMat_wang_subset),
cellTypes_test = list(wang = wang_cellTypes),
tree = "HOPACH",
algorithm = "WKNN",
selectFeatures = c("limma"),
similarity = c("pearson", "cosine"),
weighted_ensemble = FALSE,
returnList = FALSE,
verbose = FALSE)
We can compare the two base classifiers predictions as below.
table(scClassify_res_ensemble$testRes$wang$pearson_WKNN_limma$predRes,
scClassify_res_ensemble$testRes$wang$cosine_WKNN_limma$predRes)
#>
#> alpha beta beta_delta_gamma delta gamma unassigned
#> alpha 208 0 0 0 0 0
#> beta 0 119 0 0 0 0
#> beta_delta_gamma 0 0 22 0 0 3
#> delta 0 0 0 10 0 0
#> gamma 0 0 0 0 19 0
#> unassigned 1 3 15 0 0 101
Now, check the final ensemble results:
table(scClassify_res_ensemble$testRes$wang$ensembleRes$cellTypes,
wang_cellTypes)
#> wang_cellTypes
#> acinar alpha beta delta ductal gamma stellate
#> alpha 0 206 0 0 0 2 0
#> beta 0 0 118 0 1 0 0
#> beta_delta_gamma 0 0 0 0 25 0 0
#> delta 0 0 0 10 0 0 0
#> gamma 0 0 0 0 0 19 0
#> unassigned 5 0 0 0 70 0 45
You can also train your own model scClassifyTrainModel
using
train_scClassify()
. Note that by setting weightsCal = TRUE
,
we will calculate the training error of the reference data as
the weights for the individual classifiers.
Here, we illustrate the training function with gene selection methods based on differential expression (“limma”) and biomodal distribution (“BI”).
trainClass <- train_scClassify(exprsMat_train = exprsMat_xin_subset,
cellTypes_train = xin_cellTypes,
selectFeatures = c("limma", "BI"),
returnList = FALSE
)
#> after filtering not expressed genes
#> [1] 980 674
#> [1] "Feature Selection..."
#> [1] "Number of genes selected to construct HOPACH tree 160"
#> [1] "Constructing tree ..."
#> [1] "Training...."
#> [1] "=== selecting features by: limma ===="
#> [1] "=== selecting features by: BI ===="
trainClass
#> Class: scClassifyTrainModel
#> Model name: training
#> Feature selection methods: limma BI
#> Number of cells in the training data: 674
#> Number of cell types in the training data: 4
sessionInfo()
#> R version 4.4.0 RC (2024-04-16 r86468)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [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
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] scClassify_1.17.0 BiocStyle_2.33.1
#>
#> loaded via a namespace (and not attached):
#> [1] gridExtra_2.3 rlang_1.1.4
#> [3] magrittr_2.0.3 matrixStats_1.3.0
#> [5] compiler_4.4.0 mgcv_1.9-1
#> [7] DelayedMatrixStats_1.27.1 vctrs_0.6.5
#> [9] reshape2_1.4.4 stringr_1.5.1
#> [11] pkgconfig_2.0.3 crayon_1.5.3
#> [13] fastmap_1.2.0 magick_2.8.3
#> [15] XVector_0.45.0 labeling_0.4.3
#> [17] ggraph_2.2.1 utf8_1.2.4
#> [19] rmarkdown_2.27 UCSC.utils_1.1.0
#> [21] tinytex_0.51 purrr_1.0.2
#> [23] xfun_0.45 zlibbioc_1.51.1
#> [25] cachem_1.1.0 GenomeInfoDb_1.41.1
#> [27] jsonlite_1.8.8 highr_0.11
#> [29] rhdf5filters_1.17.0 DelayedArray_0.31.3
#> [31] Rhdf5lib_1.27.0 BiocParallel_1.39.0
#> [33] tweenr_2.0.3 parallel_4.4.0
#> [35] cluster_2.1.6 R6_2.5.1
#> [37] bslib_0.7.0 stringi_1.8.4
#> [39] limma_3.61.2 diptest_0.77-1
#> [41] GenomicRanges_1.57.1 jquerylib_0.1.4
#> [43] Rcpp_1.0.12 bookdown_0.39
#> [45] SummarizedExperiment_1.35.0 knitr_1.47
#> [47] mixtools_2.0.0 IRanges_2.39.0
#> [49] Matrix_1.7-0 splines_4.4.0
#> [51] igraph_2.0.3 tidyselect_1.2.1
#> [53] abind_1.4-5 yaml_2.3.8
#> [55] hopach_2.65.0 viridis_0.6.5
#> [57] codetools_0.2-20 minpack.lm_1.2-4
#> [59] Cepo_1.11.1 lattice_0.22-6
#> [61] tibble_3.2.1 plyr_1.8.9
#> [63] Biobase_2.65.0 withr_3.0.0
#> [65] evaluate_0.24.0 survival_3.7-0
#> [67] proxy_0.4-27 polyclip_1.10-6
#> [69] kernlab_0.9-32 pillar_1.9.0
#> [71] BiocManager_1.30.23 MatrixGenerics_1.17.0
#> [73] stats4_4.4.0 plotly_4.10.4
#> [75] generics_0.1.3 S4Vectors_0.43.0
#> [77] ggplot2_3.5.1 sparseMatrixStats_1.17.2
#> [79] munsell_0.5.1 scales_1.3.0
#> [81] glue_1.7.0 lazyeval_0.2.2
#> [83] proxyC_0.4.1 tools_4.4.0
#> [85] data.table_1.15.4 graphlayouts_1.1.1
#> [87] tidygraph_1.3.1 rhdf5_2.49.0
#> [89] grid_4.4.0 tidyr_1.3.1
#> [91] colorspace_2.1-0 SingleCellExperiment_1.27.2
#> [93] nlme_3.1-165 GenomeInfoDbData_1.2.12
#> [95] patchwork_1.2.0 ggforce_0.4.2
#> [97] HDF5Array_1.33.3 cli_3.6.3
#> [99] fansi_1.0.6 segmented_2.1-0
#> [101] S4Arrays_1.5.1 viridisLite_0.4.2
#> [103] dplyr_1.1.4 gtable_0.3.5
#> [105] sass_0.4.9 digest_0.6.36
#> [107] BiocGenerics_0.51.0 SparseArray_1.5.10
#> [109] ggrepel_0.9.5 htmlwidgets_1.6.4
#> [111] farver_2.1.2 memoise_2.0.1
#> [113] htmltools_0.5.8.1 lifecycle_1.0.4
#> [115] httr_1.4.7 statmod_1.5.0
#> [117] MASS_7.3-61