Getting Start with Assisted CIDER (asCIDER)

Introduction

This vignette performs asCIDER, a meta-clustering method, on a cross-species pancreas dataset. AsCIDER is aimed to achieve the clustering task in data confounded by unwanted variables. In this scenario, the unwanted variable is species effects.

asCIDER is short for assisted CIDER, assisted by the prior batch-specific annotations or clustering results.

Set up

In addition to CIDER, we will load the following packages:

library(CIDER)
library(Seurat)
library(pheatmap)
library(ggplot2)
library(cowplot)
library(viridis)

Pancreas data

The example data can be downloaded from https://figshare.com/s/d5474749ca8c711cc205.

This data\(^1\) contain single-cell RNA-seq pancreatic data from human (8241 cells) and mouse (1886 cells) with 12474 rows (genes).

# Download the data
data_url <- "https://figshare.com/ndownloader/files/31469387"
data_file <- file.path(tempdir(), "pancreas_counts.RData")

if (!file.exists(data_file)) {
  message("Downloading count data to temporary directory...")
  download.file(data_url, destfile = data_file, mode = "wb")
}
#> Downloading count data to temporary directory...
# Load counts matrix and metadata
data_env <- new.env()
loaded_objects <- load(file = data_file, envir = data_env)
pancreas_counts <- data_env[[loaded_objects[1]]]

load("../data/pancreas_meta.RData") # meta data/cell information
seu <- CreateSeuratObject(counts = pancreas_counts, meta.data = pancreas_meta)
table(seu$Batch)
#> 
#> human mouse 
#>  8241  1886

Exam if data are confounded

Prior to use CIDER (or other integration methods for clustering) it is important to exam if clustering will be confounded by the cross-species factors.

Her we first perform the conventional Seurat\(^2\) clustering pipeline.

seu <- NormalizeData(seu, verbose = FALSE)
seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
seu <- ScaleData(seu, verbose = FALSE)
seu <- RunPCA(seu, npcs = 20, verbose = FALSE)
seu <- RunTSNE(seu, reduction = "pca", dims = 1:12)

Confounded dimension reduction

The PCA and t-SNE plots both showed that the data are confounded by species.

The scatterPlot function is used to generate dimension reduction figures. It takes input of a Seurat object (seu here), the name of reduction (pca and tsne here), the variable deciding dot colours (Batch here) and the title of plots. See more information by ?scatterPlot

p1 <- scatterPlot(seu, "pca",colour.by = "Batch", title = "PCA") 
p2 <- scatterPlot(seu, "tsne",colour.by = "Batch", title = "t-SNE") 
plot_grid(p1, p2)

Confounded clustering results

seu <- FindNeighbors(seu, dims = 1:12)
#> Computing nearest neighbor graph
#> Computing SNN
seu <- FindClusters(seu)
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#> 
#> Number of nodes: 10127
#> Number of edges: 352494
#> 
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.9149
#> Number of communities: 22
#> Elapsed time: 1 seconds

scatterPlot(seu, "tsne",colour.by = "seurat_clusters", title = "t-SNE") 

res <- data.frame(table(seu$seurat_clusters, seu$Batch))
ggplot(res, aes(fill=Var2, y=Freq, x=Var1)) + 
    geom_bar(position="stack", stat="identity") + xlab("CIDER_cluster") + ylab("Proportions")

asCIDER

Prepare initial clusters

asCIDER uses existing within-batch clustering results. Here we concatenate the batch ID and the within-batch cluster ID to obtain cluster-specific groups (i.e. initial clusters).

seu$initial_cluster <- paste(seu$Group, seu$Batch, sep = "_")
table(seu$initial_cluster)
#> 
#>             acinar_human activated_stellate_human activated_stellate_mouse 
#>                      932                      275                       14 
#>              alpha_human              alpha_mouse             b_cell_mouse 
#>                     2241                      191                       10 
#>               beta_human               beta_mouse              delta_human 
#>                     2455                      894                      592 
#>              delta_mouse             ductal_human             ductal_mouse 
#>                      218                     1033                      275 
#>        endothelial_human        endothelial_mouse            epsilon_human 
#>                      214                      139                       17 
#>              gamma_human              gamma_mouse       immune_other_mouse 
#>                      241                       41                        8 
#>         macrophage_human         macrophage_mouse               mast_human 
#>                       39                       36                       25 
#> quiescent_stellate_human quiescent_stellate_mouse            schwann_human 
#>                      160                       47                       10 
#>            schwann_mouse             t_cell_human             t_cell_mouse 
#>                        6                        7                        7

Calculate of IDER similarity matrix

The function getIDEr calculate the IDER-based distance matrix.

By default, it will use the column called “initial_cluster” as initial clusters, and “Batch” as batch. If you are using columns other than this two, please revise these two parameters.

For this step, you can choose to use parallel computation by setting use.parallel = TRUE or not (default use.parallel = FALSE). The default number of cores used for parallel computation is detectCores(logical = FALSE) - 1.

ider <- getIDEr(seu, 
                group.by.var = "initial_cluster",
                batch.by.var = "Batch",
                downsampling.size = 35, 
                use.parallel = FALSE, verbose = FALSE)

Visualise the similarity matrix

groups <- c("alpha","beta","delta", "gamma","ductal","endothelial", "activated_stellate", "quiescent_stellate", "macrophage")
idx1 <- paste0(groups, "_human")
idx2 <- paste0(groups, "_mouse")

pheatmap::pheatmap(
  ider[[1]][idx1, idx2],
  color = inferno(10),
  border_color = NA,
  display_numbers = TRUE,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  width = 7,
  height = 5,
  cellwidth = 22,
  cellheight = 22
)

Perform final Clustering

Next we put both Seurat object and the similarity matrix list into finalClustering. You can either cut trees by height (set cutree.by = 'h) or by k (set cutree.by = 'k). The default is cutting by the height of 0.45.

The final clustering results are stored in Colume cider_final_cluster of Seurat object metadata and can be extracted using seu$cider_final_cluster.

seu <- finalClustering(seu, ider, cutree.h = 0.45)
head(seu@meta.data)
#>                             orig.ident nCount_RNA nFeature_RNA Batch   Group
#> mouse1_lib1.final_cell_0001     mouse1      11003         2994 mouse    beta
#> mouse1_lib1.final_cell_0002     mouse1      10701         3565 mouse  ductal
#> mouse1_lib1.final_cell_0003     mouse1       9445         2401 mouse   delta
#> mouse1_lib1.final_cell_0004     mouse1       8193         2865 mouse schwann
#> mouse1_lib1.final_cell_0005     mouse1       7323         2355 mouse   delta
#> mouse1_lib1.final_cell_0006     mouse1       8504         2637 mouse    beta
#>                             Sample RNA_snn_res.0.8 seurat_clusters
#> mouse1_lib1.final_cell_0001 mouse1               9               9
#> mouse1_lib1.final_cell_0002 mouse1              13              13
#> mouse1_lib1.final_cell_0003 mouse1               7               7
#> mouse1_lib1.final_cell_0004 mouse1              21              21
#> mouse1_lib1.final_cell_0005 mouse1               7               7
#> mouse1_lib1.final_cell_0006 mouse1               8               8
#>                             initial_cluster CIDER_cluster
#> mouse1_lib1.final_cell_0001      beta_mouse             1
#> mouse1_lib1.final_cell_0002    ductal_mouse             2
#> mouse1_lib1.final_cell_0003     delta_mouse             3
#> mouse1_lib1.final_cell_0004   schwann_mouse             4
#> mouse1_lib1.final_cell_0005     delta_mouse             3
#> mouse1_lib1.final_cell_0006      beta_mouse             1

Visualise clustering results

plot_list <- list()
plot_list[[1]] <- scatterPlot(seu, "tsne", colour.by = "CIDER_cluster", title = "asCIDER clustering results") 
plot_list[[2]] <- scatterPlot(seu, "tsne", colour.by = "Group", title = "Ground truth of cell populations") 
plot_grid(plotlist = plot_list, ncol = 2)

res <- data.frame(table(seu$CIDER_cluster, seu$Batch))
ggplot(res, aes(fill=Var2, y=Freq, x=Var1)) + 
    geom_bar(position="stack", stat="identity") + xlab("CIDER_cluster") + ylab("Proportions")

Reproducibility

sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Monterey 12.5.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Europe/London
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] viridis_0.6.5      viridisLite_0.4.2  cowplot_1.1.3      ggplot2_3.5.1     
#> [5] pheatmap_1.0.12    Seurat_5.1.0       SeuratObject_5.0.2 sp_2.1-4          
#> [9] CIDER_0.99.4      
#> 
#> loaded via a namespace (and not attached):
#>   [1] RColorBrewer_1.1-3     rstudioapi_0.16.0      jsonlite_1.8.8        
#>   [4] magrittr_2.0.3         spatstat.utils_3.1-0   farver_2.1.2          
#>   [7] rmarkdown_2.27         vctrs_0.6.5            ROCR_1.0-11           
#>  [10] spatstat.explore_3.3-2 htmltools_0.5.8.1      sass_0.4.9            
#>  [13] sctransform_0.4.1      parallelly_1.38.0      KernSmooth_2.23-24    
#>  [16] bslib_0.7.0            htmlwidgets_1.6.4      ica_1.0-3             
#>  [19] plyr_1.8.9             plotly_4.10.4          zoo_1.8-12            
#>  [22] cachem_1.1.0           igraph_2.0.3           mime_0.12             
#>  [25] lifecycle_1.0.4        iterators_1.0.14       pkgconfig_2.0.3       
#>  [28] Matrix_1.7-0           R6_2.5.1               fastmap_1.2.0         
#>  [31] fitdistrplus_1.2-1     future_1.34.0          shiny_1.9.1           
#>  [34] digest_0.6.37          colorspace_2.1-1       patchwork_1.2.0       
#>  [37] tensor_1.5             RSpectra_0.16-2        irlba_2.3.5.1         
#>  [40] labeling_0.4.3         progressr_0.14.0       fansi_1.0.6           
#>  [43] spatstat.sparse_3.1-0  httr_1.4.7             polyclip_1.10-7       
#>  [46] abind_1.4-5            compiler_4.4.1         withr_3.0.1           
#>  [49] doParallel_1.0.17      fastDummies_1.7.4      highr_0.11            
#>  [52] MASS_7.3-61            tools_4.4.1            lmtest_0.9-40         
#>  [55] httpuv_1.6.15          future.apply_1.11.2    goftest_1.2-3         
#>  [58] glue_1.7.0             dbscan_1.2.2           nlme_3.1-165          
#>  [61] promises_1.3.0         grid_4.4.1             Rtsne_0.17            
#>  [64] cluster_2.1.6          reshape2_1.4.4         generics_0.1.3        
#>  [67] gtable_0.3.5           spatstat.data_3.1-2    tidyr_1.3.1           
#>  [70] data.table_1.16.0      utf8_1.2.4             spatstat.geom_3.3-2   
#>  [73] RcppAnnoy_0.0.22       ggrepel_0.9.5          RANN_2.6.2            
#>  [76] foreach_1.5.2          pillar_1.9.0           stringr_1.5.1         
#>  [79] spam_2.10-0            RcppHNSW_0.6.0         limma_3.60.6          
#>  [82] later_1.3.2            splines_4.4.1          dplyr_1.1.4           
#>  [85] lattice_0.22-6         survival_3.7-0         deldir_2.0-4          
#>  [88] tidyselect_1.2.1       locfit_1.5-9.10        miniUI_0.1.1.1        
#>  [91] pbapply_1.7-2          knitr_1.48             gridExtra_2.3         
#>  [94] edgeR_4.2.2            scattermore_1.2        xfun_0.46             
#>  [97] statmod_1.5.0          matrixStats_1.4.1      stringi_1.8.4         
#> [100] lazyeval_0.2.2         yaml_2.3.10            evaluate_0.24.0       
#> [103] codetools_0.2-20       kernlab_0.9-33         tibble_3.2.1          
#> [106] cli_3.6.3              uwot_0.2.2             xtable_1.8-4          
#> [109] reticulate_1.39.0      munsell_0.5.1          jquerylib_0.1.4       
#> [112] Rcpp_1.0.13            globals_0.16.3         spatstat.random_3.3-1 
#> [115] png_0.1-8              spatstat.univar_3.0-1  parallel_4.4.1        
#> [118] dotCall64_1.1-1        listenv_0.9.1          scales_1.3.0          
#> [121] ggridges_0.5.6         leiden_0.4.3.1         purrr_1.0.2           
#> [124] rlang_1.1.4

References

  1. Baron, M. et al. A Single-Cell Transcriptomic Map of the Human and Mouse Pancreas Reveals Inter- and Intra-cell Population Structure. Cell Syst 3, 346–360.e4 (2016).
  2. Satija R, et al. Spatial reconstruction of single-cell gene expression data. Nature Biotechnology 33, 495-502 (2015).
  3. The count matrix and sample information were downloaded from NCBI GEO accession GSE84133.