# scClassify Model Building and Prediction

Yingxin Lin1,2

1School of Mathematics and Statistics, The University of Sydney, Australia
2Charles Perkins Centre, The University of Sydney, Australia

#### 1 November 2022

# 1 Introduction

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")

# 2 Setting up the data

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.

# 3 scClassify

## 3.1 Non-ensemble scClassify

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

## 3.2 Ensemble Classify

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:

1. `WKNN` + `DE` + `pearson`
2. `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

# 4 Train your own model

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

# 5 Session Info

    sessionInfo()
    #> R version 4.2.1 (2022-06-23 ucrt)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server x64 (build 20348)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.utf8
    #> [3] LC_MONETARY=English_United States.utf8
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.utf8
    #>
    #> attached base packages:
    #> [1] stats     graphics  grDevices utils     datasets  methods   base
    #>
    #> other attached packages:
    #> [1] scClassify_1.10.0 BiocStyle_2.26.0
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] segmented_1.6-0             nlme_3.1-160
    #>   [3] bitops_1.0-7                matrixStats_0.62.0
    #>   [5] hopach_2.58.0               GenomeInfoDb_1.34.0
    #>   [7] tools_4.2.1                 bslib_0.4.0
    #>   [9] utf8_1.2.2                  R6_2.5.1
    #>  [11] HDF5Array_1.26.0            mgcv_1.8-41
    #>  [13] DBI_1.1.3                   BiocGenerics_0.44.0
    #>  [15] colorspace_2.0-3            rhdf5filters_1.10.0
    #>  [17] withr_2.5.0                 tidyselect_1.2.0
    #>  [19] gridExtra_2.3               proxyC_0.3.3
    #>  [21] compiler_4.2.1              cli_3.4.1
    #>  [23] Biobase_2.58.0              DelayedArray_0.24.0
    #>  [25] labeling_0.4.2              bookdown_0.29
    #>  [27] sass_0.4.2                  diptest_0.76-0
    #>  [29] scales_1.2.1                proxy_0.4-27
    #>  [31] stringr_1.4.1               digest_0.6.30
    #>  [33] mixtools_1.2.0              rmarkdown_2.17
    #>  [35] XVector_0.38.0              pkgconfig_2.0.3
    #>  [37] htmltools_0.5.3             sparseMatrixStats_1.10.0
    #>  [39] Cepo_1.4.0                  MatrixGenerics_1.10.0
    #>  [41] highr_0.9                   fastmap_1.1.0
    #>  [43] limma_3.54.0                rlang_1.0.6
    #>  [45] DelayedMatrixStats_1.20.0   jquerylib_0.1.4
    #>  [47] generics_0.1.3              farver_2.1.1
    #>  [49] jsonlite_1.8.3              BiocParallel_1.32.0
    #>  [51] dplyr_1.0.10                RCurl_1.98-1.9
    #>  [53] magrittr_2.0.3              GenomeInfoDbData_1.2.9
    #>  [55] patchwork_1.1.2             Matrix_1.5-1
    #>  [57] Rcpp_1.0.9                  munsell_0.5.0
    #>  [59] S4Vectors_0.36.0            Rhdf5lib_1.20.0
    #>  [61] fansi_1.0.3                 viridis_0.6.2
    #>  [63] lifecycle_1.0.3             stringi_1.7.8
    #>  [65] yaml_2.3.6                  ggraph_2.1.0
    #>  [67] MASS_7.3-58.1               SummarizedExperiment_1.28.0
    #>  [69] zlibbioc_1.44.0             rhdf5_2.42.0
    #>  [71] plyr_1.8.7                  grid_4.2.1
    #>  [73] parallel_4.2.1              ggrepel_0.9.1
    #>  [75] lattice_0.20-45             splines_4.2.1
    #>  [77] graphlayouts_0.8.3          magick_2.7.3
    #>  [79] knitr_1.40                  pillar_1.8.1
    #>  [81] igraph_1.3.5                GenomicRanges_1.50.0
    #>  [83] reshape2_1.4.4              codetools_0.2-18
    #>  [85] stats4_4.2.1                glue_1.6.2
    #>  [87] evaluate_0.17               RcppParallel_5.1.5
    #>  [89] BiocManager_1.30.19         vctrs_0.5.0
    #>  [91] tweenr_2.0.2                gtable_0.3.1
    #>  [93] purrr_0.3.5                 polyclip_1.10-4
    #>  [95] tidyr_1.2.1                 kernlab_0.9-31
    #>  [97] assertthat_0.2.1            cachem_1.0.6
    #>  [99] ggplot2_3.3.6               xfun_0.34
    #> [101] ggforce_0.4.1               tidygraph_1.2.2
    #> [103] survival_3.4-0              viridisLite_0.4.1
    #> [105] minpack.lm_1.2-2            SingleCellExperiment_1.20.0
    #> [107] tibble_3.1.8                IRanges_2.32.0
    #> [109] cluster_2.1.4               statmod_1.4.37