DuckDBArray 0.99.7
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:
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("DuckDBArray")library(DuckDBArray)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):
rowNnzs(mat)[1:5]
#> gene
#> Gene1 Gene2 Gene3 Gene4 Gene5
#> 20 18 19 20 20
rowDeviances(mat, family = "binomial")[1:5]
#> [1] 14.46593 26.66845 25.31237 19.76327 20.57861Because 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:
is_sparse(mat2)
#> [1] TRUE
as(mat2[1:5, 1:4], "dgCMatrix")
#> 5 x 4 sparse Matrix of class "dgCMatrix"
#> index2
#> index1 Cell1 Cell2 Cell3 Cell4
#> Gene1 2 1 . 1
#> Gene2 3 3 . .
#> Gene3 2 . 3 2
#> Gene4 1 1 2 1
#> Gene5 . 3 . .A DuckDBMatrix is a DelayedArray, so delayed operations
accumulate without touching the data, and materialize only on demand:
transformed <- log1p(mat) * 2
class(transformed)
#> [1] "DuckDBMatrix"
#> attr(,"package")
#> [1] "DuckDBArray"
as.matrix(transformed[1:3, 1:3])
#> cell
#> gene Cell1 Cell2 Cell3
#> Gene1 2.197225 3.218876 2.197225
#> Gene2 2.197225 2.197225 2.197225
#> Gene3 2.772589 2.197225 2.772589Most 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 1The example above already knows keycols because it wrote arr_path itself.
Given a COO Parquet file produced by another tool, use arrow to
find the key and data column names first: the schema names the columns and
their types.
open_dataset(arr_path)$schema
#> Schema
#> i: int32
#> j: int32
#> k: int32
#> value: int32
#>
#> See $metadata for additional Schema metadatascanKeycols() then scans the named key columns for their distinct values,
sorted, and returns them already in the keycols shape the
DuckDBArray()/DuckDBMatrix() constructors expect, so there’s no need to
repeat the distinct()/arrange()/collect() steps by hand for each key
column:
scanKeycols(arr_path, keycols = c("i", "j", "k"))
#> $i
#> [1] 1 2 3 4
#>
#> $j
#> [1] 1 2 3
#>
#> $k
#> [1] 1 2arr2 <- DuckDBArray(arr_path, datacol = "value",
keycols = scanKeycols(arr_path, c("i", "j", "k")))
dim(arr2)
#> [1] 4 3 2DuckDBArray 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 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.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] arrow_25.0.0 DuckDBArray_0.99.7 DelayedArray_0.39.5
#> [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-6
#> [10] DuckDBDataFrame_0.99.20 IRanges_2.47.2 S4Vectors_0.51.6
#> [13] BiocGenerics_0.59.12 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.23-1 digest_0.6.39
#> [4] magrittr_2.0.5 sparseMatrixStats_1.25.0 evaluate_1.0.5
#> [7] grid_4.6.1 bookdown_0.47 blob_1.3.0
#> [10] fastmap_1.2.0 jsonlite_2.0.0 DBI_1.3.0
#> [13] BiocManager_1.30.27 purrr_1.2.2 jquerylib_0.1.4
#> [16] duckdb_1.5.5 cli_3.6.6 rlang_1.3.0
#> [19] dbplyr_2.6.0 XVector_0.53.0 withr_3.0.3
#> [22] cachem_1.1.0 yaml_2.3.12 otel_0.2.0
#> [25] tools_4.6.1 dplyr_1.2.1 assertthat_0.2.1
#> [28] vctrs_0.7.3 R6_2.6.1 lifecycle_1.0.5
#> [31] bit_4.6.0 pkgconfig_2.0.3 bslib_0.12.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] knitr_1.51 htmltools_0.5.9 rmarkdown_2.31
#> [43] compiler_4.6.1