Genomic workflows routinely produce interval sets that are awkward to hold in memory: a gnomAD-scale variant catalogue is hundreds of millions of positions, a multi-sample scATAC-seq study is millions of peaks, a comprehensive gene model is millions of exons. Loading the whole set into an R session, alongside the intermediate objects an analysis creates, is wasteful when a workflow only ever touches a region, a chromosome, or the ranges passing a metadata filter.
GenomicRanges
is the standard Bioconductor representation for such data, and its API
(seqnames(), start(), restrict(),
findOverlaps()) is what most genomic code is written
against. What varies is where the data lives: an ordinary
GRanges holds every range in memory.
DuckDBGRanges
provides a GRanges (and GRangesList) backed by
columnar Parquet queried through DuckDB. It exposes the full GenomicRanges
API, but the coordinates stay on disk and operations
are recorded as lazy SQL queries: you can filter by region, restrict,
subset, and compose pipelines without loading the ranges into
memory. The R object stays a small, constant size (~175 KB)
whether it fronts a thousand ranges or a hundred million.
It builds on DuckDBDataFrame
(the tabular foundation of the BiocDuckDB suite) and is
designed for the filter-first shape of most genomic work: use
DuckDB to cut a large interval set down to the ranges of interest, then
materialize that manageable subset to an ordinary GRanges
for the interval-tree operations R does best.
This vignette is a practical introduction. Two companion vignettes go further:
GRanges, on
scATAC-seq and variant-scale data.A DuckDBGRanges is backed by a Parquet file with columns
for the genomic coordinates. For a small example we build a
GRanges, write it to Parquet with arrow, and open
it as a DuckDBGRanges.
library(arrow)
gr <- GRanges(
seqnames = rep(c("chr1", "chr2"), c(3, 2)),
ranges = IRanges(start = c(100, 200, 300, 150, 250),
end = c(150, 275, 360, 230, 295)),
strand = c("+", "-", "+", "-", "+"),
gene_id = paste0("GENE", 1:5),
score = c(10, 20, 15, 25, 18))
gr_df <- as.data.frame(gr)
gr_df$range_id <- paste0("range", seq_len(nrow(gr_df)))
gr_path <- tempfile(fileext = ".parquet")
write_parquet(gr_df, gr_path)We construct the object by naming the coordinate columns; extra
columns become metadata columns (mcols). The
keycol argument names the column that identifies each
range.
gr_ddb <- DuckDBGRanges(gr_path,
seqnames = "seqnames", start = "start",
end = "end", strand = "strand",
keycol = list(range_id = gr_df$range_id),
mcols = c("gene_id", "score"))
gr_ddb
#> DuckDBGRanges object with 5 ranges and 2 metadata columns:
#> seqnames start end width strand | gene_id
#> <character> <integer> <integer> <integer> <character> | <character>
#> range1 chr1 100 150 51 + | GENE1
#> range2 chr1 200 275 76 - | GENE2
#> range3 chr1 300 360 61 + | GENE3
#> range4 chr2 150 230 81 - | GENE4
#> range5 chr2 250 295 46 + | GENE5
#> score
#> <double>
#> range1 10
#> range2 20
#> range3 15
#> range4 25
#> range5 18
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengthsIt looks and behaves like a GRanges, but the ranges stay
on disk:
All the standard GRanges accessors work and return the
same kinds of values, queried from disk on demand:
start(gr_ddb)
#> DuckDBColumn of length 5
#> range1 range2 range3 range4 range5
#> 100 200 300 150 250
width(gr_ddb)
#> DuckDBColumn of length 5
#> range1 range2 range3 range4 range5
#> 51 76 61 81 46
ranges(gr_ddb)
#> DuckDBDataFrame with 5 rows and 3 columns
#> start end width
#> <integer> <integer> <integer>
#> range1 100 150 51
#> range2 200 275 76
#> range3 300 360 61
#> range4 150 230 81
#> range5 250 295 46
gr_ddb$gene_id
#> DuckDBColumn of length 5
#> range1 range2 range3 range4 range5
#> GENE1 GENE2 GENE3 GENE4 GENE5Ranges can be selected by index, by a logical condition on a
coordinate or metadata column (which becomes a SQL WHERE
clause), or by overlap with a query region:
gr_ddb[1:3]
#> DuckDBGRanges object with 3 ranges and 2 metadata columns:
#> seqnames start end width strand | gene_id
#> <character> <integer> <integer> <integer> <character> | <character>
#> range1 chr1 100 150 51 + | GENE1
#> range2 chr1 200 275 76 - | GENE2
#> range3 chr1 300 360 61 + | GENE3
#> score
#> <double>
#> range1 10
#> range2 20
#> range3 15
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths
gr_ddb[seqnames(gr_ddb) == "chr1"]
#> DuckDBGRanges object with 3 ranges and 2 metadata columns:
#> seqnames start end width strand | gene_id
#> <character> <integer> <integer> <integer> <character> | <character>
#> range1 chr1 100 150 51 + | GENE1
#> range2 chr1 200 275 76 - | GENE2
#> range3 chr1 300 360 61 + | GENE3
#> score
#> <double>
#> range1 10
#> range2 20
#> range3 15
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths
gr_ddb[gr_ddb$score > 15]
#> DuckDBGRanges object with 3 ranges and 2 metadata columns:
#> seqnames start end width strand | gene_id
#> <character> <integer> <integer> <integer> <character> | <character>
#> range4 chr2 150 230 81 - | GENE4
#> range2 chr1 200 275 76 - | GENE2
#> range5 chr2 250 295 46 + | GENE5
#> score
#> <double>
#> range4 25
#> range2 20
#> range5 18
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths
gr_ddb[gr_ddb %over% GRanges("chr1:150-250")]
#> DuckDBGRanges object with 2 ranges and 2 metadata columns:
#> seqnames start end width strand | gene_id
#> <character> <integer> <integer> <integer> <character> | <character>
#> range1 chr1 100 150 51 + | GENE1
#> range2 chr1 200 275 76 - | GENE2
#> score
#> <double>
#> range1 10
#> range2 20
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengthsThe intra- and inter-range operations of GenomicRanges are implemented as SQL, so the familiar functions work directly:
shift(gr_ddb, 50L)
#> DuckDBGRanges object with 5 ranges and 2 metadata columns:
#> seqnames start end width strand | gene_id
#> <character> <integer> <integer> <integer> <character> | <character>
#> range1 chr1 150 200 51 + | GENE1
#> range2 chr1 250 325 76 - | GENE2
#> range3 chr1 350 410 61 + | GENE3
#> range4 chr2 200 280 81 - | GENE4
#> range5 chr2 300 345 46 + | GENE5
#> score
#> <double>
#> range1 10
#> range2 20
#> range3 15
#> range4 25
#> range5 18
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths
restrict(gr_ddb, start = 200L, end = 300L)
#> DuckDBGRanges object with 4 ranges and 0 metadata columns:
#> seqnames start end width strand
#> <character> <integer> <integer> <integer> <character>
#> 1 chr1 200 275 76 -
#> 2 chr1 300 300 1 +
#> 3 chr2 200 230 31 -
#> 4 chr2 250 295 46 +
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths
reduce(gr_ddb)
#> DuckDBGRanges object with 5 ranges and 0 metadata columns:
#> seqnames start end width strand
#> <character> <integer> <integer> <integer> <character>
#> 1 chr1 100 150 51 +
#> 2 chr1 300 360 61 +
#> 3 chr1 200 275 76 -
#> 4 chr2 250 295 46 +
#> 5 chr2 150 230 81 -
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengthscoverage() computes per-base depth the same way,
entirely in SQL, and returns an ordinary (small) RleList,
one run-length-encoded vector per seqname:
The idiomatic workflow uses DuckDB to shrink a large set to the
ranges of interest, then coerces to an ordinary GRanges for
the interval-tree operations (findOverlaps(),
nearest()) that in-memory algorithms do best:
## lazy, on-disk filter down to a manageable subset...
sub <- gr_ddb[seqnames(gr_ddb) == "chr1" & gr_ddb$score > 12]
## ...then materialize the small result for overlap work
sub_gr <- as(sub, "GRanges")
class(sub_gr)
#> [1] "GRanges"
#> attr(,"package")
#> [1] "GenomicRanges"
findOverlaps(sub_gr, GRanges("chr1:100-320"))
#> Hits object with 2 hits and 0 metadata columns:
#> queryHits subjectHits
#> <integer> <integer>
#> [1] 1 1
#> [2] 2 1
#> -------
#> queryLength: 2 / subjectLength: 1Grouped features (a transcript’s exons, a gene’s peaks) are held by a
DuckDBGRangesList, which stores each group’s coordinates as
DuckDB LIST[] columns, one row per group, and behaves like
a GRangesList:
grl_data <- data.frame(
transcript_id = c("ENST001", "ENST002"),
gene_id = c("GENE1", "GENE2"),
seqnames = I(list(rep("chr1", 3), rep("chr2", 2))),
start = I(list(c(100L, 200L, 300L), c(150L, 250L))),
width = I(list(c(50L, 75L, 60L), c(80L, 45L))),
strand = I(list(c("+", "-", "+"), c("-", "+"))))
grl_path <- tempfile(fileext = ".parquet")
write_parquet(grl_data, grl_path)
grl_ddb <- DuckDBGRangesList(grl_path,
seqnames = "seqnames", start = "start",
width = "width", strand = "strand",
mcols = "gene_id",
keycol = list(transcript_id = grl_data$transcript_id))
elementNROWS(grl_ddb)
#> ENST001 ENST002
#> 3 2
grl_ddb[[1]]
#> DuckDBGRanges object with 3 ranges and 0 metadata columns:
#> seqnames start end width strand
#> <character> <integer> <integer> <integer> <character>
#> 1 chr1 100 149 50 +
#> 2 chr1 200 274 75 -
#> 3 chr1 300 359 60 +
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengthsA good fit when the interval set is larger than
memory (or you want to keep memory free), when the workload is
filter-heavy, restricting to regions, chromosomes, or
metadata thresholds before the expensive step, or when the data already
lives on disk as Parquet that other tools (Python,
Julia, cloud query engines) can read. Because the R object is a constant
~175 KB, a DuckDBGRanges also lets you reference a
huge annotation without paying for it until a query touches it.
An in-memory GRanges remains preferable when the data
fits comfortably in RAM, and for overlap enumeration
(findOverlaps(), subsetByOverlaps()), where
the in-memory interval trees are hard to beat. (Filtering and
distanceToNearest(), by contrast, are faster on the DuckDB
backend, see the benchmarking vignette.) The two compose well: the
filter, then materialize pattern above uses each where it is
strongest.
For a quantitative comparison see Benchmarking DuckDBGRanges; for the class structure and SQL translation see Design and extension of DuckDBGRanges.
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 DuckDBGRanges_0.99.5 GenomicRanges_1.65.1
#> [4] Seqinfo_1.3.0 DuckDBDataFrame_0.99.23 IRanges_2.47.2
#> [7] S4Vectors_0.51.6 BiocGenerics_0.59.12 generics_0.1.4
#> [10] bit64_4.8.2 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.10 SparseArray_1.13.2 lattice_0.23-1
#> [4] digest_0.6.39 magrittr_2.0.5 evaluate_1.0.5
#> [7] grid_4.6.1 blob_1.3.0 fastmap_1.2.0
#> [10] jsonlite_2.0.0 Matrix_1.7-6 DBI_1.3.0
#> [13] BiocManager_1.30.27 purrr_1.2.2 jquerylib_0.1.4
#> [16] duckdb_1.5.5 abind_1.4-8 cli_3.6.6
#> [19] rlang_1.3.0 dbplyr_2.6.0 XVector_0.53.0
#> [22] withr_3.0.3 cachem_1.1.0 DelayedArray_0.39.6
#> [25] yaml_2.3.12 otel_0.2.0 S4Arrays_1.13.0
#> [28] tools_4.6.1 dplyr_1.2.1 assertthat_0.2.1
#> [31] buildtools_1.0.0 vctrs_0.7.3 R6_2.6.1
#> [34] matrixStats_1.5.0 lifecycle_1.0.5 bit_4.6.0
#> [37] pkgconfig_2.0.3 bslib_0.12.0 pillar_1.11.1
#> [40] glue_1.8.1 xfun_0.60 tibble_3.3.1
#> [43] tidyselect_1.2.1 sys_3.4.3 MatrixGenerics_1.25.0
#> [46] knitr_1.51 htmltools_0.5.9 rmarkdown_2.31
#> [49] maketools_1.3.2 compiler_4.6.1