Modern single-cell and genomic assays routinely produce matrices that
do not fit in memory. The 10x Genomics 1.3 million brain cell dataset,
for example, is a 27,998 x 1,306,127 matrix of UMI counts —
over 36 billion entries, or more than 130 GB as an ordinary dense
integer matrix. Even stored sparsely, datasets like this are awkward to
hold in an R session alongside the intermediate objects an analysis
creates.
DelayedArray
solves half of this problem: it lets an array live on disk
behind the ordinary R array API, so dim(), [,
rowSums(), and arithmetic all work without loading the data
into memory. What remains is the choice of backend —
the on-disk representation that actually stores and serves the data.
HDF5Array
provides an HDF5 backend and TileDBArray
a TileDB one.
DuckDBArray
provides a backend built on columnar Parquet files queried
through DuckDB. Rather than
treating a count matrix as a grid of blocks to be read and processed
chunk by chunk, DuckDBArray treats it as a database table and
lets DuckDB’s SQL engine do the work: a rowSums() becomes a
SUM ... GROUP BY, a rowVars() a
VAR_SAMP, and the query planner handles memory and
multi-threading for you. As with any DelayedArray
backend, none of these operations load the matrix into
memory — they read only what they need from disk.
This vignette is a practical introduction: how to create a
DuckDBMatrix, how the familiar Bioconductor array
operations behave on it, and when the DuckDB backend is a good fit. Two
companion vignettes go further:
A DuckDBMatrix is backed by a Parquet file in
coordinate (COO) form: one row per stored entry, with
key columns giving the position and a data column giving the value. For
a small example we write such a file with arrow, storing
a gene-by-cell count matrix as (gene, cell, count)
triples.
library(arrow)
counts_df <- expand.grid(gene = paste0("Gene", 1:100),
cell = paste0("Cell", 1:20),
stringsAsFactors = FALSE)
set.seed(1L)
counts_df$count <- rpois(nrow(counts_df), lambda = 3)
parquet_path <- tempfile(fileext = ".parquet")
write_parquet(counts_df, parquet_path)We construct the matrix by naming the data column and the key columns that index each dimension. The key columns carry the dimension names and their order:
mat <- DuckDBMatrix(parquet_path,
datacol = "count",
keycols = list(gene = paste0("Gene", 1:100),
cell = paste0("Cell", 1:20)))
mat
#> <100 x 20> sparse DuckDBMatrix object of type "integer":
#> Cell1 Cell2 Cell3 Cell4 ... Cell17 Cell18 Cell19 Cell20
#> Gene1 2 4 2 4 . 1 4 1 5
#> Gene2 2 2 2 1 . 0 4 2 1
#> Gene3 3 2 3 3 . 0 4 3 3
#> Gene4 5 8 2 3 . 6 1 3 5
#> Gene5 2 3 1 2 . 5 1 3 4
#> ... . . . . . . . . .
#> Gene96 4 3 4 2 . 2 3 4 4
#> Gene97 3 1 1 4 . 0 0 2 1
#> Gene98 2 5 1 2 . 4 1 0 2
#> Gene99 4 2 1 2 . 3 2 4 2
#> Gene100 3 4 6 1 . 1 2 3 2It looks and behaves like a matrix, but the data stays on disk:
dim(mat)
#> [1] 100 20
rowSums(mat)[1:5]
#> gene
#> Gene1 Gene2 Gene3 Gene4 Gene5
#> 65 54 65 78 61
colMeans(mat)[1:5]
#> cell
#> Cell1 Cell2 Cell3 Cell4 Cell5
#> 3.05 3.09 2.67 3.13 3.08Nothing above read the whole matrix into memory: each summary is evaluated as a DuckDB query over the Parquet file.
In practice you usually start from an existing (possibly sparse)
matrix. Use writeCoordArray() to serialize it to Parquet,
and createDimTables() to build the dimension lookup tables
that let DuckDB prune to the rows and columns a query touches. Storing
the matrix transposed puts features in columns, which Parquet’s
columnar layout makes cheap to scan feature-by-feature.
library(Matrix)
m <- Matrix(rpois(200 * 50, lambda = 1), nrow = 200, ncol = 50, sparse = TRUE)
rownames(m) <- paste0("Gene", seq_len(nrow(m)))
colnames(m) <- paste0("Cell", seq_len(ncol(m)))
path <- tempfile()
writeCoordArray(m, path)
dimtbls <- createDimTables(m)
mat2 <- DuckDBMatrix(path, datacol = "value",
keycols = list(index1 = setNames(seq_len(nrow(m)), rownames(m)),
index2 = setNames(seq_len(ncol(m)), colnames(m))),
dimtbls = dimtbls)
dim(mat2)
#> [1] 200 50Subsetting is delayed: it records the selection and is resolved only when values are needed.
sub <- mat[1:5, 1:4]
sub
#> <5 x 4> sparse DuckDBMatrix object of type "integer":
#> cell
#> gene Cell1 Cell2 Cell3 Cell4
#> Gene1 2 4 2 4
#> Gene2 2 2 2 1
#> Gene3 3 2 3 3
#> Gene4 5 8 2 3
#> Gene5 2 3 1 2
as.matrix(sub)
#> cell
#> gene Cell1 Cell2 Cell3 Cell4
#> Gene1 2 4 2 4
#> Gene2 2 2 2 1
#> Gene3 3 2 3 3
#> Gene4 5 8 2 3
#> Gene5 2 3 1 2DuckDBArray implements the MatrixGenerics row/column summaries as SQL aggregations, so the usual functions work directly:
library(MatrixGenerics)
rowSums(mat)[1:5]
#> gene
#> Gene1 Gene2 Gene3 Gene4 Gene5
#> 65 54 65 78 61
rowVars(mat)[1:5]
#> gene
#> Gene1 Gene2 Gene3 Gene4 Gene5
#> 2.197368 3.168421 3.881579 3.884211 3.418421
colSds(mat)[1:5]
#> cell
#> Cell1 Cell2 Cell3 Cell4 Cell5
#> 1.465943 1.639752 1.620918 1.862197 1.835013Two helpers are added for count-based feature selection —
rowNnzs() / colNnzs() count non-zero entries
(detection rate), and rowDeviances() computes
binomial/Poisson deviance for highly-variable-gene selection (Townes et
al. 2019):
Because only stored entries live in the Parquet file, sparse matrices
stay compact end to end. A DuckDBMatrix reports its
sparsity and coerces to the familiar sparse containers when you need
them in memory:
A DuckDBMatrix is a DelayedArray,
so delayed operations accumulate without touching the data, and
materialize only on demand:
Most work never materializes the matrix: reductions and grouped
statistics push down to SQL (see Matrix statistics), so they
stay memory-bounded regardless of size. Whole-object coercions like
as.matrix() or as(x, "dgCMatrix"), by
contrast, build the entire result at once — ideal when
it fits, but not when it doesn’t.
When you need to consume a result too large for memory — to feed a
chunked/minibatch step, or to re-write it to another on-disk backend —
realize it a block at a time instead. Because a
DuckDBMatrix is a DelayedArray,
the standard block machinery applies: blockApply() reads
and processes the matrix one block at a time, and each block is a
separate, pruned query (only that block’s entries are
read), so peak memory is bounded by a single block rather than the full
result:
# one block in memory at a time; each block is a separate, pruned query
block_sums <- blockApply(mat, function(block) sum(block))
sum(unlist(block_sums))
#> [1] 5960realize(x, BACKEND = "HDF5Array") uses the same
block-streaming path to write a large result to another on-disk format
without ever holding it whole, and
DelayedArray::setAutoBlockSize() tunes how much is read per
block.
The same backend serves arrays of any dimension through
DuckDBArray(); a DuckDBMatrix is simply the
two-dimensional case. Here the key columns index three dimensions:
arr_df <- expand.grid(i = 1:4, j = 1:3, k = 1:2)
arr_df$value <- rpois(nrow(arr_df), lambda = 2)
arr_path <- tempfile(fileext = ".parquet")
write_parquet(arr_df, arr_path)
arr <- DuckDBArray(arr_path, datacol = "value",
keycols = list(i = 1:4, j = 1:3, k = 1:2))
dim(arr)
#> [1] 4 3 2
arr[1:2, 1:2, 1]
#> <2 x 2> sparse DuckDBArray object of type "integer":
#> j
#> i 1 2
#> 1 1 2
#> 2 0 1DuckDBArray is a good fit when:
Other backends remain the better choice when data fits comfortably in
memory (an in-memory dgCMatrix is hard to beat for small
data), or for workloads built around dense linear algebra or heavy
random element access.
For a head-to-head comparison with in-memory, HDF5Array, and TileDBArray backends — including where each one wins — see the Benchmarking DuckDBArray vignette. To understand or extend the backend itself, see Implementing the DuckDBArray backend.
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] arrow_25.0.0 DuckDBArray_0.99.2 DelayedArray_0.39.3
#> [4] SparseArray_1.13.2 S4Arrays_1.13.0 abind_1.4-8
#> [7] MatrixGenerics_1.25.0 matrixStats_1.5.0 Matrix_1.7-5
#> [10] DuckDBDataFrame_0.99.6 IRanges_2.47.2 S4Vectors_0.51.5
#> [13] BiocGenerics_0.59.10 generics_0.1.4 bit64_4.8.2
#> [16] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.10 lattice_0.22-9 digest_0.6.39
#> [4] magrittr_2.0.5 sparseMatrixStats_1.25.0 evaluate_1.0.5
#> [7] grid_4.6.1 blob_1.3.0 fastmap_1.2.0
#> [10] jsonlite_2.0.0 DBI_1.3.0 BiocManager_1.30.27
#> [13] purrr_1.2.2 jquerylib_0.1.4 duckdb_1.5.4.3
#> [16] cli_3.6.6 rlang_1.3.0 dbplyr_2.6.0
#> [19] XVector_0.53.0 withr_3.0.3 cachem_1.1.0
#> [22] yaml_2.3.12 otel_0.2.0 tools_4.6.1
#> [25] dplyr_1.2.1 assertthat_0.2.1 buildtools_1.0.0
#> [28] vctrs_0.7.3 R6_2.6.1 lifecycle_1.0.5
#> [31] bit_4.6.0 pkgconfig_2.0.3 bslib_0.11.0
#> [34] pillar_1.11.1 Rcpp_1.1.2 glue_1.8.1
#> [37] xfun_0.60 tibble_3.3.1 tidyselect_1.2.1
#> [40] sys_3.4.3 knitr_1.51 htmltools_0.5.9
#> [43] rmarkdown_2.31 maketools_1.3.2 compiler_4.6.1