Doublet origins

Abstract

Functions to identify the origins of doublets, i.e. the cell types composing them, and perform analysis on these origins.

Introduction

This vignette addresses the question of doublet origins, i.e. identifying which combination of cell-types or clusters generated a given doublet. The first part of the vignette concerns the identifying doublet origins, while the second part covers test for enrichment in specific kinds of doublets.

Identifying doublet origins

As a preliminary, note that doublet origins called in versions prior to 1.27.8 were mostly wrong.

kNN-based approach

Identifying doublet origins requires running scDblFinder() on your data passing your cell type labels, i.e. using the clusters argument of scDblFinder() to pass cluster labels (the cluster labels of doublets will be ignored, so no need to worry about those). For example:

# we first generate example data:
sce <- mockDoubletSCE(ncells=c(100,150,200), ngenes=500, dbl.rate=0.15)
# we call doublets (somewhat quickly):
sce <- scDblFinder(sce, artificialDoublets = 300, clusters="cluster", threshold=0.5)

## 3 clusters

## Creating ~300 artificial doublets...

## Dimensional reduction

## Evaluating kNN...

## Training model...

## iter=0, 16 cells excluded from training.

## iter=1, 8 cells excluded from training.

## iter=2, 7 cells excluded from training.

## Threshold found:0.5

## 11 (2.2%) doublets called

Ran in this fashion, the output already includes a rough kNN-based guess of the doublet origin:

table(call=sce$scDblFinder.class,
      mostLikelyOrigin=sce$scDblFinder.mostLikelyOrigin)

##          mostLikelyOrigin
## call      cluster1+cluster2 cluster1+cluster3 cluster2+cluster3
##   singlet                16                15                13
##   doublet                 5                 4                 2

Note that since this is independent of the classifier (and based purely on the kNN), some non-doublets may have a combination of clusters as their mostLikelyOrigin, we recommend ignoring those. (scDblFinder.originAmbiguous in colData also contains information about whether the call is ambiguous or not)

We can compare this to the ground truth (focusing on actual doublets):

isDoublet <- which(sce$type=="doublet")
table(prediction=sce$scDblFinder.mostLikelyOrigin[isDoublet],
      truth=sce$origin[isDoublet])

##                    truth
## prediction          cluster1+cluster2 cluster1+cluster3 cluster2+cluster3
##   cluster1+cluster2                17                 0                 0
##   cluster1+cluster3                 0                14                 0
##   cluster2+cluster3                 0                 0                15

In this easy case, all true doublet origins are accurately identified. However, a slightly more powerful approach is to identify doublets using a specific multi-label classifier.

Multi-label classifier

The identifyDoubletOrigins() function trains a multi-label classifier on doublet origins from artificial doublets, which can then be used to predict doublet origins.

# to keep things quick here, with cap it to 20 rounds: 
clf <- identifyDoubletOrigins(sce, clusters="cluster", max_rounds = 20, verbose = FALSE)

Here the function uses the sce$scDblFinder.class doublets, but one could specify manually for which cells predictions should be made.

We can compare the predictions to the ground truth:

table(prediction=clf$calls, truth=colData(sce)[names(clf$calls), "origin"])

##                    truth
## prediction          cluster1+cluster2 cluster1+cluster3 cluster2+cluster3
##   cluster1+cluster2                 5                 0                 0
##   cluster1+cluster3                 0                 4                 0
##   cluster2+cluster3                 0                 0                 2

Note that in real data, due to the very large variations in library sizes, it is not uncommon for doublets to contain a large fraction of reads from one cell type, and conversely a small fraction from the other cell type. In those circumstances, while one of the two originating cell type is clear, the other can however be mis-assigned.

Doublet enrichment analysis

Doublet enrichment analysis aims to investigate whether certain doublet types are found more often than would be expected by chance. Such enrichment could, for instance, indicate cell-to-cell interactions.

When ran with the clusters argument, scDblFinder() saves some stats on the doublet origins in metadata(sce)$scDblFinder.stats. Similarly, the output of identifyDoubletOrigins() contains a stats slot:

clf$stats

##         combination observed expected deviation prop.deviation        FNR
## 1 cluster1+cluster2        5  0.34584   4.65416       3.599394 0.00000000
## 2 cluster1+cluster3        4  0.41920   3.58080       2.769288 0.01075269
## 3 cluster2+cluster3        2  0.52800   1.47200       1.138403 0.00000000
##   difficulty
## 1 0.05194089
## 2 0.05148621
## 3 0.04395504

This table contains the observed and expected types of doublets, along with the difficulty of their identification. The plotDoubletMap() function can be used to visualize enrichment. Most importantly, these compiled statistics can then be used to test for enrichment.

Specifically, we define two forms of doublet enrichment, each with their own testing function:

The clusterStickiness() function tests whether each cluster forms more doublet than would be expected given its abundance, by default using a single quasi-binomial model fitted across all doublet types. (Only applicable with >=4 clusters)
The doubletPairwiseEnrichment() function separately tests whether each specific doublet type (i.e. combination of clusters) is more abundant than expected, by default using a poisson model.

For example:

doubletPairwiseEnrichment(clf$stats)

## theta=0.188806055483142

##         combination  log2enrich    p.value       FDR
## 1 cluster1+cluster2  0.55183959 0.09357226 0.2807168
## 2 cluster1+cluster3  0.06689129 0.32700281 0.6540056
## 3 cluster2+cluster3 -1.08032686 1.00000000 1.0000000

This shows no significant enrichment in any of the doublet types, which makes sense because they were randomly generated.

For more detail, see the appropriate section of Germain et al., 2022.

Session information

sessionInfo()