1 Introduction

This tutorial shows how to apply CDI to select optimal clustering labels among different candidate cell labels for UMI counts of single-cell RNA-sequencing dataset. Sections 1-3 are for datasets from one batch; Section 4 is for datasets from multiple batches. Datasets used in this tutorial were simulated for illustration purpose.

2 Installation

CDI can be installed from Bioconductor:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("CDI")

The latest version can be installed from Github:

if (!requireNamespace("remotes", quietly = TRUE))
    install.packages("remotes")
remotes::install_github("jichunxie/CDI", build_vignettes = TRUE) 

3 Load datasets

CDI provides simulated toy single-cell RNA-seq UMI count matrices for illustration purpose. To use other datasets, please change the code in the following two blocks. Filtering of cells/genes can be applied before feature gene selection; normalization/imputation is not applicable here, as CDI models the raw UMI counts.

Inputs:

  • one_batch_matrix: a single-cell RNA-seq UMI count matrix from one batch (each row represents a gene and each column represent a cell). To apply CDI to datasets from multiple batches, please check section 4.
  • one_batch_matrix_label_df: the label sets of the one_batch_matrix, where each row represents a cell, and each column represents a set of cell labels generated by existing clustering methods such as K-Means.
  • one_batch_matrix_celltype: a vector of characters indicating the benchmark cell types of cells in one_batch_matrix (not necessary except for evaluation purpose).
data(one_batch_matrix, package = "CDI")
dim(one_batch_matrix)
## [1] 3000 2000
data(one_batch_matrix_celltype, package = "CDI")
table(one_batch_matrix_celltype)
## one_batch_matrix_celltype
## type1 type2 type3 type4 type5 
##   400   400   400   400   400

Here we provide 12 label sets generated from PCA\(+\)KMeans and Seurat v3.1.5, where the number of clusters range from 2 to 7.

For better visualization and comparison, the column names of this data frame indicate the clustering method and the number of clusters. For example, “KMeans_k5” refers to the set of cell labels generated by KMeans with 5 clusters.

data(one_batch_matrix_label_df, package = "CDI")
knitr::kable(head(one_batch_matrix_label_df[,c("KMeans_k2", "KMeans_k4", "Seurat_k2", "Seurat_k3")], 3))
KMeans_k2 KMeans_k4 Seurat_k2 Seurat_k3
2 3 1 2
2 3 1 2
2 3 1 2

4 Select feature genes with feature_gene_selection

Feature genes are those genes that express differently across cell types. Several methods are available for feature gene selection. Here we propose a new method called working dispersion score (WDS). We select genes with highest WDS as the feature genes.

feature_gene_indx <- feature_gene_selection(
    X = one_batch_matrix, 
    batch_label = NULL, 
    method = "wds", 
    nfeature = 500)
sub_one_batch_matrix <- one_batch_matrix[feature_gene_indx,]

5 calculate_CDI

First, calculate the size factor of each gene with size_factor. The output of size_factor will be one of the inputs of calculate_CDI.

one_batch_matrix_size_factor <- size_factor(X = one_batch_matrix)

Second, calculate CDI for each candidate set of cell labels:

start_time <- Sys.time()
CDI_return1 <- calculate_CDI(
    X = sub_one_batch_matrix, 
    cand_lab_df = one_batch_matrix_label_df, 
    batch_label = NULL, 
    cell_size_factor = one_batch_matrix_size_factor)
end_time <- Sys.time()
print(difftime(end_time, start_time))
## Time difference of 1.003124 mins

6 Select the optimal label set with minimum CDI-AIC/CDI-BIC

knitr::kable(CDI_return1)
Label_name Cluster_method CDI_AIC CDI_BIC neg_llk_val N_cluster
KMeans_k2 KMeans_k2 KMeans 1135287 1146489 565643.4 2
KMeans_k3 KMeans_k3 KMeans 1120096 1136899 557047.9 3
KMeans_k4 KMeans_k4 KMeans 1112561 1134964 552280.4 4
KMeans_k5 KMeans_k5 KMeans 1101970 1129974 545984.9 5
KMeans_k6 KMeans_k6 KMeans 1103150 1136756 545575.1 6
KMeans_k7 KMeans_k7 KMeans 1105022 1144229 545511.2 7
Seurat_k2 Seurat_k2 Seurat 1129063 1140264 562531.3 2
Seurat_k3 Seurat_k3 Seurat 1120002 1136805 557001.2 3
Seurat_k4 Seurat_k4 Seurat 1109115 1131518 550557.3 4
Seurat_k5 Seurat_k5 Seurat 1102321 1130325 546160.4 5
Seurat_k6 Seurat_k6 Seurat 1103152 1136751 545575.9 6
Seurat_k7 Seurat_k7 Seurat 1104099 1143292 545049.6 7

CDI counts the number of unique characters in each label set as the “N_cluster”. If the column names of label set data frame are provided with the format "[ClusteringMethod]_k[NumberOfClusters]" such as “KMeans_K5, calculate_CDI will extract the”[ClusteringMethod]" as the Cluster_method. The clustering method can also be provided in the argument clustering_method for each label set.

The outputs of calculate_CDI include CDI-AIC and CDI-BIC. CDI calculates AIC and BIC of cell-type-specific gene-specific NB model for UMI counts, where the cell types are based on each candidate label set, and only the selected subset of genes are considered. Whether to use CDI-AIC or CDI-BIC depend on the goals. We suggest using CDI-BIC to select optimal main cell types and using CDI-AIC to select optimal subtypes, because BIC puts more penalty on the complexity of models (number of clusters).

Output visualization

The outputs of CDI are demonstrated with a lineplot. The x-axis is for the number of clusters. The y-axis is for the CDI values. Different colors represent different clustering methods: The orange line represents label sets from Seurat; the blue line represents label sets from K-Means. The red triangle represents the optimal clustering result corresponding to the smallest CDI value.

The cell population in this simulated dataset doesn’t have a hierarchical structure. We use CDI-BIC to select the optimal set of cell labels. The label set “KMeans_k5” has the smallest CDI-BIC and is selected as the optimal.

CDI_lineplot(cdi_dataframe = CDI_return1, cdi_type = "CDI_BIC")

A benchmark label set refers to human annotated cell type label set from biomarkers for real data or true label set for simulated data. If such label set is available, we can also compare the optimal clustering label set selected by CDI with it. Here we provide the heatmap of contingency table between benchmark label set and the optimal clustering label set selected by CDI. Each row represents a cell type in the benchmark label set, and each column represents a cluster in the optimal clustering label set selected by CDI. Each rectangle is color-scaled by the proportion of the cells in the given cluster coming from the benchmark types. Each column sums to 1.

contingency_heatmap(benchmark_label = one_batch_matrix_celltype,
    candidate_label = one_batch_matrix_label_df$KMeans_k5,
    rename_candidate_clusters = TRUE,
    candidate_cluster_names = paste0('cluster', seq_len(length(unique(one_batch_matrix_label_df$KMeans_k5)))))

We can also calculate the CDI for the benchmark label set as the following:

benchmark_return1 <- calculate_CDI(X = sub_one_batch_matrix,
    cand_lab_df = one_batch_matrix_celltype, 
    batch_label = NULL, 
    cell_size_factor = one_batch_matrix_size_factor)

The results of benchmark label set can be added to the CDI lineplot:

CDI_lineplot(cdi_dataframe = CDI_return1,
    cdi_type = "CDI_BIC",
    benchmark_celltype_cdi = benchmark_return1,
    benchmark_celltype_ncluster = length(unique(one_batch_matrix_celltype)))

The purple star represents the CDI value for the benchmark label set. The result shows that the benchmark label set achieves lowest CDI-BIC among all label sets. Therefore, if the benchmark label set is available as one candidate label set, it will be selected by CDI-BIC as the optimal.

7 Apply CDI to datasets from multiple batches

CDI accepts datasets from multiple batches. The candidate label sets can be obtained after batch effect correction, but CDI will be calculated based on raw UMI count matrices. The procedure of applying CDI to multi-batch datasets is similar to that of one batch, except for one extra input of batch labels.

7.1 Inputs:

  • two_batch_matrix: a single-cell RNA-seq UMI count matrix from the multiple batches (each row represents a gene and each column represent a cell).
  • two_batch_matrix_batch: a vector of characters indicating the batch labels of cells in two_batch_matrix. The number of batches can be greater than 2.
  • two_batch_matrix_label_df: the label sets of the two_batch_matrix, where each row represents a cell, and each column represents a label set generated by existing clustering methods such as K-Means after batch effect correction.
  • two_batch_matrix_celltype: a vector of characters indicating the benchmark cell types of cells in two_batch_matrix (not necessary unless for evaluation purpose).

7.2 Data loading

data(two_batch_matrix_celltype, package = "CDI")
table(two_batch_matrix_celltype)
## two_batch_matrix_celltype
## type1 type2 type3 type4 type5 
##   400   400   400   400   400
data(two_batch_matrix_batch, package = "CDI")
table(two_batch_matrix_batch)
## two_batch_matrix_batch
## batch1 batch2 
##   1000   1000
data(two_batch_matrix, package = "CDI")
dim(two_batch_matrix)
## [1] 3000 2000
data(two_batch_matrix_label_df, package = "CDI")
knitr::kable(head(two_batch_matrix_label_df[,c("KMeans_k5", "KMeans_k6", "Seurat_k5", "Seurat_k6")], 3))
KMeans_k5 KMeans_k6 Seurat_k5 Seurat_k6
3 3 2 2
3 3 2 2
3 3 2 2

7.3 Feature gene selection

feature_gene_indx <- feature_gene_selection(
    X = two_batch_matrix,
    batch_label = two_batch_matrix_batch,
    method = "wds",
    nfeature = 500)
sub_two_batch_matrix <- two_batch_matrix[feature_gene_indx,]

7.4 Calculate CDI

two_batch_matrix_size_factor <- size_factor(two_batch_matrix)

start_time <- Sys.time()
CDI_return2 <- calculate_CDI(
    X = sub_two_batch_matrix,
    cand_lab_df = two_batch_matrix_label_df, 
    batch_label = two_batch_matrix_batch,
    cell_size_factor = two_batch_matrix_size_factor)
end_time <- Sys.time()
print(difftime(end_time, start_time))
## Time difference of 5.816393 mins
knitr::kable(CDI_return2)
Label_name Cluster_method CDI_AIC CDI_BIC neg_llk_val N_cluster
KMeans_k3 KMeans_k3 KMeans 1116010 1133295 554919.1 3
KMeans_k4 KMeans_k4 KMeans 1090471 1114398 540963.7 4
KMeans_k5 KMeans_k5 KMeans 1094696 1124000 542115.9 5
KMeans_k6 KMeans_k6 KMeans 1084735 1119774 536111.6 6
KMeans_k7 KMeans_k7 KMeans 1085682 1125975 535646.8 7
KMeans_k8 KMeans_k8 KMeans 1086246 1132224 534910.8 8
KMeans_k9 KMeans_k9 KMeans 1086179 1137386 533943.5 9
KMeans_k10 KMeans_k10 KMeans 1089714 1146373 534741.1 10
KMeans_k11 KMeans_k11 KMeans 1090071 1151995 533979.7 11
KMeans_k12 KMeans_k12 KMeans 1091465 1159185 533640.6 12
KMeans_k13 KMeans_k13 KMeans 1093719 1167259 533729.7 13
KMeans_k14 KMeans_k14 KMeans 1093773 1172530 532822.4 14
KMeans_k15 KMeans_k15 KMeans 1097256 1182247 533451.8 15
Seurat_k3 Seurat_k3 Seurat 1098939 1117153 546217.3 3
Seurat_k4 Seurat_k4 Seurat 1090520 1114458 540986.0 4
Seurat_k5 Seurat_k5 Seurat 1083559 1113412 536449.6 5
Seurat_k6 Seurat_k6 Seurat 1084476 1119538 535978.0 6
Seurat_k7 Seurat_k7 Seurat 1084354 1124707 534970.9 7
Seurat_k8 Seurat_k8 Seurat 1086408 1132145 535037.9 8
Seurat_k9 Seurat_k9 Seurat 1087299 1138390 534527.4 9
Seurat_k10 Seurat_k10 Seurat 1087046 1143337 533471.2 10
Seurat_k11 Seurat_k11 Seurat 1089098 1150966 533501.1 11
Seurat_k12 Seurat_k12 Seurat 1090061 1157517 532984.3 12
Seurat_k13 Seurat_k13 Seurat 1090988 1164020 532449.9 13
Seurat_k14 Seurat_k14 Seurat 1091273 1170009 531572.7 14
Seurat_k15 Seurat_k15 Seurat 1091733 1176071 530800.4 15

Calculate the benchmark label set (if available) CDI:

benchmark_return <- calculate_CDI(
    X = sub_two_batch_matrix,
    cand_lab_df = two_batch_matrix_celltype,
    batch_label = two_batch_matrix_batch,
    cell_size_factor = two_batch_matrix_size_factor)

Output visualization:

CDI_lineplot(cdi_dataframe = CDI_return2,
    cdi_type = "CDI_BIC",
    benchmark_celltype_cdi = benchmark_return,
    benchmark_celltype_ncluster = length(unique(two_batch_matrix_celltype))) 

Compare the optimal clustering label set selected by CDI with the benchmark cell type labels:

contingency_heatmap(
    benchmark_label = two_batch_matrix_celltype,
    candidate_label = two_batch_matrix_label_df$Seurat_k5,
    rename_candidate_clusters = TRUE,
    candidate_cluster_names = paste0('cluster', seq_len(length(unique(one_batch_matrix_label_df$Seurat_k5)))))

8 Session Information

sessionInfo()
## R version 4.3.2 Patched (2023-11-13 r85521)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.18-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] CDI_1.0.2        BiocStyle_2.30.0
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3          jsonlite_1.8.8             
##   [3] magrittr_2.0.3              magick_2.8.1               
##   [5] spatstat.utils_3.0-4        farver_2.1.1               
##   [7] rmarkdown_2.25              zlibbioc_1.48.0            
##   [9] vctrs_0.6.5                 ROCR_1.0-11                
##  [11] spatstat.explore_3.2-5      RCurl_1.98-1.13            
##  [13] S4Arrays_1.2.0              htmltools_0.5.7            
##  [15] SparseArray_1.2.2           sass_0.4.7                 
##  [17] sctransform_0.4.1           parallelly_1.36.0          
##  [19] KernSmooth_2.23-22          bslib_0.6.1                
##  [21] htmlwidgets_1.6.3           ica_1.0-3                  
##  [23] plyr_1.8.9                  plotly_4.10.3              
##  [25] zoo_1.8-12                  cachem_1.0.8               
##  [27] igraph_1.5.1                mime_0.12                  
##  [29] lifecycle_1.0.4             pkgconfig_2.0.3            
##  [31] Matrix_1.6-4                R6_2.5.1                   
##  [33] fastmap_1.1.1               GenomeInfoDbData_1.2.11    
##  [35] MatrixGenerics_1.14.0       fitdistrplus_1.1-11        
##  [37] future_1.33.0               shiny_1.8.0                
##  [39] digest_0.6.33               colorspace_2.1-0           
##  [41] patchwork_1.1.3             S4Vectors_0.40.2           
##  [43] Seurat_5.0.1                tensor_1.5                 
##  [45] RSpectra_0.16-1             irlba_2.3.5.1              
##  [47] GenomicRanges_1.54.1        labeling_0.4.3             
##  [49] progressr_0.14.0            fansi_1.0.5                
##  [51] spatstat.sparse_3.0-3       httr_1.4.7                 
##  [53] polyclip_1.10-6             abind_1.4-5                
##  [55] compiler_4.3.2              withr_2.5.2                
##  [57] BiocParallel_1.36.0         fastDummies_1.7.3          
##  [59] highr_0.10                  MASS_7.3-60                
##  [61] DelayedArray_0.28.0         ggsci_3.0.0                
##  [63] tools_4.3.2                 lmtest_0.9-40              
##  [65] httpuv_1.6.12               future.apply_1.11.0        
##  [67] goftest_1.2-3               glue_1.6.2                 
##  [69] nlme_3.1-164                promises_1.2.1             
##  [71] grid_4.3.2                  Rtsne_0.16                 
##  [73] cluster_2.1.6               reshape2_1.4.4             
##  [75] generics_0.1.3              gtable_0.3.4               
##  [77] spatstat.data_3.0-3         tidyr_1.3.0                
##  [79] data.table_1.14.8           XVector_0.42.0             
##  [81] sp_2.1-2                    utf8_1.2.4                 
##  [83] BiocGenerics_0.48.1         spatstat.geom_3.2-7        
##  [85] RcppAnnoy_0.0.21            ggrepel_0.9.4              
##  [87] RANN_2.6.1                  pillar_1.9.0               
##  [89] stringr_1.5.1               spam_2.10-0                
##  [91] RcppHNSW_0.5.0              later_1.3.1                
##  [93] splines_4.3.2               dplyr_1.1.4                
##  [95] lattice_0.22-5              survival_3.5-7             
##  [97] deldir_2.0-2                tidyselect_1.2.0           
##  [99] SingleCellExperiment_1.24.0 miniUI_0.1.1.1             
## [101] pbapply_1.7-2               knitr_1.45                 
## [103] gridExtra_2.3               bookdown_0.37              
## [105] IRanges_2.36.0              SummarizedExperiment_1.32.0
## [107] scattermore_1.2             stats4_4.3.2               
## [109] xfun_0.41                   Biobase_2.62.0             
## [111] matrixStats_1.1.0           stringi_1.8.2              
## [113] lazyeval_0.2.2              yaml_2.3.7                 
## [115] evaluate_0.23               codetools_0.2-19           
## [117] tibble_3.2.1                BiocManager_1.30.22        
## [119] cli_3.6.1                   uwot_0.1.16                
## [121] xtable_1.8-4                reticulate_1.34.0          
## [123] munsell_0.5.0               jquerylib_0.1.4            
## [125] GenomeInfoDb_1.38.1         Rcpp_1.0.11                
## [127] globals_0.16.2              spatstat.random_3.2-2      
## [129] png_0.1-8                   parallel_4.3.2             
## [131] ellipsis_0.3.2              ggplot2_3.4.4              
## [133] dotCall64_1.1-1             bitops_1.0-7               
## [135] listenv_0.9.0               viridisLite_0.4.2          
## [137] scales_1.3.0                ggridges_0.5.4             
## [139] crayon_1.5.2                SeuratObject_5.0.1         
## [141] leiden_0.4.3.1              purrr_1.0.2                
## [143] rlang_1.1.2                 cowplot_1.1.1