ccImpute Package Manual
2024-10-29
Contents
1 Introduction
Single-cell RNA sequencing (scRNA-seq) is a powerful technique, but its analysis is hampered by dropout events. These events occur when expressed genes are missed and recorded as zero. This makes distinguishing true zero expression from low expression difficult, affecting downstream analyses like cell type classification. To address this challenge, we introduce ccImpute (Malec, Kurban, and Dalkilic 2022), an R package that leverages consensus clustering. This approach measures cell similarity effectively, allowing ccImpute to identify and impute the most probable dropout events. Compared to existing methods, ccImpute excels in two ways: it delivers superior performance and introduces minimal additional noise, as evidenced by improved clustering on datasets with known cell identities.
1.1 Installation
To install this package, start R (version "4.2") and enter:
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ccImpute")2 Data Pre-Processing
ccImpute is an imputation tool that does not provide functions for
pre-processing the data. This tool expects the user to pre-process the data
before using it. The input data is expected to be in a log-normalized format
and accessible through the SingleCellExperiment object logcounts method. This
manual includes sample minimal pre-processing of a dataset from
scRNAseq database using the
scater tool.
3 Sample Usage
3.1 Required libraries
library(scRNAseq)
library(scater)
library(ccImpute)
library(SingleCellExperiment)
library(stats)
library(mclust)3.2 Input Data
The code below loads the raw mouse neuron dataset from (Usoskin et al. 2015) and performs preprocessing steps to facilitate meaningful analysis, including the computation of log-transformed normalized counts.
sce <- UsoskinBrainData()
X <- cpm(sce)
labels <- colData(sce)$"Level 1"
#Filter bad cells
filt <- !grepl("Empty well", labels) &
!grepl("NF outlier", labels) &
!grepl("TH outlier", labels) &
!grepl("NoN outlier", labels) &
!grepl("NoN", labels) &
!grepl("Central, unsolved", labels) &
!grepl(">1 cell", labels) &
!grepl("Medium", labels)
labels <-labels[filt]
X <- as.matrix(X[,filt])
#Remove genes that are not expressed in any cells:
X <- X[rowSums(X)>0,]
#Recreate the SingleCellExperiment and add log-transformed data:
ann <- data.frame(cell_id = labels)
sce <- SingleCellExperiment(assays = list(normcounts = as.matrix(X)),
colData = ann)
logcounts(sce) <- log(normcounts(sce) + 1)3.3 Pre-processing data
A user may consider performing feature selection before running the imputation. ccImpute only imputes the most probable dropout events and is unlikely to benefit from the presence of scarcely expressed genes nor make any corrections to their expression.
3.4 Adjusted Rand Index (ARI)
Adjusted Rand Index measures the similarity between two data clusterings adjusted for the chance grouping of elements. This measure allows us to evaluate the performance of the clustering algorithm as a similarity to the optimal clustering assignments derived from cell labels.
3.5 Compute Adjusted Rand Index (ARI) without imputation.
# Set seed for reproducibility purposes.
set.seed(0)
# Compute PCA reduction of the dataset
reducedDims(sce) <- list(PCA=prcomp(t(logcounts(sce)))$x)
# Get an actual number of cell types
k <- length(unique(colData(sce)$cell_id))
# Cluster the PCA reduced dataset and store the assignments
set.seed(0)
assgmts <- kmeans(reducedDim(sce, "PCA"), centers = k, iter.max = 1e+09,
nstart = 1000)$cluster
# Use ARI to compare the k-means assignments to label assignments
adjustedRandIndex(assgmts, colData(sce)$cell_id)
#> [1] 0.23190023.6 Perform the imputation with 2 CPU cores and fill in the ‘imputed’ assay.
library(BiocParallel)
BPPARAM = MulticoreParam(2)
sce <- ccImpute(sce, BPPARAM = BPPARAM)
#> Running ccImpute on dataset (622 cells) with 2
#> cores.
#> [Elapsed time: 12.67s] Distance matrix has been computed.
#> [Elapsed time: 12.98s] SVD completed.
#> Warning in runKM(logX, v, maxSets, k, consMin, kmNStart, kmMax, BPPARAM): For potentially better imputation results, please specify the
#> number of clusters (k). Currently estimating k using the
#> Tracy-Widom bound.
#> [Elapsed time: 31.63s] Clustering completed.
#> [Elapsed time: 40.01s] Dropouts identified.
#> [Elapsed time: 40.61s] Dropouts imputed.3.7 Re-compute Adjusted Rand Index (ARI) with imputation.
# Recompute PCA reduction of the dataset
reducedDim(sce, "PCA_imputed") <- prcomp(t(assay(sce, "imputed")))$x
# Cluster the PCA reduced dataset and store the assignments
assgmts <- kmeans(reducedDim(sce, "PCA_imputed"), centers = k,
iter.max = 1e+09, nstart = 1000)$cluster
# Use ARI to compare the k-means assignments to label assignments
adjustedRandIndex(assgmts, colData(sce)$cell_id)
#> [1] 0.89671784 ccImpute Algorithm Overview
ccImpute’s takes the following steps:
(1) Input: ccImpute starts with a log-normalized expression matrix.
(2) Distance calculation: Next, ccImpute first computes the weighted Spearman distance between all the cells in the data.
(3) Dimensionality reduction: This is followed by a single value decomposition (SVD) to reduce the distance matrix to the top \(l\) most informative singular vectors (typically \(l = 0.08 \times min(n,2000)\)).
(4) Clustering: The algorithm then runs multiple instances of the k-means clustering algorithm in parallel (default: eight runs) on different subsets of these singular vectors with the results form to create a consensus matrix.
(5) Dropout identification: Using the consensus matrix and modified expression matrix, ccImpute identifies the most likely dropout events that need imputed.
(6) Imputation: Finally, ccImpute imputes the dropouts using a weighted mean approach. It considers the influence of surrounding values, assigning them weights from the consensus matrix. There are two options for handling dropout values: either they are included in the weighting calculation (Method I-II), or their influence is deliberately skipped (Method III).
5 Key Analytical Choices
In the preceding section, we utilized the ccImpute method by providing only one argument: the SingleCellExperiment object containing the scRNA-seq data. Nevertheless, numerous parameters can be explicitly specified instead of relying on default values. Here, we present the invocation of the ccImpute method, including all the parameters with default values assigned:
ccImpute(sce, dist, nCeil = 2000, svdMaxRatio = 0.08,
maxSets = 8, k, consMin=0.75, kmNStart, kmMax=1000,
fastSolver = TRUE, BPPARAM=bpparam(), verbose = TRUE)This function can be decomposed into a series of steps, providing finer control over the execution of the ccImpute algorithm:
cores <- 2
BPPARAM = MulticoreParam(cores)
w <- rowVars_fast(logcounts(sce), cores)
corMat <- getCorM("spearman", logcounts(sce), w, cores)
v <- doSVD(corMat, svdMaxRatio=.08, nCeil=2000, nCores=cores)
consMtx <- runKM(logX, v, maxSets = 8, k, consMin=0.75, kmNStart, kmMax=1000,
BPPARAM=bpparam())
dropIds <- findDropouts(logX, consMtx)
iLogX <- computeDropouts(consMtx, logX, dropIds,
fastSolver=TRUE, nCores=cores)
assay(sce, "imputed") <- iLogXIn the following sections, we will examine these parameters in more detail and clarify their influence on the imputation performance. For a detailed description of each input argument, please consult the reference manual.
5.1 Distance/Similarity Measures
By default, if the dist parameter is not specified in the ccImpute function, the ccImpute algorithm employs a weighted Spearman correlation measure between cells, with weights corresponding to gene variances. However, any distance or correlation matrix can be utilized for this parameter. Furthermore, the package provides the getCorM function, which can efficiently compute correlations and weighted correlations. Here is an example of using this method to compute Pearson correlation in parallel using 2 cores:
corMat <- getCorM("pearson", logcounts(sce), nCores=2)5.2 Singular Value Decomposition (SVD)
In the singular value decomposition step that follows the computing of the distance matrix, the \(nCeil\) parameter specifies the maximum number of cells used to compute the range of top singular vectors. The number of singular vectors used in clustering runs is in \([.5N \times svdMaxRatio, N \times svdMaxRatio]\). However, with high enough N, the imputation performance drops due to increased noise that is introduced. Experimental data suggest that parameters \(NCeil = 2000\), and \(svdMaxRatio = 0.08\) results in optimal performance. However, different values may work better depending on the distance measure used and how other parameters are modified.
5.3 Clustering: k, kmMaxIter, kmNStart
The clustering step employs the k-means algorithm, which requires choosing the number of clusters (k), the maximum number of iterations (kmMaxIter), and the number of times the algorithm is repeated (kmNStart). Ideally, k should match the actual number of cell types in the dataset. However, overestimating k can still result in good imputation quality. If k is not specified, it is estimated using the Tracy-Widom bound but should be reviewed for correctness.
kmMaxIter sets the maximum iterations for each k-means run and defaults to 1000. kmNStart determines how many times the k-means algorithm is repeated. Repetition allows for different initial cluster selections, which can improve clustering quality. K-means is sensitive to initial centroid choices, typically random data points, and repetition mitigates the risk of suboptimal results. If kmNStart is not set, k-means runs \(1000\) times for \(N <= 2000\) and \(50\) times otherwise.
The bpparam parameter controls parallel execution within this package, leveraging the BiocParallel package’s capabilities. BiocParallel manages parallel execution for k-means clustering and determines the number of cores available for subsequent computations using the openMP library (if ccImpute function is used) or a user-specified value. It’s recommended to set the number of cores parameters to match the number of physical cores on your system for optimal performance. This is an example of threaded parallel computation with 4 cores:
BPPARAM = MulticoreParam(4)
sce <- ccImpute(sce, BPPARAM = BPPARAM)Additionally, using parallel and fast BLAS libraries linked to R can significantly speed up the imputation. The RhpcBLASctl library allows you to control the number of threads used by the BLAS library. For example, to use 4 cores:
library(RhpcBLASctl)
blas_set_num_threads(4)6 Runtime Performance
The latest release of ccImpute demonstrates substantial performance improvements over its initial version, particularly when handling larger datasets. To quantify these enhancements, we measured runtime (RT) in minutes across various combinations of processing cores and dataset sizes (denoted by n).
A heatmap visualization reveals that runtime is indicated by red shades, with lighter shades signifying faster performance. The benchmark was conducted on a system equipped with an Intel Xeon Platinum 8268 24-core CPU, 1.5 TB of RAM, and utilized R 4.3.3 with OpenBlas v0.3.20 (pthread variant).
The results unequivocally demonstrate that the current release of ccImpute significantly outperforms its predecessor, especially with larger datasets. This performance gain is consistent across all tested core configurations, underscoring the effectiveness of the enhancements in optimizing ccImpute for practical, real-world applications.
7 R session information.
#> ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
#> setting value
#> version R Under development (unstable) (2024-10-21 r87258)
#> os Ubuntu 24.04.1 LTS
#> system x86_64, linux-gnu
#> ui X11
#> language (EN)
#> collate C
#> ctype en_US.UTF-8
#> tz America/New_York
#> date 2024-10-29
#> pandoc 3.1.3 @ /usr/bin/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> abind 1.4-8 2024-09-12 [2] CRAN (R 4.5.0)
#> alabaster.base 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> alabaster.matrix 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> alabaster.ranges 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> alabaster.sce 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> alabaster.schemas 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> alabaster.se 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> AnnotationDbi 1.69.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> AnnotationFilter 1.31.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> AnnotationHub 3.15.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> beachmat 2.23.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> beeswarm 0.4.0 2021-06-01 [2] CRAN (R 4.5.0)
#> Biobase * 2.67.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocFileCache 2.15.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocGenerics * 0.53.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocIO 1.17.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocManager 1.30.25 2024-08-28 [2] CRAN (R 4.5.0)
#> BiocNeighbors 2.1.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocParallel * 1.41.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocSingular 1.23.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocStyle * 2.35.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> BiocVersion 3.21.1 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> Biostrings 2.75.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> bit 4.5.0 2024-09-20 [2] CRAN (R 4.5.0)
#> bit64 4.5.2 2024-09-22 [2] CRAN (R 4.5.0)
#> bitops 1.0-9 2024-10-03 [2] CRAN (R 4.5.0)
#> blob 1.2.4 2023-03-17 [2] CRAN (R 4.5.0)
#> bookdown 0.41 2024-10-16 [2] CRAN (R 4.5.0)
#> bslib 0.8.0 2024-07-29 [2] CRAN (R 4.5.0)
#> cachem 1.1.0 2024-05-16 [2] CRAN (R 4.5.0)
#> ccImpute * 1.9.0 2024-10-29 [1] Bioconductor 3.21 (R 4.5.0)
#> cli 3.6.3 2024-06-21 [2] CRAN (R 4.5.0)
#> codetools 0.2-20 2024-03-31 [3] CRAN (R 4.5.0)
#> colorspace 2.1-1 2024-07-26 [2] CRAN (R 4.5.0)
#> crayon 1.5.3 2024-06-20 [2] CRAN (R 4.5.0)
#> curl 5.2.3 2024-09-20 [2] CRAN (R 4.5.0)
#> DBI 1.2.3 2024-06-02 [2] CRAN (R 4.5.0)
#> dbplyr 2.5.0 2024-03-19 [2] CRAN (R 4.5.0)
#> DelayedArray 0.33.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> digest 0.6.37 2024-08-19 [2] CRAN (R 4.5.0)
#> dplyr 1.1.4 2023-11-17 [2] CRAN (R 4.5.0)
#> ensembldb 2.31.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> evaluate 1.0.1 2024-10-10 [2] CRAN (R 4.5.0)
#> ExperimentHub 2.15.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> fansi 1.0.6 2023-12-08 [2] CRAN (R 4.5.0)
#> fastmap 1.2.0 2024-05-15 [2] CRAN (R 4.5.0)
#> filelock 1.0.3 2023-12-11 [2] CRAN (R 4.5.0)
#> generics 0.1.3 2022-07-05 [2] CRAN (R 4.5.0)
#> GenomeInfoDb * 1.43.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> GenomeInfoDbData 1.2.13 2024-10-23 [2] Bioconductor
#> GenomicAlignments 1.43.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> GenomicFeatures 1.59.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> GenomicRanges * 1.59.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> ggbeeswarm 0.7.2 2023-04-29 [2] CRAN (R 4.5.0)
#> ggplot2 * 3.5.1 2024-04-23 [2] CRAN (R 4.5.0)
#> ggrepel 0.9.6 2024-09-07 [2] CRAN (R 4.5.0)
#> glue 1.8.0 2024-09-30 [2] CRAN (R 4.5.0)
#> gridExtra 2.3 2017-09-09 [2] CRAN (R 4.5.0)
#> gtable 0.3.6 2024-10-25 [2] CRAN (R 4.5.0)
#> gypsum 1.3.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> HDF5Array 1.35.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> htmltools 0.5.8.1 2024-04-04 [2] CRAN (R 4.5.0)
#> httr 1.4.7 2023-08-15 [2] CRAN (R 4.5.0)
#> httr2 1.0.5 2024-09-26 [2] CRAN (R 4.5.0)
#> IRanges * 2.41.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> irlba 2.3.5.1 2022-10-03 [2] CRAN (R 4.5.0)
#> jquerylib 0.1.4 2021-04-26 [2] CRAN (R 4.5.0)
#> jsonlite 1.8.9 2024-09-20 [2] CRAN (R 4.5.0)
#> KEGGREST 1.47.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> knitr 1.48 2024-07-07 [2] CRAN (R 4.5.0)
#> lattice 0.22-6 2024-03-20 [3] CRAN (R 4.5.0)
#> lazyeval 0.2.2 2019-03-15 [2] CRAN (R 4.5.0)
#> lifecycle 1.0.4 2023-11-07 [2] CRAN (R 4.5.0)
#> magrittr 2.0.3 2022-03-30 [2] CRAN (R 4.5.0)
#> Matrix 1.7-1 2024-10-18 [3] CRAN (R 4.5.0)
#> MatrixGenerics * 1.19.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> matrixStats * 1.4.1 2024-09-08 [2] CRAN (R 4.5.0)
#> mclust * 6.1.1 2024-04-29 [2] CRAN (R 4.5.0)
#> memoise 2.0.1 2021-11-26 [2] CRAN (R 4.5.0)
#> munsell 0.5.1 2024-04-01 [2] CRAN (R 4.5.0)
#> pillar 1.9.0 2023-03-22 [2] CRAN (R 4.5.0)
#> pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.5.0)
#> png 0.1-8 2022-11-29 [2] CRAN (R 4.5.0)
#> ProtGenerics 1.39.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> R6 2.5.1 2021-08-19 [2] CRAN (R 4.5.0)
#> rappdirs 0.3.3 2021-01-31 [2] CRAN (R 4.5.0)
#> Rcpp 1.0.13 2024-07-17 [2] CRAN (R 4.5.0)
#> RCurl 1.98-1.16 2024-07-11 [2] CRAN (R 4.5.0)
#> restfulr 0.0.15 2022-06-16 [2] CRAN (R 4.5.0)
#> rhdf5 2.51.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> rhdf5filters 1.19.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> Rhdf5lib 1.29.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> rjson 0.2.23 2024-09-16 [2] CRAN (R 4.5.0)
#> rlang 1.1.4 2024-06-04 [2] CRAN (R 4.5.0)
#> rmarkdown 2.28 2024-08-17 [2] CRAN (R 4.5.0)
#> Rsamtools 2.23.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> RSQLite 2.3.7 2024-05-27 [2] CRAN (R 4.5.0)
#> rsvd 1.0.5 2021-04-16 [2] CRAN (R 4.5.0)
#> rtracklayer 1.67.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> S4Arrays 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> S4Vectors * 0.45.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> sass 0.4.9 2024-03-15 [2] CRAN (R 4.5.0)
#> ScaledMatrix 1.15.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> scales 1.3.0 2023-11-28 [2] CRAN (R 4.5.0)
#> scater * 1.35.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> scRNAseq * 2.19.1 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> scuttle * 1.17.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> sessioninfo * 1.2.2 2021-12-06 [2] CRAN (R 4.5.0)
#> SingleCellExperiment * 1.29.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> SparseArray 1.7.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> sparseMatrixStats 1.19.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> SummarizedExperiment * 1.37.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> tibble 3.2.1 2023-03-20 [2] CRAN (R 4.5.0)
#> tidyselect 1.2.1 2024-03-11 [2] CRAN (R 4.5.0)
#> UCSC.utils 1.3.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> utf8 1.2.4 2023-10-22 [2] CRAN (R 4.5.0)
#> vctrs 0.6.5 2023-12-01 [2] CRAN (R 4.5.0)
#> vipor 0.4.7 2023-12-18 [2] CRAN (R 4.5.0)
#> viridis 0.6.5 2024-01-29 [2] CRAN (R 4.5.0)
#> viridisLite 0.4.2 2023-05-02 [2] CRAN (R 4.5.0)
#> withr 3.0.2 2024-10-28 [2] CRAN (R 4.5.0)
#> xfun 0.48 2024-10-03 [2] CRAN (R 4.5.0)
#> XML 3.99-0.17 2024-06-25 [2] CRAN (R 4.5.0)
#> XVector 0.47.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#> yaml 2.3.10 2024-07-26 [2] CRAN (R 4.5.0)
#> zlibbioc 1.53.0 2024-10-29 [2] Bioconductor 3.21 (R 4.5.0)
#>
#> [1] /tmp/Rtmp4Fek0s/Rinst8439051d269d4
#> [2] /home/biocbuild/bbs-3.21-bioc/R/site-library
#> [3] /home/biocbuild/bbs-3.21-bioc/R/library
#>
#> ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────