Doublet detection with doubletrs

doubletrs detects doublets (droplets containing two cells) in single-cell RNA sequencing data with the DoubletFinder method of McGinnis, Murrow and Gartner (2019). The memory-heavy steps (artificial doublet simulation, nearest-neighbour search and pANN scoring) are written in Rust. DoubletFinder builds a dense cells x cells distance matrix; doubletrs never does.

doubletrs has the same functions and arguments as DoubletFinder, so an existing workflow switches by replacing library(DoubletFinder) with library(doubletrs). Parent cells for artificial doublets are drawn with the same R random number calls, so after the same set.seed() the results are identical to DoubletFinder’s.

A DoubletFinder workflow

The workflow functions need Seurat. This example uses pbmc_small, an 80-cell dataset from SeuratObject that has already been normalised, scaled and run through PCA.

seu <- SeuratObject::pbmc_small

1. Choose pK

paramSweep() computes pANN over a grid of pN and pK values. summarizeSweep() measures how bimodal each pANN distribution is, and find.pK() combines those into BCmvn. The pK with the highest BCmvn is DoubletFinder’s recommended choice.

set.seed(1)
sweep <- paramSweep(seu, PCs = 1:5)
stats <- summarizeSweep(sweep)
bcmvn <- find.pK(stats)

pK <- as.numeric(as.character(bcmvn$pK[which.max(bcmvn$BCmetric)]))
pK
#> [1] 0.24

2. Estimate the number of doublets

The expected number of doublets comes from the loading density of the experiment; here we assume 7.5%. Doublets of two cells of the same type (homotypic) are hard to detect, so the estimate is scaled down by the expected homotypic fraction from cluster labels.

homotypic <- modelHomotypic(seu$RNA_snn_res.1)
nExp <- round(0.075 * ncol(seu) * (1 - homotypic))
nExp
#> [1] 4

3. Call doublets

set.seed(1)
seu <- doubletFinder(seu, PCs = 1:5, pN = 0.25, pK = pK, nExp = nExp)
calls <- grep("^DF.classifications", colnames(seu@meta.data), value = TRUE)
table(seu@meta.data[[calls]])
#> 
#> Doublet Singlet 
#>       4      76

The metadata columns are named as DoubletFinder names them: pANN_<pN>_<pK>_<nExp> and DF.classifications_<pN>_<pK>_<nExp>.

Lower-level functions

The steps are also available without Seurat, for other preprocessing pipelines. simulate_doublets() averages random pairs of cells from a count matrix. After merging real cells and doublets, and normalising and running PCA on the merged matrix in the pipeline of your choice, compute_pann() scores each real cell.

set.seed(1)
counts <- abs(Matrix::rsparsematrix(200, 100, density = 0.2))
doublets <- simulate_doublets(counts, n = 33)
dim(doublets)
#> [1] 200  33

# Stand-in for PCs of the merged matrix: real cells first, then doublets.
pcs <- matrix(rnorm(133 * 5), ncol = 5)
pann <- compute_pann(pcs, n_real = 100, k = 10)
summary(pann)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   0.000   0.100   0.200   0.199   0.300   0.500

find_neighbors() returns the neighbour indices themselves.