Doublet prediction using {scDblFinder}

Dependencies

Install Bioconductor packages using BiocManager::install().

if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("scDblFinder")
BiocManager::install("TENxIO")

Load libraries.

suppressPackageStartupMessages(library(scDblFinder))
suppressPackageStartupMessages(library(SingleCellExperiment))
suppressPackageStartupMessages(library(TENxIO))

Data

Download data.

outdir <- 'data/'
filtered_h5_url <- 'https://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v3/pbmc_1k_v3_filtered_feature_bc_matrix.h5'
filtered_h5 <- paste0(outdir, basename(filtered_h5_url))

pbmc5k_1_url <- "https://cf.10xgenomics.com/samples/cell-exp/9.0.0/5k_Human_Donor1_PBMC_3p_gem-x_5k_Human_Donor1_PBMC_3p_gem-x/5k_Human_Donor1_PBMC_3p_gem-x_5k_Human_Donor1_PBMC_3p_gem-x_count_sample_filtered_feature_bc_matrix.h5"
pbmc5k_2_url <- "https://cf.10xgenomics.com/samples/cell-exp/9.0.0/5k_Human_Donor2_PBMC_3p_gem-x_5k_Human_Donor2_PBMC_3p_gem-x/5k_Human_Donor2_PBMC_3p_gem-x_5k_Human_Donor2_PBMC_3p_gem-x_count_sample_filtered_feature_bc_matrix.h5"
pbmc5k_1_h5 <- paste0(outdir, basename(pbmc5k_1_url))
pbmc5k_2_h5 <- paste0(outdir, basename(pbmc5k_2_url))

download_file <- function(url, outfile){
  fn <- basename(url)
  if(file.exists(outfile) == FALSE){
    download.file(url, destfile = outfile)
  } else {
    message(paste0(outfile, " already exists"))
  }
}

download_file(filtered_h5_url, filtered_h5)

data/pbmc_1k_v3_filtered_feature_bc_matrix.h5 already exists

download_file(pbmc5k_1_url, pbmc5k_1_h5)

data/5k_Human_Donor1_PBMC_3p_gem-x_5k_Human_Donor1_PBMC_3p_gem-x_count_sample_filtered_feature_bc_matrix.h5 already exists

download_file(pbmc5k_2_url, pbmc5k_2_h5)

data/5k_Human_Donor2_PBMC_3p_gem-x_5k_Human_Donor2_PBMC_3p_gem-x_count_sample_filtered_feature_bc_matrix.h5 already exists

Create objects

Create SingleCellExperiment files.

create_sce_obj <- function(h5){
  con <- TENxH5(h5)
  import(con)
}

pbmc1k <- create_sce_obj(filtered_h5)

Warning in rhdf5::h5read(file, name = paste0(group, ranges), index = list(1L),
: Object 'matrix/features/interval' does not exist in this HDF5 file.

pbmc1k

class: SingleCellExperiment 
dim: 33538 1222 
metadata(1): TENxFile
assays(1): counts
rownames(33538): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
  ENSG00000268674
rowData names(3): ID Symbol Type
colnames(1222): AAACCCAAGGAGAGTA-1 AAACGCTTCAGCCCAG-1 ...
  TTTGGTTGTAGAATAC-1 TTTGTTGCAATTAGGA-1
colData names(0):
reducedDimNames(0):
mainExpName: Gene Expression
altExpNames(0):

Predict doublets

Run scDblFinder().

set.seed(1984)
scDblFinder(pbmc1k) |>
  colData() |>
  as.data.frame() -> pbmc1k.pred

Creating ~1500 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 163 cells excluded from training.

iter=1, 139 cells excluded from training.

iter=2, 128 cells excluded from training.

Threshold found:0.354

28 (2.3%) doublets called

set.seed(1984)
scDblFinder(pbmc1k) |>
  colData() |>
  as.data.frame() -> pbmc1k.pred2

Creating ~1500 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 163 cells excluded from training.

iter=1, 139 cells excluded from training.

iter=2, 128 cells excluded from training.

Threshold found:0.354

28 (2.3%) doublets called

head(pbmc1k.pred)

                   scDblFinder.class scDblFinder.score scDblFinder.weighted
AAACCCAAGGAGAGTA-1           singlet       0.045929316            0.4312649
AAACGCTTCAGCCCAG-1           singlet       0.004174906            0.2353700
AAAGAACAGACGACTG-1           singlet       0.012409605            0.4142709
AAAGAACCAATGGCAG-1           singlet       0.003421359            0.2900510
AAAGAACGTCTGCAAT-1           singlet       0.008659453            0.2425525
AAAGGATAGTAGACAT-1           singlet       0.016740253            0.1398486
                   scDblFinder.cxds_score
AAACCCAAGGAGAGTA-1            0.005120514
AAACGCTTCAGCCCAG-1            0.065461443
AAAGAACAGACGACTG-1            0.078274898
AAAGAACCAATGGCAG-1            0.056915421
AAAGAACGTCTGCAAT-1            0.114158769
AAAGGATAGTAGACAT-1            0.094354705

head(pbmc1k.pred2)

                   scDblFinder.class scDblFinder.score scDblFinder.weighted
AAACCCAAGGAGAGTA-1           singlet       0.045929316            0.4312649
AAACGCTTCAGCCCAG-1           singlet       0.004174906            0.2353700
AAAGAACAGACGACTG-1           singlet       0.012409605            0.4142709
AAAGAACCAATGGCAG-1           singlet       0.003421359            0.2900510
AAAGAACGTCTGCAAT-1           singlet       0.008659453            0.2425525
AAAGGATAGTAGACAT-1           singlet       0.016740253            0.1398486
                   scDblFinder.cxds_score
AAACCCAAGGAGAGTA-1            0.005120514
AAACGCTTCAGCCCAG-1            0.065461443
AAAGAACAGACGACTG-1            0.078274898
AAAGAACCAATGGCAG-1            0.056915421
AAAGAACGTCTGCAAT-1            0.114158769
AAAGGATAGTAGACAT-1            0.094354705

Plot scores per label.

boxplot(scDblFinder.score ~ scDblFinder.class, data = pbmc1k.pred, pch = 16)

Past versions of boxplot-1.png

Compare runs.

x <- pbmc1k.pred$scDblFinder.score
y <- pbmc1k.pred2$scDblFinder.score

plot(x, y, pch = 16)

Data subset

Since doublets are artifically generated using the available data, providing a different input should affect the scoring and prediction.

set.seed(1984)

n <- ceiling(ncol(pbmc1k)*0.9)
bcs <- sample(x = colnames(pbmc1k), size = n)

pbmc1k_subset <- pbmc1k[, bcs]
pbmc1k_subset

class: SingleCellExperiment 
dim: 33538 1100 
metadata(1): TENxFile
assays(1): counts
rownames(33538): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
  ENSG00000268674
rowData names(3): ID Symbol Type
colnames(1100): AGGACGAAGACCTTTG-1 TTTATGCGTTGATCGT-1 ...
  GTCACTCTCTTGAACG-1 TTCAGGAGTCTACAAC-1
colData names(0):
reducedDimNames(0):
mainExpName: Gene Expression
altExpNames(0):

Predict doublets on the data subset.

pbmc1k_subset.pred <- scDblFinder(pbmc1k_subset)

Creating ~1500 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 136 cells excluded from training.

iter=1, 120 cells excluded from training.

iter=2, 118 cells excluded from training.

Threshold found:0.46

33 (3%) doublets called

colData(pbmc1k_subset.pred) |>
  as.data.frame() -> pbmc1k_subset.pred

idx <- match(
  row.names(pbmc1k_subset.pred),
  row.names(pbmc1k.pred)
)

x <- pbmc1k.pred[idx, ]$scDblFinder.score
y <- pbmc1k_subset.pred$scDblFinder.score

plot(x, y, pch = 16, xlab = 'Full dataset', ylab = 'Subset')

Past versions of pbmc1k_subset_pred-1.png

Running in batches

Create single SingleCellExperiment object.

pbmc5k_1 <- create_sce_obj(pbmc5k_1_h5)

Warning in rhdf5::h5read(file, name = paste0(group, ranges), index = list(1L),
: Object 'matrix/features/interval' does not exist in this HDF5 file.

pbmc5k_2 <- create_sce_obj(pbmc5k_2_h5)

Warning in rhdf5::h5read(file, name = paste0(group, ranges), index = list(1L),
: Object 'matrix/features/interval' does not exist in this HDF5 file.

colnames(pbmc5k_1) <- paste0(colnames(pbmc5k_1), "_1")
colnames(pbmc5k_2) <- paste0(colnames(pbmc5k_2), "_2")

if(any(colnames(pbmc5k_1) %in% colnames(pbmc5k_2)) == TRUE){
  stop("Cell barcodes are not unique")
}

cbind(
  pbmc5k_1,
  pbmc5k_2
) -> pbmc5k

coldata <- DataFrame(
  donor = paste0('donor', sub(pattern = ".*_", replacement = "", x = colnames(pbmc5k)))
)
row.names(coldata) <- colnames(pbmc5k)
colData(pbmc5k) <- coldata

pbmc5k

class: SingleCellExperiment 
dim: 38606 11699 
metadata(2): TENxFile TENxFile
assays(1): counts
rownames(38606): ENSG00000290825 ENSG00000243485 ... ENSG00000278817
  ENSG00000277196
rowData names(3): ID Symbol Type
colnames(11699): AAACCAAAGGTGACGA-1_1 AAACCCTGTGACGAGT-1_1 ...
  TGTGTTAGTTAAGTGG-1_2 TGTGTTGAGCCAATGC-1_2
colData names(1): donor
reducedDimNames(0):
mainExpName: Gene Expression
altExpNames(0):

Predict on combined object but on different samples.

A vector of the same length as cells (or the name of a column of colData(x)), indicating to which sample each cell belongs. Here, a sample is understood as being processed independently. If omitted, doublets will be searched for with all cells together. If given, doublets will be searched for independently for each sample, which is preferable if they represent different captures. If your samples were multiplexed using cell hashes, what you want to give here are the different batches/wells (i.e. independent captures, since doublets cannot arise across them) rather than biological samples.

set.seed(1984)
pbmc5k.pred <- scDblFinder(pbmc5k, samples="donor")

  |                                                                              |                                                                      |   0%  |                                                                              |===================================                                   |  50%  |                                                                              |======================================================================| 100%

Predict on single batch.

set.seed(1984)
pbmc5k_1.pred <- scDblFinder(pbmc5k_1)

Creating ~4568 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 452 cells excluded from training.

iter=1, 456 cells excluded from training.

iter=2, 457 cells excluded from training.

Threshold found:0.337

263 (4.6%) doublets called

Compare results.

idx <- match(colnames(pbmc5k_1.pred), colnames(pbmc5k.pred))

x <- pbmc5k.pred[, idx]$scDblFinder.score
y <- pbmc5k_1.pred$scDblFinder.score

plot(x, y, pch = 16, xlab = "Combined but using batch", ylab = "Independent batch")

Past versions of compare_batch_run-1.png

Numbers.

table(pbmc5k.pred[, idx]$scDblFinder.class)

singlet doublet 
   5462     248 

table(pbmc5k_1.pred$scDblFinder.class)

singlet doublet 
   5447     263 

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=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [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       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] TENxIO_1.8.2                scDblFinder_1.20.2         
 [3] SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
 [5] Biobase_2.66.0              GenomicRanges_1.58.0       
 [7] GenomeInfoDb_1.42.3         IRanges_2.40.1             
 [9] S4Vectors_0.44.0            BiocGenerics_0.52.0        
[11] MatrixGenerics_1.18.1       matrixStats_1.4.1          
[13] workflowr_1.7.1            

loaded via a namespace (and not attached):
  [1] rstudioapi_0.17.1        jsonlite_1.8.9           magrittr_2.0.3          
  [4] ggbeeswarm_0.7.2         rmarkdown_2.28           fs_1.6.4                
  [7] BiocIO_1.16.0            zlibbioc_1.52.0          vctrs_0.6.5             
 [10] Rsamtools_2.22.0         RCurl_1.98-1.16          htmltools_0.5.8.1       
 [13] S4Arrays_1.6.0           BiocBaseUtils_1.8.0      curl_5.2.3              
 [16] BiocNeighbors_2.0.1      Rhdf5lib_1.28.0          xgboost_1.7.8.1         
 [19] SparseArray_1.6.1        rhdf5_2.50.2             sass_0.4.9              
 [22] bslib_0.8.0              cachem_1.1.0             GenomicAlignments_1.42.0
 [25] whisker_0.4.1            igraph_2.1.1             lifecycle_1.0.4         
 [28] pkgconfig_2.0.3          rsvd_1.0.5               Matrix_1.7-0            
 [31] R6_2.5.1                 fastmap_1.2.0            GenomeInfoDbData_1.2.13 
 [34] digest_0.6.37            colorspace_2.1-1         ps_1.8.1                
 [37] rprojroot_2.0.4          scater_1.34.0            dqrng_0.4.1             
 [40] irlba_2.3.5.1            beachmat_2.22.0          fansi_1.0.6             
 [43] httr_1.4.7               abind_1.4-8              compiler_4.4.1          
 [46] BiocParallel_1.40.0      viridis_0.6.5            highr_0.11              
 [49] HDF5Array_1.34.0         R.utils_2.12.3           MASS_7.3-60.2           
 [52] DelayedArray_0.32.0      rjson_0.2.23             bluster_1.16.0          
 [55] tools_4.4.1              vipor_0.4.7              beeswarm_0.4.0          
 [58] httpuv_1.6.15            R.oo_1.26.0              glue_1.8.0              
 [61] restfulr_0.0.15          callr_3.7.6              rhdf5filters_1.18.0     
 [64] promises_1.3.0           grid_4.4.1               getPass_0.2-4           
 [67] cluster_2.1.6            generics_0.1.3           gtable_0.3.6            
 [70] tzdb_0.4.0               R.methodsS3_1.8.2        data.table_1.16.2       
 [73] hms_1.1.3                BiocSingular_1.22.0      ScaledMatrix_1.14.0     
 [76] metapod_1.14.0           utf8_1.2.4               XVector_0.46.0          
 [79] ggrepel_0.9.6            pillar_1.9.0             stringr_1.5.1           
 [82] limma_3.62.2             later_1.3.2              dplyr_1.1.4             
 [85] lattice_0.22-6           rtracklayer_1.66.0       tidyselect_1.2.1        
 [88] locfit_1.5-9.10          Biostrings_2.74.1        scuttle_1.16.0          
 [91] knitr_1.48               git2r_0.35.0             gridExtra_2.3           
 [94] edgeR_4.4.2              xfun_0.48                statmod_1.5.0           
 [97] stringi_1.8.4            UCSC.utils_1.2.0         yaml_2.3.10             
[100] evaluate_1.0.1           codetools_0.2-20         tibble_3.2.1            
[103] cli_3.6.3                munsell_0.5.1            processx_3.8.4          
[106] jquerylib_0.1.4          Rcpp_1.0.13              XML_3.99-0.17           
[109] parallel_4.4.1           ggplot2_3.5.1            readr_2.1.5             
[112] scran_1.34.0             bitops_1.0-9             viridisLite_0.4.2       
[115] scales_1.3.0             crayon_1.5.3             rlang_1.1.4