DuckDBGRanges 0.99.4
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.if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("DuckDBGRanges")library(DuckDBGRanges)
library(GenomicRanges)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:
length(gr_ddb)
#> [1] 5
seqnames(gr_ddb)
#> DuckDBColumn of length 5
#> range1 range2 range3 range4 range5
#> chr1 chr1 chr1 chr2 chr2All 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
#> range3 chr1 300 360 61 + | GENE3
#> range2 chr1 200 275 76 - | GENE2
#> score
#> <double>
#> range1 10
#> range3 15
#> range2 20
#> -------
#> 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>
#> range5 chr2 250 295 46 + | GENE5
#> range2 chr1 200 275 76 - | GENE2
#> range4 chr2 150 230 81 - | GENE4
#> score
#> <double>
#> range5 18
#> range2 20
#> range4 25
#> -------
#> 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 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 DuckDBGRanges_0.99.4 GenomicRanges_1.65.1
#> [4] Seqinfo_1.3.0 DuckDBDataFrame_0.99.20 IRanges_2.47.2
#> [7] S4Vectors_0.51.6 BiocGenerics_0.59.11 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 bookdown_0.47 blob_1.3.0
#> [10] fastmap_1.2.0 jsonlite_2.0.0 Matrix_1.7-6
#> [13] DBI_1.3.0 BiocManager_1.30.27 purrr_1.2.2
#> [16] jquerylib_0.1.4 abind_1.4-8 duckdb_1.5.5
#> [19] cli_3.6.6 rlang_1.3.0 dbplyr_2.6.0
#> [22] XVector_0.53.0 withr_3.0.3 cachem_1.1.0
#> [25] DelayedArray_0.39.4 yaml_2.3.12 otel_0.2.0
#> [28] S4Arrays_1.13.0 tools_4.6.1 dplyr_1.2.1
#> [31] assertthat_0.2.1 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 MatrixGenerics_1.25.0 knitr_1.51
#> [46] htmltools_0.5.9 rmarkdown_2.31 compiler_4.6.1