An introduction to the methods implemented for doublet detection in single-cell ATAC-seq.
scDblFinder 1.14.0
Analyses in single-cell RNAseq are typically limited to a relatively small (e.g. one or two thousands) set of features that are most informative; these are often the genes with a higher expression (and hence more chances of being quantified). In contrast, single-cell ATACseq (scATACseq) data is considerably more sparse, with most reads being spread across hundreds of thousands of regions. In this context, selecting a subset of genes is highly ineffective, and therefore many of the methods developed for single-cell RNAseq are not easily applicable, and need to be adapted. Methods have therefore been developed specifically for scATACseq data (Granja et al. 2021; Thibodeau et al. 2021).
This vignette presents different approaches to doublet detection in single-cell ATAC-seq implemented in the package: the first is an adaptation of scDblFinder, the second a reimplementation of the AMULET method from Thibodeau et al. (2021). The latter has the advantage of capturing homotypic doublets, but does not perform well in all datasets, and especially requires the cells to have a high library size. We therefore next present two ways of combining the two.
With default parameters, the scDblFinder method performs very poorly on scATACseq data due to the increase spread of the reads across many features. Since working with all features (i.e. tiles or peaks) is computationally very expensive, an alternative approach is to begin by reducing the size of the dataset, not through selection (as in scRNAseq), but by aggregating correlated features into a relatively small set.
This has the advantage of using all information, as well as making the count data more continuous.
This method yields comparable performance to specialized single-cell ATACseq software (Germain et al., 2021).
The feature aggregation can be triggered using the aggregateFeatures=TRUE argument, which will aggregate peak or tile counts into the number of meta-features defined by the nfeatures.
If the number of meta-features is low (which we recommend), the meta-features can be directly used to calculated distances rather than going through the SVD step (which can be triggered with the processing argument). Such an example would be:
library(scDblFinder)
# we use a dummy SingleCellExperiment as example:
sce <- mockDoubletSCE(ngenes=300)
# setting low number of artificial doublets (same as ncells) just for speedup:
sce <- scDblFinder(sce, artificialDoublets=1, aggregateFeatures=TRUE, nfeatures=25, processing="normFeatures")## Aggregating features...## Creating ~526 artificial doublets...## Dimensional reduction## Evaluating kNN...## Training model...## iter=0, 23 cells excluded from training.## iter=1, 28 cells excluded from training.## iter=2, 33 cells excluded from training.## Threshold found:0.562## 22 (4.2%) doublets calledThe AMULET method from Thibodeau et al. (2021) is based on the assumption that, in a diploid cell, any given genomic region should be captured at most twice. Therefore, cells with loci covered by more than two fragments are indicative of the droplet being a doublet. Of note, this approach has the advantage of capturing homotypic doublets, which instead tend to be missed by other methods. Since it was only available in the form of a mixture of java and python scripts, we re-implemented the method in scDblFinder (see ?amulet), leading to equal or superior results to the original implementation (Germain et al. 2021).
As in the original implementation, we recommend excluding the mitochondrial and sex chromosomes, as well as repetitive regions. This can be specified with the regionsToExclude argument (see the underlying ?getFragmentOverlaps). It can be used as follows:
# here we use a dummy fragment file for example:
fragfile <- system.file("extdata", "example_fragments.tsv.gz", package="scDblFinder")
# we might also give a GRanges of repeat elements, so that these regions are excluded:
suppressPackageStartupMessages(library(GenomicRanges))
repeats <- GRanges("chr6", IRanges(1000,2000))
# it's better to combine these with mitochondrial and sex chromosomes
otherChroms <- GRanges(c("M","chrM","MT","X","Y","chrX","chrY"),IRanges(1L,width=10^8))
# here since I don't know what chromosome notation you'll be using I've just put them all,
# although this will trigger a warning when combining them:
toExclude <- suppressWarnings(c(repeats, otherChroms))
# we then launch the method
res <- amulet(fragfile, regionsToExclude=toExclude)## Fragment file is not tabix-indexed, requiring thewhole file to be imported in memory.## 18:39:33 - Splitting and subsetting barcodes...## 18:39:33 - Obtaining overlaps...res## nFrags uniqFrags nAbove2 total.nAbove2 p.value q.value
## barcode1 878 878 1 1 0.475069053 0.791781755
## barcode2 2401 2401 0 0 0.798103482 0.798103482
## barcode3 2325 2325 1 1 0.475069053 0.791781755
## barcode4 1882 1882 0 0 0.798103482 0.798103482
## barcode5 1355 1355 6 6 0.001335761 0.006678806The results is a data.frame with statistics for each barcode, including a p-value. In contrast to the scDblFinder score, a lower p-value here is indicative of the droplet being more likely to be a doublet (as in the original method).
By default, only the barcodes with a minimum number of reads are considered, but it is possible to specify the droplets for which to gather statistics using the barcodes argument.
While the package includes an implementation that works based on peak/tile count matrices (see ?amuletFromCounts), it has a much lower performance with respect to the one based directly on the fragment files (see ?amulet), and we therefore discourage its use.
The workhorse behind the amulet function is the getFragmentOverlaps, which also includes all of the relevant arguments.
If the fragment files are not Tabix-indexed, the whole fragment file will have to be loaded in memory for processing; while this ensures relatively rapid computation, it has high memory requirements. Therefore, if the fragment file is Tabix-indexed (as is for instance done as part of the ArchR pipeline), it will be read and processed per chromosome, which is a little slower due to overhead, but keeps memory requirements rather low. This behavior can be disabled by specifying fullInMemory=TRUE.
While the scDblFinder-based approach generally performs well, none of the two approach is optimal across all datasets tested. We therefore investigated two strategies for combining the rationales of each approach.
The Amulet method tends to perform best with datasets that have homotypic doublets and where cells have a high library size (i.e. median library size per cell of 10-15k reads), while the scDblFinder-based approach works better for heterotypic doublets. Until an optimal solution is found, we recommend using multiple approaches to inform decisions, in particular using the p-value combination method below.
The clamulet method (Classification-powered Amulet-like method) operates similarly to the scDblFinder method, but generates artificial doublets by operating on the fragment coverages. This has the advantage that the number of loci covered by more than two reads can be computed for artificial doublets, enabling the use of this feature (along with the kNN-based ones) in a classification scheme. It however has the disadvantage of being rather slow and memory hungry, and appears to be outperformed by a simple p-value combination of the two methods (see below).
The clamulet method uses the aforementioned aggregation approach, and its usage includes a number of arguments from both the scDblFinder and amulet method (see in particular ?getFragmentOverlaps):
# not run
d <- clamulet("path/to/fragments.tsv.gz")Since our dummy fragment file is so small (5 barcodes), here we’ll have to adjust the arguments for an example to run:
d <- clamulet(fragfile, k=2, nfeatures=3)## 18:39:35 - Reading full fragments...## 18:39:35 - Splitting and subsetting barcodes...## 18:39:35 - Computing coverages## 18:39:35 - Obtaining windows## 18:39:35 - Obtaining window counts## 18:39:35 - Aggregating features## Warning in (function (A, nv = 5, nu = nv, maxit = 1000, work = nv + 7, reorth =
## TRUE, : You're computing too large a percentage of total singular values, use a
## standard svd instead.## 18:39:35 - Computing features for artificial doublets## 18:39:35 - Counting overlaps for real cells## 18:39:35 - Counting overlaps for artificial doublets## 18:39:35 - Scoring network## 18:39:35 - Iterative training## iter=0, 0 cells excluded from training.## iter=1, 0 cells excluded from training.## 18:39:36 Done!d## total nAbove2 total.nAbove2 weighted ratio.k2 include.in.training
## barcode1 23 1 1 0.08405864 0.5 TRUE
## barcode2 17 0 0 NaN 1.0 TRUE
## barcode3 9 1 1 1.00000000 1.0 TRUE
## barcode4 12 0 0 NaN 1.0 TRUE
## barcode5 21 6 6 0.09141137 0.5 TRUE
## score
## barcode1 0.5491791
## barcode2 0.3417997
## barcode3 0.3417997
## barcode4 0.3417997
## barcode5 0.3417997The score can then be interpreted as for scDblFinder. We however note that this method proved inferior to alternatives.
The amulet and scDblFinder scores above can be simply combined by treating them as p-values and aggregating them (here using Fisher’s method from the aggregation package, but see also the metap package):
res$scDblFinder.p <- 1-colData(sce)[row.names(res), "scDblFinder.score"]
res$combined <- apply(res[,c("scDblFinder.p", "p.value")], 1, FUN=function(x){
x[x<0.001] <- 0.001 # prevent too much skew from very small or 0 p-values
suppressWarnings(aggregation::fisher(x))
})We found this to perform better than averaging the scores or their ranks, and while it is not the very best method in any of the datasets tested, it has a more robust performance overall (see Germain et al., 2021).
Jeffrey M. Granja et al., “ArchR Is a Scalable Software Package for Integrative Single-Cell Chromatin Accessibility Analysis,” Nature Genetics, February 25, 2021, 1–9, https://doi.org/10.1038/s41588-021-00790-6
Asa Thibodeau et al., “AMULET: A Novel Read Count-Based Method for Effective Multiplet Detection from Single Nucleus ATAC-Seq Data,” Genome Biology 22, no. 1 (December 2021): 252, https://doi.org/10.1186/s13059-021-02469-x
Pierre-Luc Germain et al., “Doublet Identification in Single-Cell Sequencing Data Using ScDblFinder” (F1000Research, September 28, 2021), https://doi.org/10.12688/f1000research.73600.1
sessionInfo()## R version 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## 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
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] bluster_1.10.0 scDblFinder_1.14.0
## [3] scater_1.28.0 ggplot2_3.4.2
## [5] scran_1.28.0 scuttle_1.10.0
## [7] ensembldb_2.24.0 AnnotationFilter_1.24.0
## [9] GenomicFeatures_1.52.0 AnnotationDbi_1.62.0
## [11] scRNAseq_2.13.0 SingleCellExperiment_1.22.0
## [13] SummarizedExperiment_1.30.0 Biobase_2.60.0
## [15] GenomicRanges_1.52.0 GenomeInfoDb_1.36.0
## [17] IRanges_2.34.0 S4Vectors_0.38.0
## [19] BiocGenerics_0.46.0 MatrixGenerics_1.12.0
## [21] matrixStats_0.63.0 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.4 magrittr_2.0.3
## [3] magick_2.7.4 ggbeeswarm_0.7.1
## [5] farver_2.1.1 rmarkdown_2.21
## [7] BiocIO_1.10.0 zlibbioc_1.46.0
## [9] vctrs_0.6.2 memoise_2.0.1
## [11] Rsamtools_2.16.0 DelayedMatrixStats_1.22.0
## [13] RCurl_1.98-1.12 htmltools_0.5.5
## [15] progress_1.2.2 AnnotationHub_3.8.0
## [17] curl_5.0.0 BiocNeighbors_1.18.0
## [19] xgboost_1.7.5.1 sass_0.4.5
## [21] bslib_0.4.2 cachem_1.0.7
## [23] GenomicAlignments_1.36.0 igraph_1.4.2
## [25] mime_0.12 lifecycle_1.0.3
## [27] pkgconfig_2.0.3 rsvd_1.0.5
## [29] Matrix_1.5-4 R6_2.5.1
## [31] fastmap_1.1.1 GenomeInfoDbData_1.2.10
## [33] shiny_1.7.4 digest_0.6.31
## [35] colorspace_2.1-0 dqrng_0.3.0
## [37] irlba_2.3.5.1 ExperimentHub_2.8.0
## [39] RSQLite_2.3.1 beachmat_2.16.0
## [41] labeling_0.4.2 filelock_1.0.2
## [43] fansi_1.0.4 httr_1.4.5
## [45] compiler_4.3.0 bit64_4.0.5
## [47] withr_2.5.0 BiocParallel_1.34.0
## [49] viridis_0.6.2 DBI_1.1.3
## [51] highr_0.10 biomaRt_2.56.0
## [53] MASS_7.3-59 rappdirs_0.3.3
## [55] DelayedArray_0.26.0 rjson_0.2.21
## [57] tools_4.3.0 vipor_0.4.5
## [59] beeswarm_0.4.0 interactiveDisplayBase_1.38.0
## [61] httpuv_1.6.9 glue_1.6.2
## [63] restfulr_0.0.15 promises_1.2.0.1
## [65] grid_4.3.0 Rtsne_0.16
## [67] cluster_2.1.4 generics_0.1.3
## [69] gtable_0.3.3 data.table_1.14.8
## [71] hms_1.1.3 BiocSingular_1.16.0
## [73] ScaledMatrix_1.8.0 metapod_1.8.0
## [75] xml2_1.3.3 utf8_1.2.3
## [77] XVector_0.40.0 ggrepel_0.9.3
## [79] BiocVersion_3.17.1 pillar_1.9.0
## [81] stringr_1.5.0 limma_3.56.0
## [83] later_1.3.0 dplyr_1.1.2
## [85] BiocFileCache_2.8.0 lattice_0.21-8
## [87] rtracklayer_1.60.0 bit_4.0.5
## [89] tidyselect_1.2.0 locfit_1.5-9.7
## [91] Biostrings_2.68.0 knitr_1.42
## [93] gridExtra_2.3 bookdown_0.33
## [95] ProtGenerics_1.32.0 edgeR_3.42.0
## [97] xfun_0.39 statmod_1.5.0
## [99] stringi_1.7.12 lazyeval_0.2.2
## [101] yaml_2.3.7 evaluate_0.20
## [103] codetools_0.2-19 tibble_3.2.1
## [105] BiocManager_1.30.20 cli_3.6.1
## [107] xtable_1.8-4 munsell_0.5.0
## [109] jquerylib_0.1.4 Rcpp_1.0.10
## [111] dbplyr_2.3.2 png_0.1-8
## [113] XML_3.99-0.14 parallel_4.3.0
## [115] ellipsis_0.3.2 blob_1.2.4
## [117] prettyunits_1.1.1 sparseMatrixStats_1.12.0
## [119] bitops_1.0-7 viridisLite_0.4.1
## [121] scales_1.2.1 purrr_1.0.1
## [123] crayon_1.5.2 rlang_1.1.0
## [125] cowplot_1.1.1 KEGGREST_1.40.0