---
title: "Introduction to 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{1. Introduction to 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)
has_airway <- requireNamespace("airway", quietly = TRUE)
```

# Introduction

Bioconductor's experiment containers, `SummarizedExperiment`,
`SingleCellExperiment`, `MultiAssayExperiment`, hold their assays, annotations,
and reduced dimensions in memory. For modern datasets that is increasingly
impractical: a million-cell single-cell experiment is tens of gigabytes, and a
multi-omics study multiplies that across assays.

`r Biocpkg("BiocDuckDB")` is the integration layer of the **BiocDuckDB** suite. It
provides two high-level functions,

- **`writeParquet()`** serializes a Bioconductor experiment object to a structured
  directory of columnar **Parquet** files, and
- **`readParquet()`** reads that directory back as a **DuckDB-backed** experiment
  object.

So an experiment keeps its assays on disk while presenting the *same*
Bioconductor API. Assays become `r Biocpkg("DuckDBArray")` matrices, `rowData` and
`colData` become `r Biocpkg("DuckDBDataFrame")`s, and `rowRanges` becomes a
`r Biocpkg("DuckDBGRanges")` list, and operations on them run as SQL against the
Parquet files. The on-disk format is a self-describing
[Frictionless Data Package](https://frictionlessdata.io/), readable by other tools
and languages, not just R.

This vignette covers the round-trip workflow and the operations it enables. For the
SQL-optimized `r Biocpkg("scran")`/`r Biocpkg("scuttle")` analysis methods and how
they perform at scale, see *Benchmarking BiocDuckDB*.

## Installation

```{r install, eval=FALSE}
if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("BiocDuckDB")
```

```{r load}
library(BiocDuckDB)
library(SummarizedExperiment)
```

# The round-trip

We write a `SummarizedExperiment` to Parquet and read it back as a DuckDB-backed
object. The `r Biocpkg("airway")` bulk RNA-seq dataset is a `RangedSummarizedExperiment`,
so it also exercises the genomic-coordinate path.

```{r airway-check, echo=FALSE, results='asis'}
if (!has_airway) {
    cat("_The `airway` package is not installed; the code below is shown but",
        "not evaluated._\n")
}
```

```{r roundtrip, eval=has_airway}
data(airway, package = "airway")

se_path <- file.path(tempdir(), "airway_se")
writeParquet(airway, se_path)

airway_ddb <- readParquet(se_path)
class(airway_ddb)
dim(airway_ddb)
```

## The storage layout

`writeParquet()` lays out each component as its own Parquet under a
Frictionless-style directory, described by a `datapackage.json`:

```{r layout, eval=has_airway}
head(list.files(se_path, recursive = TRUE), 10)
```

Assays are stored in **coordinate (COO) form**, one row per stored entry with
`__feature__`/`__sample__` index columns and a `value`, which is compact for the
sparse matrices typical of single-cell data, and lets DuckDB scan only what a query
touches.

## Advanced: writing and attaching resources directly

The rest of this section is for interoperating with other tools or writing data
incrementally; skip ahead to "Operations run as SQL" below if you only need the
`writeParquet()`/`readParquet()` round trip above.

The `datapackage.json` that `writeParquet()` writes follows the
[Frictionless Data Package](https://datapackage.org) specification, a widely used,
language-agnostic way to describe a directory of data files. BiocDuckDB adds a few
of its own fields on top (`dimension` and `layout`, used below), but the file
otherwise reads like any other Frictionless package. In practice that means two
things: another tool, or another language, can write COO Parquet files with a
matching `datapackage.json` and `readParquet()` will read them directly, with no R
object involved in producing them; and R code that already has COO Parquet on disk
from some other process can attach to it directly (below), instead of writing it
out again through `writeParquet()`.

To build a `datapackage.json` from pieces written separately, for example a
dataset streamed in one part at a time, or migrated from another store, collect the
descriptor each `writeParquet()` call returns and pass them all to
`writeDatapackage()`:

```{r ingest, eval=FALSE}
resources <- c(
    writeParquet(features, file.path(dir, "features"), dimension = "feature"),
    writeParquet(assays,   dir,                        dimension = "crossed"))
writeDatapackage("summarized_experiment", resources, dir)
```

### Streaming a flat resource from a cursor

The pattern above assumes each piece already exists as a real in-memory object
(a `data.frame`, an array). A source with no natural in-memory array
representation -- a database cursor, an API page, anything read in a loop --
is a different case, and `writeStreamingResource()` is the dedicated writer
for it: it owns the offset/part/append bookkeeping, types the row index wide
enough for a resource that will end up over 2^31 rows (`index_max`), and
writes an `_INCOMPLETE` marker so a crash mid-stream can never be read back as
a complete resource by mistake.

```{r stream, eval=FALSE}
cur <- DBI::dbSendQuery(con, "SELECT * FROM samples")
res <- writeStreamingResource(
    blocks = function(i) {
        chunk <- DBI::dbFetch(cur, n = 10000)
        if (nrow(chunk) == 0) NULL else chunk   # NULL signals end of stream
    },
    path = file.path(dir, "samples"), dimension = "sample",
    index_max = expected_row_count)
```

`blocks()` is called with the 1-based part index and must return the next
block as a `data.frame`, or `NULL` once exhausted; empty blocks are skipped
without consuming a part number. This covers a flat (non-array) resource --
for a hive-partitioned coordinate array (a sparse assay) streamed from the
same kind of source, see the "Streaming a hive-partitioned coordinate array"
paragraph on `?writeParquet`, which documents the equivalent
`writeParquet(slab, ..., append = TRUE, along = , offset = , group_offset = )`
recipe.

Each slab's own `grid` still defaults to a dense element budget, so for a
sparse slab call `setAutoBlockSize()`/`setAutoBlockShape("scale")` *before*
that `writeParquet()` call, not after -- `grid = NULL` sizes itself from
whatever is currently set:

```{r stream-grid, eval=FALSE}
density <- Matrix::nnzero(slab) / prod(dim(slab))
setAutoBlockSize(target_bytes_per_partition / density)
setAutoBlockShape("scale")
writeParquet(slab, path, append = TRUE, along = , offset = , group_offset = )
```

See `?writeCoordArray`'s "Sizing grid for large or sparse arrays" section for
why the default needs this correction and the hard ceiling on how far it can
be pushed.

Going the other way, the `DuckDBMatrix()` / `DuckDBArray()` / `DuckDBTable()`
constructors attach an existing COO Parquet in place, wrapping it as a lazy,
out-of-core object without copying or materializing it:

```{r attach, eval=FALSE}
mat <- DuckDBMatrix(coo_path, datacol = "value",
                    keycols = c("__feature__", "__sample__"))
```

`coo_path` is the directory holding the COO Parquet; `datacol` names the column
holding the stored values, and `keycols` names the two index columns that locate
each value in the matrix. The documented `model` values and the `dimension` x
`layout` dispatch that maps each resource to its accessor are described on
`?writeParquet` and `?readParquet`.

## Operations run as SQL

The DuckDB-backed object supports the standard API; component access, subsetting,
and summaries are pushed to DuckDB:

```{r operations, eval=has_airway}
assay(airway_ddb, "counts")
colData(airway_ddb)[, 1:3]

## subset on disk, then summarize
sub <- airway_ddb[1:1000, 1:4]
colSums(assay(sub, "counts"))
```

Because only stored values and the touched columns are read, the in-memory
footprint of the object stays small relative to the full data:

```{r memory, eval=has_airway}
c(in_memory = format(object.size(airway), units = "MB"),
  duckdb    = format(object.size(airway_ddb), units = "MB"))
```

## Genomic coordinates are preserved

`writeParquet()` stores a `RangedSummarizedExperiment`'s `rowRanges` as Parquet
`LIST[]` columns, so they come back as a `r Biocpkg("DuckDBGRanges")` list with the
genomic API intact:

```{r ranges, eval=has_airway}
rr <- rowRanges(airway_ddb)
class(rr)
elementNROWS(rr)[1:8]      # exons per gene, queried from disk
```

# Single-cell data

Single-cell experiments are where the on-disk representation pays off most: the
matrices are large and sparse, so COO storage compresses well and the DuckDB-backed
object holds essentially only metadata. A small sparse `SingleCellExperiment`
illustrates the round-trip:

```{r sce}
library(SingleCellExperiment)
library(Matrix)

set.seed(1L)
counts <- Matrix(rpois(2000 * 500, lambda = 0.3), nrow = 2000, ncol = 500,
                 sparse = TRUE)
rownames(counts) <- paste0("Gene", seq_len(nrow(counts)))
colnames(counts) <- paste0("Cell", seq_len(ncol(counts)))
sce <- SingleCellExperiment(assays = list(counts = counts))

sce_path <- file.path(tempdir(), "demo_sce")
writeParquet(sce, sce_path)
sce_ddb <- readParquet(sce_path)

## QC filtering stays on disk (SQL), then summarize
totals <- colSums(assay(sce_ddb, "counts"))
keep <- sce_ddb[, totals > median(totals)]
dim(keep)
```

# The filter, realize, analyze pattern

For data larger than memory, the idiomatic workflow uses DuckDB for the parts it
does well, filtering and summarizing on disk, and only then realizes the
manageable subset into memory for the standard analysis pipeline:

```{r pattern, eval=FALSE}
library(scran)
library(scater)

sce_ddb <- readParquet("data/raw_counts")          # DuckDB-backed, low RAM

## 1. filter on disk (SQL-optimized)
keep <- colSums(assay(sce_ddb, "counts")) > 1000
sce_ddb <- sce_ddb[, keep]
detected <- rowSums(assay(sce_ddb, "counts") > 0)
sce_ddb <- sce_ddb[detected >= 10, ]

## 2. realize the filtered subset into memory
sce <- as(sce_ddb, "SingleCellExperiment")

## 3. standard in-memory pipeline
sce <- logNormCounts(sce)
dec <- modelGeneVar(sce)
sce <- runPCA(sce[getTopHVGs(dec, n = 2000), ])

## 4. persist results back to Parquet
writeParquet(sce, "data/processed")
```

Many `r Biocpkg("scran")`/`r Biocpkg("scuttle")` steps can run *directly* on the
DuckDB-backed object as SQL, before any realization, that is the subject of
*Benchmarking BiocDuckDB*.

Step 3's `runPCA()` is a partial exception to needing step 2 first:
`r Biocpkg("BiocSingular")`'s `runSVD()` (and so `scran::fixedPCA()` /
`scater::runPCA()`) already compute a correct PCA directly on `sce_ddb`, through
the SQL-pushdown `%*%`/`crossprod` methods `r Biocpkg("DuckDBArray")` implements,
no realization required for correctness. Speed is a separate question: the
ordinary `BSPARAM = BiocSingular::IrlbaParam()` calls `%*%` once per solver
iteration, so on a `DuckDBMatrix` that cost is paid repeatedly rather than once,
covered in detail (with measured numbers) in *Benchmarking BiocDuckDB*. Passing
`BSPARAM = DuckDBIrlbaParam()` instead closes that gap: it materializes just the
matrix `runPCA()` touches (already filtered and HVG-selected here, not the whole
`SingleCellExperiment` the way step 2's `as(sce_ddb, "SingleCellExperiment")`
would) and drives `irlba` directly on it, falling back to the ordinary path
whenever that materialization doesn't apply. It is a lighter-weight, PCA-specific
alternative to full realization, not just a slower substitute for it.

# Multi-omics

For a `MultiAssayExperiment`, write each experiment to its own directory and read
them back as DuckDB-backed experiments:

```{r mae, eval=FALSE}
library(MultiAssayExperiment)

writeParquet(rna_sce,     file.path(mae_path, "rna"))
writeParquet(protein_se,  file.path(mae_path, "protein"))

mae <- MultiAssayExperiment(experiments = list(
    rna     = readParquet(file.path(mae_path, "rna")),
    protein = readParquet(file.path(mae_path, "protein"))))
```

Each modality stays on disk until accessed. Spatial experiments
(`MultiAssaySpatialExperiment`) are supported too, with their geometry written as
GeoParquet via `r Biocpkg("DuckDBSpatial")`.

# When to use BiocDuckDB

A good fit for experiments **larger than memory**, for **exploration, QC, and
filtering** of large datasets, and for a **self-describing, cross-language** on-disk
format (the Parquet files are readable by Python, Julia, and DuckDB directly).
Parquet is also a compact archival format, typically several-fold smaller than
`.rds` for sparse data, and self-describing.

Realize a subset to an ordinary in-memory object for iterative algorithms with
random access patterns or dense linear algebra; and small datasets that fit
comfortably in memory need no backend at all.

# Session information

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