Scoring potential doublets from simulated densities

tl;dr

To demonstrate, we’ll use one of the mammary gland datasets from the scRNAseq package. We will subset it down to a random set of 1000 cells for speed.

library(scRNAseq)
sce <- BachMammaryData(samples="G_1")

set.seed(1001)
sce <- sce[,sample(ncol(sce), 1000)]

For the purposes of this demonstration, we’ll perform an extremely expedited analysis. One would usually take more care here and do some quality control, create some diagnostic plots, etc., but we don’t have the space for that.

library(scuttle)
sce <- logNormCounts(sce)

library(scran)
dec <- modelGeneVar(sce)
hvgs <- getTopHVGs(dec, n=1000)

library(scater)
set.seed(1002)
sce <- runPCA(sce, ncomponents=10, subset_row=hvgs)
sce <- runTSNE(sce, dimred="PCA")

We run computeDoubletDensity() to obtain a doublet score for each cell based on the density of simulated doublets around it. We log this to get some better dynamic range.

set.seed(1003)
library(scDblFinder)
scores <- computeDoubletDensity(sce, subset.row=hvgs)
plotTSNE(sce, colour_by=I(log1p(scores)))

Algorithm overview

We use a fairly simple approach in doubletCells that involves creating simulated doublets from the original data set:

Perform a PCA on the log-normalized expression for all cells in the dataset.
Randomly select two cells and add their count profiles together. Compute the log-normalized profile and project it into the PC space.
Repeat 2 to obtain $N_s$ simulated doublet cells.
For each cell, compute the local density of simulated doublets, scaled by the density of the original cells. This is used as the doublet score.

Size factor handling

Normalization size factors

We allow specification of two sets of size factors for different purposes. The first set is the normalization set: division of counts by these size factors yields expression values to be compared across cells. This is necessary to compute log-normalized expression values for the PCA.

These size factors are usually computed from some method that assumes most genes are not DE. We default to library size normalization though any arbitrary set of size factors can be used. The size factor for each doublet is computed as the sum of size factors for the individual cells, based on the additivity of scaling biases.

RNA content size factors

The second set is the RNA content set: division of counts by these size factors yields expression values that are proportional to absolute abundance across cells. This affects the creation of simulated doublets by controlling the scaling of the count profiles for the individual cells. These size factors would normally be estimated with spike-ins, but in their absence we default to using unity for all cells.

The use of unity values implies that the library size for each cell is a good proxy for total RNA content. This is unlikely to be true: technical biases mean that the library size is an imprecise relative estimate of the content. Saturation effects and composition biases also mean that the expected library size for each population is not an accurate estimate of content. The imprecision will spread out the simulated doublets while the inaccuracy will result in a systematic shift from the location of true doublets.

Arguably, such problems exist for any doublet estimation method without spike-in information. We can only hope that the inaccuracies have only minor effects on the creation of simulated cells. Indeed, the first effect does mitigate the second to some extent by ensuring that some simulated doublets will occupy the neighbourhood of the true doublets.

Interactions between them

These two sets of size factors play different roles so it is possible to specify both of them. We use the following algorithm to accommodate non-unity values for the RNA content size factors:

The RNA content size factors are used to scale the counts first. This ensures that RNA content has the desired effect in step 2 of Section @ref(overview).
The normalization size factors are also divided by the content size factors. This ensures that normalization has the correct effect, see below.
The rest of the algorithm proceeds as if the RNA content size factors were unity. Addition of count profiles is done without further scaling, and normalized expression values are computed with the rescaled normalization size factors.

To understand the correctness of the rescaled normalization size factors, consider a non-DE gene with abundance $\lambda_g$ . The expected count in each cell is $\lambda_g s_i$ for scaling bias $s_i$ (i.e., normalization size factor). The rescaled count is $\lambda_g s_i c_i^{-1}$ for some RNA content size factor $c_i$ . The rescaled normalization size factor is $s_i c_i^{-1}$ , such that normalization yields $\lambda_g$ as desired. This also holds for doublets where the scaling biases and size factors are additive.

Doublet score calculations

We assume that the simulation accurately mimics doublet creation - amongst other things, we assume that doublets are equally likely to form between any cell populations and any differences in total RNA between subpopulations are captured or negligible. If these assumptions hold, then at any given region in the expression space, the number of doublets among the real cells is proportional to the number of simulated doublets lying in the same region. Thus, the probability that a cell is a doublet is proportional to the ratio of the number of neighboring simulated doublets to the number of neighboring real cells.

A mild additional challenge here is that the number of simulated cells $N_s$ can vary. Ideally, we would like the expected output of the function to be the same regardless of the user’s choice of $N_s$ , i.e., the chosen value should only affect the precision/speed trade-off. Many other doublet-based methods take a $k$ -nearest neighbours approach to compute densities; but if $N_s$ is too large relative to the number of real cells, all of the $k$ nearest neighbours will be simulated, while if $N_s$ is too small, all of the nearest neighbors will be original cells.

Thus, we use a modified version of the $k$ NN approach whereby we identify the distance from each cell to its $k$ -th nearest neighbor. This defines a hypersphere around that cell in which we count the number of simulated cells. We then compute the odds ratio of the number of simulated cells in the hypersphere to $N_s$ , divided by the ratio of $k$ to the total number of cells in the dataset. This score captures the relative frequency of simulated cells to real cells while being robust to changes to $N_s$ .

Session information

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] bluster_1.14.0              scDblFinder_1.19.6         
##  [3] scater_1.32.1               ggplot2_3.5.1              
##  [5] scran_1.32.0                scuttle_1.14.0             
##  [7] ensembldb_2.28.1            AnnotationFilter_1.28.0    
##  [9] GenomicFeatures_1.56.0      AnnotationDbi_1.66.0       
## [11] scRNAseq_2.18.0             SingleCellExperiment_1.26.0
## [13] SummarizedExperiment_1.34.0 Biobase_2.64.0             
## [15] GenomicRanges_1.56.1        GenomeInfoDb_1.40.1        
## [17] IRanges_2.38.1              S4Vectors_0.42.1           
## [19] BiocGenerics_0.50.0         MatrixGenerics_1.16.0      
## [21] matrixStats_1.4.1           BiocStyle_2.32.1           
## 
## loaded via a namespace (and not attached):
##   [1] jsonlite_1.8.8            magrittr_2.0.3           
##   [3] ggbeeswarm_0.7.2          gypsum_1.0.1             
##   [5] farver_2.1.2              rmarkdown_2.28           
##   [7] fs_1.6.4                  BiocIO_1.14.0            
##   [9] zlibbioc_1.50.0           ragg_1.3.3               
##  [11] vctrs_0.6.5               memoise_2.0.1            
##  [13] Rsamtools_2.20.0          DelayedMatrixStats_1.26.0
##  [15] RCurl_1.98-1.16           htmltools_0.5.8.1        
##  [17] S4Arrays_1.4.1            AnnotationHub_3.12.0     
##  [19] curl_5.2.2                BiocNeighbors_1.22.0     
##  [21] xgboost_1.7.8.1           Rhdf5lib_1.26.0          
##  [23] SparseArray_1.4.8         rhdf5_2.48.0             
##  [25] sass_0.4.9                alabaster.base_1.4.2     
##  [27] bslib_0.8.0               desc_1.4.3               
##  [29] alabaster.sce_1.4.0       httr2_1.0.4              
##  [31] cachem_1.1.0              GenomicAlignments_1.40.0 
##  [33] igraph_2.0.3              mime_0.12                
##  [35] lifecycle_1.0.4           pkgconfig_2.0.3          
##  [37] rsvd_1.0.5                Matrix_1.7-0             
##  [39] R6_2.5.1                  fastmap_1.2.0            
##  [41] GenomeInfoDbData_1.2.12   digest_0.6.37            
##  [43] colorspace_2.1-1          dqrng_0.4.1              
##  [45] irlba_2.3.5.1             ExperimentHub_2.12.0     
##  [47] textshaping_0.4.0         RSQLite_2.3.7            
##  [49] beachmat_2.20.0           labeling_0.4.3           
##  [51] filelock_1.0.3            fansi_1.0.6              
##  [53] httr_1.4.7                abind_1.4-8              
##  [55] compiler_4.4.1            bit64_4.0.5              
##  [57] withr_3.0.1               BiocParallel_1.38.0      
##  [59] viridis_0.6.5             DBI_1.2.3                
##  [61] highr_0.11                HDF5Array_1.32.1         
##  [63] alabaster.ranges_1.4.2    alabaster.schemas_1.4.0  
##  [65] MASS_7.3-60.2             rappdirs_0.3.3           
##  [67] DelayedArray_0.30.1       rjson_0.2.23             
##  [69] tools_4.4.1               vipor_0.4.7              
##  [71] beeswarm_0.4.0            glue_1.7.0               
##  [73] restfulr_0.0.15           rhdf5filters_1.16.0      
##  [75] grid_4.4.1                Rtsne_0.17               
##  [77] cluster_2.1.6             generics_0.1.3           
##  [79] gtable_0.3.5              data.table_1.16.0        
##  [81] metapod_1.12.0            BiocSingular_1.20.0      
##  [83] ScaledMatrix_1.12.0       utf8_1.2.4               
##  [85] XVector_0.44.0            ggrepel_0.9.6            
##  [87] BiocVersion_3.19.1        pillar_1.9.0             
##  [89] limma_3.60.4              dplyr_1.1.4              
##  [91] BiocFileCache_2.12.0      lattice_0.22-6           
##  [93] rtracklayer_1.64.0        bit_4.0.5                
##  [95] tidyselect_1.2.1          locfit_1.5-9.10          
##  [97] Biostrings_2.72.1         knitr_1.48               
##  [99] gridExtra_2.3             bookdown_0.40            
## [101] ProtGenerics_1.36.0       edgeR_4.2.1              
## [103] xfun_0.47                 statmod_1.5.0            
## [105] UCSC.utils_1.0.0          lazyeval_0.2.2           
## [107] yaml_2.3.10               evaluate_1.0.0           
## [109] codetools_0.2-20          tibble_3.2.1             
## [111] alabaster.matrix_1.4.2    BiocManager_1.30.25      
## [113] cli_3.6.3                 systemfonts_1.1.0        
## [115] munsell_0.5.1             jquerylib_0.1.4          
## [117] Rcpp_1.0.13               dbplyr_2.5.0             
## [119] png_0.1-8                 XML_3.99-0.17            
## [121] parallel_4.4.1            pkgdown_2.1.1            
## [123] blob_1.2.4                sparseMatrixStats_1.16.0 
## [125] bitops_1.0-8              viridisLite_0.4.2        
## [127] alabaster.se_1.4.1        scales_1.3.0             
## [129] purrr_1.0.2               crayon_1.5.3             
## [131] rlang_1.1.4               KEGGREST_1.44.1

Aaron Lun

2024-09-19