BiocDuckDB 0.99.22
The Introduction to BiocDuckDB vignette shows how an experiment can keep its assays on disk while presenting the standard Bioconductor API. This vignette asks the follow-up question: can the standard single-cell analysis methods run directly on that on-disk representation, and how fast?
BiocDuckDB implements the common scuttle and
scran generics for DuckDBMatrix as SQL-optimized queries, so QC,
normalization, variance modelling, and marker detection run on the Parquet-backed
matrix without realizing it into memory. We compare those against the same generics
on an in-memory dgCMatrix and on HDF5Array.
The headline results below were produced offline on the 10x Genomics 1.3 million brain-cell dataset (see Benchmark setup) and are rendered here from a bundled results file, so this vignette builds quickly. The A small, live comparison section runs a miniature version at build time.
Each method is implemented as per-gene / per-group SQL aggregation on the
DuckDBMatrix, so the scan and the arithmetic happen in DuckDB and only the
(small) result crosses back into R.
scuttle, QC, normalization, pseudo-bulk:
| Function | Description |
|---|---|
perCellQCMetrics / perFeatureQCMetrics |
library size, detected genes / mean, detection rate |
librarySizeFactors / normalizeCounts |
size factors and normalization |
summarizeAssayByGroup |
pseudo-bulk aggregation (GROUP BY) |
scran, variance modelling and markers:
| Function | Description |
|---|---|
modelGeneVar / modelGeneVarByPoisson |
decompose technical vs biological variance |
correlatePairs |
pairwise gene correlations |
pairwiseTTests / findMarkers |
pairwise DE and candidate markers |
To show the mechanics without a large download, we build a small sparse matrix and
run one QC metric on both an in-memory dgCMatrix and a DuckDBMatrix, confirming
the results agree.
library(BiocDuckDB)
library(DuckDBArray)
library(Matrix)
library(scuttle)
set.seed(1L)
m <- as(Matrix(rpois(2000 * 400, lambda = 0.3), nrow = 2000, ncol = 400,
sparse = TRUE), "dgCMatrix")
rownames(m) <- paste0("Gene", seq_len(nrow(m)))
colnames(m) <- paste0("Cell", seq_len(ncol(m)))
path <- tempfile()
writeParquet(t(m), path)
mt <- t(m)
mat <- DuckDBMatrix(path, datacol = "value",
keycols = list(index2 = setNames(seq_len(ncol(mt)), colnames(mt)),
index1 = setNames(seq_len(nrow(mt)), rownames(mt))),
dimtbls = createDimTables(mt))
## same answer, one in memory and one queried from disk
qc_mem <- perCellQCMetrics(m)
qc_ddb <- perCellQCMetrics(mat)
all.equal(qc_mem$sum, qc_ddb$sum)
#> [1] TRUEAt this size the in-memory matrix is faster; there is nothing to gain from going to disk. The advantage appears at scale, which is what the offline benchmark measures.
The full benchmark uses the 10x Genomics 1.3 million brain-cell dataset (available
through ExperimentHub, accession EH1039), subset to 12,500 cells,
the size used by the original comparison. Each operation runs on three
backends: an in-memory Matrix dgCMatrix, an
HDF5Array, and a DuckDBMatrix. The in-memory and HDF5 backends run
single-threaded; DuckDBMatrix autotunes DuckDB’s internal threads up to the core
budget. The variance and marker operations run on log-normalized counts produced by
each backend.
The table below reports the offline benchmark’s measured timings, bundled with the
package as inst/scripts/benchmark_results.rds. Regenerate it on your own
hardware with inst/scripts/run_scran_scuttle_benchmarks.R (see that script’s
header).
| Operation | In-memory (s) | HDF5Array (s) | DuckDB (s) | vs HDF5Array (x) | vs in-memory (x) |
|---|---|---|---|---|---|
| perCellQCMetrics | 0.56 | 3.09 | 0.10 | 30.6 | 5.5 |
| perFeatureQCMetrics | 0.77 | 4.44 | 0.09 | 48.3 | 8.4 |
| summarizeAssayByGroup | 0.11 | 1.25 | 0.58 | 2.2 | 0.2 |
| normalizeCounts | 0.78 | 1.95 | 0.15 | 13.0 | 5.2 |
| modelGeneVar | 0.28 | 10.42 | 1.02 | 10.3 | 0.3 |
| correlatePairs | 54.98 | 77.66 | 0.61 | 127.1 | 90.0 |
| pairwiseTTests | 6.98 | 15.95 | 7.11 | 2.2 | 1.0 |
| findMarkers | 21.67 | 31.17 | 22.16 | 1.4 | 1.0 |
Configuration: 27,998 genes x 12,500 cells, 16-core budget. In-memory: dgCMatrix. HDF5Array: 10x/HDF5 backend. DuckDB: DuckDBMatrix over Parquet (autotuned threads). All backends run the same scran/scuttle generics.
perCellQCMetrics(), modelGeneVar(), findMarkers(), and the rest of the
scran/scuttle generics used here dispatch to
SQL-optimized methods for DuckDBMatrix. That means existing analysis code that
calls these generics keeps working when the matrix happens to be Parquet-backed,
with no code changes and no realization into memory.
Against HDF5Array, the fair comparison since both keep the matrix on
disk, DuckDB is faster on every operation measured, from about 1.4x on
findMarkers to well over 100x on correlatePairs. perCellQCMetrics,
perFeatureQCMetrics, and normalizeCounts are pure SUM/AVG aggregations,
which DuckDB pushes straight into a columnar scan; that lets them beat even the
in-memory dgCMatrix by several times, without ever loading the matrix.
correlatePairs is the standout: its sparse-aware SQL avoids the dense
intermediate matrices the other backends build for a pairwise correlation, which
is what turns the usual disk-vs-memory gap into a roughly two-order-of-magnitude
speedup over both HDF5Array and in-memory.
Some steps still favor in-memory at this scale. summarizeAssayByGroup and
modelGeneVar are faster on an in-memory dgCMatrix (though DuckDB still beats
HDF5Array on both); a dense in-memory group-by can outrun even a
single disk scan while the matrix comfortably fits in RAM, and the DuckDB
advantage should grow as the data outgrows it. Marker detection
(pairwiseTTests, findMarkers) lands on par with in-memory, since the per-gene
statistics computed after the SQL GROUP BY are identical work regardless of
which backend did the aggregation.
The table in the Results section above has the exact measured numbers. There is
no per-call tuning needed to get them: DuckDBMatrix autotunes DuckDB’s internal
threads up to the available core budget on its own. The lever that matters is
structural, covered next.
The value is compounding: because these methods run on the DuckDB-backed object directly, an analysis can go from raw counts through QC, normalization, feature selection, and marker detection before ever realizing the matrix into memory, realizing only the small, filtered result it actually needs. That is what makes the filter, realize, analyze pattern from the introduction practical on datasets far larger than RAM.
inst/scripts/run_scran_scuttle_benchmarks.R reproduces these numbers on your
own hardware (BENCH_NCELLS, BENCH_CORES; set BENCH_SYNTHETIC=1 to
smoke-test without the EH1039 download). It writes benchmark_results.rds to
your working directory, in the same format inst/scripts/make_timings_table.R
reads to build the table in the Results section above. Running the script does
not, by itself, change that table: every installation of a given package
release bundles the same fixed, precomputed inst/scripts/benchmark_results.rds,
regardless of who built it or on what hardware. To have your own run replace the
bundled numbers, copy your output over inst/scripts/benchmark_results.rds in
the package source and rebuild the package.
For the lower-level matrix operations (colSums, rowVars, rowDeviances), see
the DuckDBArray benchmarking vignette.
sessionInfo()
#> R version 4.6.1 Patched (2026-06-24 r90190)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Ventura 13.7.8
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] scuttle_1.23.2 SingleCellExperiment_1.35.2
#> [3] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [5] BiocDuckDB_0.99.22 DuckDBGRanges_0.99.8
#> [7] GenomicRanges_1.65.4 Seqinfo_1.3.2
#> [9] DuckDBArray_0.99.8 DelayedArray_0.39.6
#> [11] SparseArray_1.13.2 S4Arrays_1.13.0
#> [13] abind_1.4-8 MatrixGenerics_1.25.0
#> [15] matrixStats_1.5.0 Matrix_1.7-6
#> [17] DuckDBDataFrame_0.99.26 IRanges_2.47.5
#> [19] S4Vectors_0.51.9 BiocGenerics_0.59.12
#> [21] generics_0.1.4 bit64_4.8.6
#> [23] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 blob_1.3.0
#> [3] dplyr_1.2.1 arrow_25.0.1
#> [5] fastmap_1.2.0 bluster_1.23.1
#> [7] duckdb_1.5.5 digest_0.6.39
#> [9] rsvd_1.0.5 lifecycle_1.0.5
#> [11] cluster_2.1.8.3 sf_1.1-2
#> [13] statmod_1.5.2 magrittr_2.0.5
#> [15] compiler_4.6.1 rlang_1.3.0
#> [17] sass_0.4.10 tools_4.6.1
#> [19] igraph_2.3.3 yaml_2.3.12
#> [21] knitr_1.52 dqrng_0.4.1
#> [23] bit_4.6.0 classInt_0.4-11
#> [25] BiocParallel_1.47.0 KernSmooth_2.23-27
#> [27] withr_3.0.3 purrr_1.2.2
#> [29] grid_4.6.1 beachmat_2.29.2
#> [31] e1071_1.7-17 edgeR_4.99.4
#> [33] MultiAssayExperiment_1.39.1 cli_3.6.6
#> [35] rmarkdown_2.32 otel_0.2.0
#> [37] metapod_1.21.0 rjson_0.2.23
#> [39] DBI_1.3.0 cachem_1.1.0
#> [41] proxy_0.4-29 assertthat_0.2.1
#> [43] parallel_4.6.1 BiocManager_1.30.27
#> [45] XVector_0.53.0 vctrs_0.7.3
#> [47] jsonlite_2.0.0 bookdown_0.48
#> [49] BiocSingular_1.29.1 BiocNeighbors_2.7.3
#> [51] irlba_2.3.7 magick_2.9.1
#> [53] locfit_1.5-9.12 limma_3.99.0
#> [55] jquerylib_0.1.4 units_1.0-1
#> [57] glue_1.8.1 codetools_0.2-20
#> [59] ScaledMatrix_1.21.0 tibble_3.3.1
#> [61] pillar_1.11.1 htmltools_0.5.9
#> [63] MultiAssaySpatialExperiment_0.99.12 R6_2.6.1
#> [65] dbplyr_2.6.0 evaluate_1.0.5
#> [67] lattice_0.23-1 SpatialExperiment_1.23.0
#> [69] scran_1.41.1 bslib_0.12.0
#> [71] class_7.3-24 Rcpp_1.1.2
#> [73] xfun_0.60 pkgconfig_2.0.3