scDblFinder 1.12.0
See the relevant section of the OSCA book for an example of the recoverDoublets() function in action on real data.
A toy example is also provided in ?recoverDoublets.
Consider any two cell states \(C_1\) and \(C_2\) forming a doublet population \(D_{12}\). We will focus on the relative frequency of inter-sample to intra-sample doublets in \(D_{12}\). Given a vector \(\vec p_X\) containing the proportion of cells from each sample in state \(X\), and assuming that doublets form randomly between pairs of samples, the expected proportion of intra-sample doublets in \(D_{12}\) is \(\vec p_{C_1} \cdot \vec p_{C_2}\). Subtracting this from 1 gives us the expected proportion of inter-sample doublets \(q_{D_{12}}\). Similarly, the expected proportion of inter-sample doublets in \(C_1\) is just \(q_{C_1} =1 - \| \vec p_{C_1} \|_2^2\).
Now, let’s consider the observed proportion of events \(r_X\) in each state \(X\) that are known doublets. We have \(r_{D_{12}} = q_{D_{12}}\) as there are no other events in \(D_{12}\) beyond actual doublets. On the other hand, we expect that \(r_{C_1} \ll q_{C_1}\) due to presence of a large majority of non-doublet cells in \(C_1\) (same for \(C_2\)). If we assume that \(q_{D_{12}} \ge q_{C_1}\) and \(q_{C_2}\), the observed proportion \(r_{D_{12}}\) should be larger than \(r_{C_1}\) and \(r_{C_2}\). (The last assumption is not always true but the \(\ll\) should give us enough wiggle room to be robust to violations.)
The above reasoning motivates the use of the proportion of known doublet neighbors as a “doublet score” to identify events that are most likely to be themselves doublets.
recoverDoublets() computes the proportion of known doublet neighbors for each cell by performing a \(k\)-nearest neighbor search against all other cells in the dataset.
It is then straightforward to calculate the proportion of neighboring cells that are marked as known doublets, representing our estimate of \(r_X\) for each cell.
While the proportions are informative, there comes a time when we need to convert these into explicit doublet calls.
This is achieved with \(\vec S\), the vector of the proportion of cells from each sample across the entire dataset (i.e., samples).
We assume that all cell states contributing to doublet states have proportion vectors equal to \(\vec S\), such that the expected proportion of doublets that occur between cells from the same sample is \(\| \vec S\|_2^2\).
We then solve
\[ \frac{N_{intra}}{(N_{intra} + N_{inter}} = \| \vec S\|_2^2 \]
for \(N_{intra}\), where \(N_{inter}\) is the number of observed inter-sample doublets. The top \(N_{intra}\) events with the highest scores (and, obviously, are not already inter-sample doublets) are marked as putative intra-sample doublets.
The rate and manner of doublet formation is (mostly) irrelevant as we condition on the number of events in \(D_{12}\). This means that we do not have to make any assumptions about the relative likelihood of doublets forming between pairs of cell types, especially when cell types have different levels of “stickiness” (or worse, stick specifically to certain other cell types). Such convenience is only possible because of the known doublet calls that allow us to focus on the inter- to intra-sample ratio.
The most problematic assumption is that required to obtain \(N_{intra}\) from \(\vec S\). Obtaining a better estimate would require, at least, the knowledge of the two parent states for each doublet population. This can be determined with some simulation-based heuristics but it is likely to be more trouble than it is worth.
In this theoretical framework, we can easily spot a case where our method fails.
If both \(C_1\) and \(C_2\) are unique to a given sample, all events in \(D_{12}\) will be intra-sample doublets.
This means that no events in \(D_{12}\) will ever be detected as inter-sample doublets, which precludes their detection as intra-sample doublets by recoverDoublets.
The computational remedy is to augment the predictions with simulation-based methods (e.g., scDblFinder()) while the experimental remedy is to ensure that multiplexed samples include technical or biological replicates.
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
##
## 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
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] bluster_1.8.0 scDblFinder_1.12.0
## [3] scater_1.26.0 ggplot2_3.3.6
## [5] scran_1.26.0 scuttle_1.8.0
## [7] ensembldb_2.22.0 AnnotationFilter_1.22.0
## [9] GenomicFeatures_1.50.0 AnnotationDbi_1.60.0
## [11] scRNAseq_2.11.0 SingleCellExperiment_1.20.0
## [13] SummarizedExperiment_1.28.0 Biobase_2.58.0
## [15] GenomicRanges_1.50.0 GenomeInfoDb_1.34.0
## [17] IRanges_2.32.0 S4Vectors_0.36.0
## [19] BiocGenerics_0.44.0 MatrixGenerics_1.10.0
## [21] matrixStats_0.62.0 BiocStyle_2.26.0
##
## loaded via a namespace (and not attached):
## [1] AnnotationHub_3.6.0 BiocFileCache_2.6.0
## [3] igraph_1.3.5 lazyeval_0.2.2
## [5] BiocParallel_1.32.0 digest_0.6.30
## [7] htmltools_0.5.3 magick_2.7.3
## [9] viridis_0.6.2 fansi_1.0.3
## [11] magrittr_2.0.3 memoise_2.0.1
## [13] ScaledMatrix_1.6.0 cluster_2.1.4
## [15] limma_3.54.0 Biostrings_2.66.0
## [17] prettyunits_1.1.1 colorspace_2.0-3
## [19] ggrepel_0.9.1 blob_1.2.3
## [21] rappdirs_0.3.3 xfun_0.34
## [23] dplyr_1.0.10 crayon_1.5.2
## [25] RCurl_1.98-1.9 jsonlite_1.8.3
## [27] glue_1.6.2 gtable_0.3.1
## [29] zlibbioc_1.44.0 XVector_0.38.0
## [31] DelayedArray_0.24.0 BiocSingular_1.14.0
## [33] scales_1.2.1 DBI_1.1.3
## [35] edgeR_3.40.0 Rcpp_1.0.9
## [37] viridisLite_0.4.1 xtable_1.8-4
## [39] progress_1.2.2 dqrng_0.3.0
## [41] bit_4.0.4 rsvd_1.0.5
## [43] metapod_1.6.0 httr_1.4.4
## [45] ellipsis_0.3.2 farver_2.1.1
## [47] pkgconfig_2.0.3 XML_3.99-0.12
## [49] sass_0.4.2 dbplyr_2.2.1
## [51] locfit_1.5-9.6 utf8_1.2.2
## [53] labeling_0.4.2 tidyselect_1.2.0
## [55] rlang_1.0.6 later_1.3.0
## [57] munsell_0.5.0 BiocVersion_3.16.0
## [59] tools_4.2.1 cachem_1.0.6
## [61] xgboost_1.6.0.1 cli_3.4.1
## [63] generics_0.1.3 RSQLite_2.2.18
## [65] ExperimentHub_2.6.0 evaluate_0.17
## [67] stringr_1.4.1 fastmap_1.1.0
## [69] yaml_2.3.6 knitr_1.40
## [71] bit64_4.0.5 purrr_0.3.5
## [73] KEGGREST_1.38.0 sparseMatrixStats_1.10.0
## [75] mime_0.12 xml2_1.3.3
## [77] biomaRt_2.54.0 compiler_4.2.1
## [79] beeswarm_0.4.0 filelock_1.0.2
## [81] curl_4.3.3 png_0.1-7
## [83] interactiveDisplayBase_1.36.0 tibble_3.1.8
## [85] statmod_1.4.37 bslib_0.4.0
## [87] stringi_1.7.8 highr_0.9
## [89] lattice_0.20-45 ProtGenerics_1.30.0
## [91] Matrix_1.5-1 vctrs_0.5.0
## [93] pillar_1.8.1 lifecycle_1.0.3
## [95] BiocManager_1.30.19 jquerylib_0.1.4
## [97] BiocNeighbors_1.16.0 cowplot_1.1.1
## [99] data.table_1.14.4 bitops_1.0-7
## [101] irlba_2.3.5.1 httpuv_1.6.6
## [103] rtracklayer_1.58.0 R6_2.5.1
## [105] BiocIO_1.8.0 bookdown_0.29
## [107] promises_1.2.0.1 gridExtra_2.3
## [109] vipor_0.4.5 codetools_0.2-18
## [111] MASS_7.3-58.1 assertthat_0.2.1
## [113] rjson_0.2.21 withr_2.5.0
## [115] GenomicAlignments_1.34.0 Rsamtools_2.14.0
## [117] GenomeInfoDbData_1.2.9 parallel_4.2.1
## [119] hms_1.1.2 grid_4.2.1
## [121] beachmat_2.14.0 rmarkdown_2.17
## [123] DelayedMatrixStats_1.20.0 Rtsne_0.16
## [125] shiny_1.7.3 ggbeeswarm_0.6.0
## [127] restfulr_0.0.15