Introduction

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.

What BiocDuckDB optimizes

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:

scran, variance modelling and markers:

Dimensionality reduction is a different kind of case. BiocSingular’s runSVD(), and so scran’s fixedPCA() and scater’s runPCA()/calculatePCA(), already compute a correct PCA directly on a DuckDBMatrix with the ordinary BSPARAM = BiocSingular::IrlbaParam(), through the SQL-pushdown %*%/crossprod methods DuckDBArray implements. But unlike the aggregation-only functions above, PCA is not itself pure SQL: irlba’s Lanczos iteration calls %*% once per solver step, so the query-planning cost of each call is paid repeatedly rather than once, and that path is markedly slower on DuckDBMatrix than in memory.

BSPARAM = DuckDBIrlbaParam() closes that gap: it materializes the (already HVG-subsetted) matrix once into an in-memory sparse matrix and drives irlba directly on it, the same size-gated materialization loadIntoMemory() uses elsewhere in this package, instead of paying the SQL-pushdown cost once per iteration. It falls back to ordinary IrlbaParam behavior whenever it doesn’t apply (not a DuckDBMatrix, too large for "memory_limit"), so it is safe to use as a drop-in replacement for IrlbaParam() wherever a DuckDBMatrix might show up. The benchmark below uses DuckDBIrlbaParam() for exactly this reason: to measure the path that is actually meant to be fast, not the one that only demonstrates correctness.

A small, live comparison

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] TRUE

At 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.

Benchmark setup

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, and calculatePCA (10 components, top 200 HVGs, BSPARAM = DuckDBIrlbaParam()), run on log-normalized counts produced by each backend.

Results

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).

Configuration: 27,998 genes x 12,500 cells, 19-core budget. In-memory: dgCMatrix. HDF5Array: 10x/HDF5 backend. DuckDB: DuckDBMatrix over Parquet (autotuned threads). All backends run the same scran/scuttle generics.

Summary

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.

calculatePCA is run here with BSPARAM = DuckDBIrlbaParam() (see What BiocDuckDB optimizes), not the ordinary BSPARAM = IrlbaParam(). The distinction matters: with the ordinary BSPARAM, each Lanczos iteration re-pays SQL query-planning overhead that a single-scan aggregation only pays once, measured at about 50x slower than in-memory on a real 12,500-cell dataset, markedly worse than every other operation in this table. DuckDBIrlbaParam() closes that gap by materializing the matrix once and driving irlba directly on it, so the number in the table above reflects that fast path, not the structural limit described for the ordinary one.

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.

When this matters

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.

Running your own benchmarks

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.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 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.32.so;  LAPACK version 3.12.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] scuttle_1.23.2              SingleCellExperiment_1.35.2
#>  [3] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [5] BiocDuckDB_0.99.23          DuckDBGRanges_0.99.8       
#>  [7] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [9] DuckDBArray_0.99.9          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] DBI_1.3.0                           rlang_1.3.0                        
#>  [3] magrittr_2.0.5                      otel_0.2.0                         
#>  [5] e1071_1.7-17                        compiler_4.6.1                     
#>  [7] vctrs_0.7.3                         pkgconfig_2.0.3                    
#>  [9] SpatialExperiment_1.23.0            fastmap_1.2.0                      
#> [11] dbplyr_2.6.0                        magick_2.9.1                       
#> [13] XVector_0.53.0                      rmarkdown_2.32                     
#> [15] purrr_1.2.2                         bit_4.6.0                          
#> [17] xfun_0.60                           MultiAssayExperiment_1.39.1        
#> [19] bluster_1.23.1                      cachem_1.1.0                       
#> [21] beachmat_2.29.2                     jsonlite_2.0.0                     
#> [23] blob_1.3.0                          BiocParallel_1.47.0                
#> [25] irlba_2.3.7                         parallel_4.6.1                     
#> [27] cluster_2.1.8.3                     R6_2.6.1                           
#> [29] MultiAssaySpatialExperiment_0.99.12 bslib_0.12.0                       
#> [31] limma_3.99.0                        jquerylib_0.1.4                    
#> [33] Rcpp_1.1.2                          assertthat_0.2.1                   
#> [35] knitr_1.52                          igraph_2.3.3                       
#> [37] tidyselect_1.2.1                    yaml_2.3.12                        
#> [39] codetools_0.2-20                    lattice_0.23-1                     
#> [41] tibble_3.3.1                        withr_3.0.3                        
#> [43] evaluate_1.0.5                      sf_1.1-2                           
#> [45] units_1.0-1                         proxy_0.4-29                       
#> [47] pillar_1.11.1                       BiocManager_1.30.27                
#> [49] KernSmooth_2.23-27                  class_7.3-24                       
#> [51] glue_1.8.1                          metapod_1.21.0                     
#> [53] maketools_1.3.2                     tools_4.6.1                        
#> [55] BiocNeighbors_2.7.3                 sys_3.4.3                          
#> [57] ScaledMatrix_1.21.0                 locfit_1.5-9.12                    
#> [59] buildtools_1.0.0                    scran_1.41.1                       
#> [61] grid_4.6.1                          edgeR_4.99.5                       
#> [63] duckdb_1.5.5                        BiocSingular_1.29.1                
#> [65] cli_3.6.6                           rsvd_1.0.5                         
#> [67] arrow_25.0.1                        dplyr_1.2.1                        
#> [69] sass_0.4.10                         digest_0.6.39                      
#> [71] classInt_0.4-11                     dqrng_0.4.1                        
#> [73] rjson_0.2.23                        htmltools_0.5.9                    
#> [75] lifecycle_1.0.5                     statmod_1.5.2