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>
#> range2 chr1 200 275 76 - | GENE2
#> range3 chr1 300 360 61 + | GENE3
#> range1 chr1 100 150 51 + | GENE1
#> score
#> <double>
#> range2 20
#> range3 15
#> range1 10
#> -------
#> 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 seqlengthsThe 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.1 GenomicRanges_1.65.1
#> [4] Seqinfo_1.3.0 DuckDBDataFrame_0.99.6 IRanges_2.47.2
#> [7] S4Vectors_0.51.5 BiocGenerics_0.59.10 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.22-9
#> [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-5 DBI_1.3.0
#> [13] BiocManager_1.30.27 purrr_1.2.2 jquerylib_0.1.4
#> [16] duckdb_1.5.4.3 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.3
#> [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.11.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