---
title: "Benchmarking BiocDuckDB"
author:
- name: Patrick Aboyoun
  email: aboyounp@gene.com
  affiliation: Genentech, Inc.
date: "Compiled: `r BiocStyle::doc_date()`; Modified: 6 July 2026"
package: BiocDuckDB
vignette: >
  %\VignetteIndexEntry{2. Benchmarking BiocDuckDB}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    number_sections: true
    toc: true
    toc_depth: 2
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      error = FALSE, warning = FALSE, message = FALSE)
```

# 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?**

`r Biocpkg("BiocDuckDB")` implements the common `r Biocpkg("scuttle")` and
`r Biocpkg("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 `r Biocpkg("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.

**`r Biocpkg("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`) |

**`r Biocpkg("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 |

**Dimensionality reduction is a different kind of case.** `r Biocpkg("BiocSingular")`'s
`runSVD()`, and so `r Biocpkg("scran")`'s `fixedPCA()` and `r Biocpkg("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 `r Biocpkg("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.

```{r live-demo}
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)
```

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 `r Biocpkg("ExperimentHub")`, accession `EH1039`), subset to 12,500 cells,
the size used by the original comparison. Each operation runs on three
backends: an in-memory `r CRANpkg("Matrix")` `dgCMatrix`, an
`r Biocpkg("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).

```{r results, echo=FALSE, results='asis'}
helper <- system.file("scripts", "make_timings_table.R", package = "BiocDuckDB")
if (nzchar(helper)) {
    source(helper)
    res <- load_vignette_timings()
} else {
    res <- NULL
}
if (is.null(res)) {
    cat("_Precomputed benchmark results are not available in this build; ",
        "generate them with `inst/scripts/run_scran_scuttle_benchmarks.R`._\n", sep = "")
} else {
    print(make_timings_table(
        caption = "Elapsed seconds per operation. Speedups > 1 favor DuckDB.",
        results = res))
    cat("\n\n")
    timings_config_note(res)
}
```

# Summary

`perCellQCMetrics()`, `modelGeneVar()`, `findMarkers()`, and the rest of the
`r Biocpkg("scran")`/`r Biocpkg("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 `r Biocpkg("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 `r Biocpkg("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
`r Biocpkg("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 `r Biocpkg("DuckDBArray")` benchmarking vignette.

# Session information

```{r sessioninfo}
sessionInfo()
```