Skip to contents

tl;dr

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.

Mathematical background

Consider any two cell states C1C_1 and C2C_2 forming a doublet population D12D_{12}. We will focus on the relative frequency of inter-sample to intra-sample doublets in D12D_{12}. Given a vector pX\vec p_X containing the proportion of cells from each sample in state XX, and assuming that doublets form randomly between pairs of samples, the expected proportion of intra-sample doublets in D12D_{12} is pC1pC2\vec p_{C_1} \cdot \vec p_{C_2}. Subtracting this from 1 gives us the expected proportion of inter-sample doublets qD12q_{D_{12}}. Similarly, the expected proportion of inter-sample doublets in C1C_1 is just qC1=1pC122q_{C_1} =1 - \| \vec p_{C_1} \|_2^2.

Now, let’s consider the observed proportion of events rXr_X in each state XX that are known doublets. We have rD12=qD12r_{D_{12}} = q_{D_{12}} as there are no other events in D12D_{12} beyond actual doublets. On the other hand, we expect that rC1qC1r_{C_1} \ll q_{C_1} due to presence of a large majority of non-doublet cells in C1C_1 (same for C2C_2). If we assume that qD12qC1q_{D_{12}} \ge q_{C_1} and qC2q_{C_2}, the observed proportion rD12r_{D_{12}} should be larger than rC1r_{C_1} and rC2r_{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 kk-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 rXr_X for each cell.

Obtaining explicit calls

While the proportions are informative, there comes a time when we need to convert these into explicit doublet calls. This is achieved with S\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 S\vec S, such that the expected proportion of doublets that occur between cells from the same sample is S22\| \vec S\|_2^2. We then solve

Nintra(Nintra+Ninter=S22 \frac{N_{intra}}{(N_{intra} + N_{inter}} = \| \vec S\|_2^2

for NintraN_{intra}, where NinterN_{inter} is the number of observed inter-sample doublets. The top NintraN_{intra} events with the highest scores (and, obviously, are not already inter-sample doublets) are marked as putative intra-sample doublets.

Discussion

The rate and manner of doublet formation is (mostly) irrelevant as we condition on the number of events in D12D_{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 NintraN_{intra} from S\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 C1C_1 and C2C_2 are unique to a given sample, all events in D12D_{12} will be intra-sample doublets. This means that no events in D12D_{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.

Session information

## 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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] BiocStyle_2.32.1
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.37       desc_1.4.3          R6_2.5.1           
##  [4] bookdown_0.40       fastmap_1.2.0       xfun_0.47          
##  [7] cachem_1.1.0        knitr_1.48          htmltools_0.5.8.1  
## [10] rmarkdown_2.28      lifecycle_1.0.4     cli_3.6.3          
## [13] sass_0.4.9          pkgdown_2.1.1       textshaping_0.4.0  
## [16] jquerylib_0.1.4     systemfonts_1.1.0   compiler_4.4.1     
## [19] tools_4.4.1         ragg_1.3.3          bslib_0.8.0        
## [22] evaluate_1.0.0      yaml_2.3.10         BiocManager_1.30.25
## [25] jsonlite_1.8.8      rlang_1.1.4         fs_1.6.4