knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(GenomicRanges)
library(superintervalsr)This small synthetic example demonstrates closed, end-inclusive
genomic coordinates. A map indexes one chromosome and ignores strand,
equivalent to ignore.strand = TRUE in the comparison below.
Zero-width GRanges are not supported. Use separate maps for distinct
chromosomes.
genes <- GRanges(c("chr1", "chr1", "chr1", "chr2"),
IRanges(c(20, 1, 8, 1), c(30, 10, 18, 10)),
gene_id = c("GENE_C", "GENE_A", "GENE_B", "GENE_D"))
peaks <- GRanges("chr1", IRanges(c(5, 15, 25, 40), c(8, 22, 28, 45)))
chr1_rows <- which(as.character(seqnames(genes)) == "chr1")
chr1_genes <- genes[chr1_rows]
gene_map <- IntervalMap(chr1_genes, value_column = "gene_id")
counts <- count.batch(gene_map, start(peaks), end(peaks))
reference <- countOverlaps(peaks, chr1_genes, ignore.strand = TRUE)
stopifnot(identical(counts, reference))
mcols(peaks)$gene_count <- counts
peaks
#> GRanges object with 4 ranges and 1 metadata column:
#> seqnames ranges strand | gene_count
#> <Rle> <IRanges> <Rle> | <integer>
#> [1] chr1 5-8 * | 2
#> [2] chr1 15-22 * | 2
#> [3] chr1 25-28 * | 1
#> [4] chr1 40-45 * | 0
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengthsbuild() sorts coordinates. Search indices refer to map
positions, so they must not be used directly to subset an unsorted or
chromosome-filtered input. Instead, retain original row IDs as
values:
row_map <- IntervalMap(start(chr1_genes), ends = end(chr1_genes),
values = chr1_rows)
row_ids <- unlist(search_values(row_map, 5, 22), use.names = FALSE)
genes[row_ids]
#> GRanges object with 3 ranges and 1 metadata column:
#> seqnames ranges strand | gene_id
#> <Rle> <IRanges> <Rle> | <character>
#> [1] chr1 20-30 * | GENE_C
#> [2] chr1 8-18 * | GENE_B
#> [3] chr1 1-10 * | GENE_A
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengthsTo return new genomic ranges containing the selected metadata column:
hits <- search_items(gene_map, 5, 22)
as_granges(hits, seqname = "chr1", valuename = "gene_id")
#> GRanges object with 3 ranges and 1 metadata column:
#> seqnames ranges strand | gene_id
#> <Rle> <IRanges> <Rle> | <list>
#> [1] chr1 20-30 * | GENE_C
#> [2] chr1 8-18 * | GENE_B
#> [3] chr1 1-10 * | GENE_A
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengthsThe conversion assigns the requested sequence name. It does not restore original strand, sequence lengths, or metadata columns omitted when constructing the map.
peak_map <- IntervalMap(peaks, value_column = "gene_count")
shared <- superintervalsr::intersection(gene_map, peak_map)
as_granges(search_items(shared, 1, 50), seqname = "chr1")
#> GRanges object with 5 ranges and 1 metadata column:
#> seqnames ranges strand | value
#> <Rle> <IRanges> <Rle> | <list>
#> [1] chr1 25-28 * | GENE_C
#> [2] chr1 20-22 * | GENE_C
#> [3] chr1 15-18 * | GENE_B
#> [4] chr1 8 * | GENE_B
#> [5] chr1 5-8 * | GENE_A
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths
# Disjoint covered regions, and uncovered regions within a specified span.
merged <- superintervalsr::merge_overlaps(gene_map)
uncovered <- superintervalsr::gaps(merged, 1, 50)
as_granges(search_items(uncovered, 1, 50), seqname = "chr1")
#> GRanges object with 2 ranges and 1 metadata column:
#> seqnames ranges strand | value
#> <Rle> <IRanges> <Rle> | <list>
#> [1] chr1 31-50 * |
#> [2] chr1 19 * |
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths
# Clamped search windows; results are already indexed.
windows <- superintervalsr::expand(peak_map, 5, 5, lo = 1, hi = 50)
search_items(windows, 1, 50)
#> $start
#> [1] 35 20 10 1
#>
#> $end
#> [1] 50 33 27 13
#>
#> $value
#> $value[[1]]
#> [1] 0
#>
#> $value[[2]]
#> [1] 1
#>
#> $value[[3]]
#> [1] 2
#>
#> $value[[4]]
#> [1] 2Intersection keeps the value of each left input interval and can emit overlapping pieces. Merge and union keep one representative value per merged region; they are coordinate operations, not an aggregation of all annotations. Adjacent intervals without shared bases remain separate. Coverage sums intersection widths with multiplicity:
Batch operations reduce R-to-C++ call overhead. Build once, query repeatedly, and rebuild after additional intervals are added. The query path deliberately does not check whether the index is up to date. Use valid whole-number query ranges with start less than or equal to end. Compare correctness before measuring performance on data representative of your analysis:
The airway experiment contains RNA-seq measurements from
human airway smooth muscle cells treated with dexamethasone. Its
SummarizedExperiment includes GRanges annotations and a count matrix.
This example asks which annotated gene bodies intersect promoter windows
for highly expressed genes. It preserves chromosome, strand, gene
identifiers and expression summaries throughout.
data("airway", package = "airway")
# The gene models are GRangesList exon annotations. Collapse each gene to its
# genomic span and keep models mapping to one contiguous sequence/strand group.
models <- SummarizedExperiment::rowRanges(airway)
gene_spans <- range(models)
single_span <- lengths(gene_spans) == 1L
features <- unlist(gene_spans[single_span], use.names = FALSE)
names(features) <- names(models)[single_span]
mcols(features)$mean_count <-
rowMeans(SummarizedExperiment::assay(airway))[single_span]
features <- features[width(features) > 0L]
# Retain a small, reproducible subset of real annotations for this demonstration.
features <- features[as.character(seqnames(features)) %in% c("1", "2")]
stopifnot(length(features) > 0L)
expressed <- features[mcols(features)$mean_count > 100]
expressed <- head(expressed, 100L)
promoter_windows <- promoters(expressed, upstream = 1000L, downstream = 100L)Build the subject index once and reuse it. Unlike the low-level IntervalMap, SuperIntervalMap accepts all chromosomes and handles strand matching. Hits address the original input rows, so annotations can be recovered directly.
genomic_index <- SuperIntervalMap(features)
hits <- findOverlaps(promoter_windows, genomic_index, ignore.strand = TRUE)
reference <- findOverlaps(promoter_windows, features, ignore.strand = TRUE)
# Compare the query/subject pairs; the order within each query is unspecified.
pairs <- function(x) sort(paste(queryHits(x), subjectHits(x)))
stopifnot(identical(pairs(hits), pairs(reference)))
stopifnot(identical(countOverlaps(promoter_windows, genomic_index),
countOverlaps(promoter_windows, features)))
annotated_features <- as_granges(genomic_index)[subjectHits(hits)]
head(annotated_features, 3)
#> GRanges object with 3 ranges and 1 metadata column:
#> seqnames ranges strand | mean_count
#> <Rle> <IRanges> <Rle> | <numeric>
#> ENSG00000000457 1 169818772-169863408 - | 242.00
#> ENSG00000000971 1 196621008-196716634 + | 6034.75
#> ENSG00000001461 1 24742284-24799466 + | 2828.88
#> -------
#> seqinfo: 722 sequences (1 circular) from an unspecified genomeThe result includes the source GRanges metadata, names, strand and
seqinfo. For a stranded analysis leave
ignore.strand = FALSE (the default); "*"
strands have the same wildcard behavior as GenomicRanges. Alternative
overlap types and thresholds delegate to GenomicRanges. Recreate the
index after changing the subject ranges, and serialize the source
GRanges rather than the external-pointer index.
annotation <- data.frame(
promoter_gene = names(expressed)[queryHits(hits)],
overlapping_gene = names(features)[subjectHits(hits)],
mean_count = mcols(annotated_features)$mean_count
)
# Exclude each promoter's own gene to inspect neighboring gene annotations.
head(annotation[annotation$promoter_gene != annotation$overlapping_gene, ])
#> promoter_gene overlapping_gene mean_count
#> 3 ENSG00000001460 ENSG00000001461 2828.875
#> 6 ENSG00000001461 ENSG00000001460 188.125
#> 12 ENSG00000003509 ENSG00000115816 530.000
#> 13 ENSG00000004455 ENSG00000142920 54.625
#> 18 ENSG00000005448 ENSG00000159239 95.750
#> 20 ENSG00000007341 ENSG00000116489 3159.875sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 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] GenomicRanges_1.65.4 Seqinfo_1.3.2 IRanges_2.47.5
#> [4] S4Vectors_0.51.10 BiocGenerics_0.59.12 generics_0.1.4
#> [7] superintervalsr_0.99.7 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-6 jsonlite_2.0.0
#> [3] compiler_4.6.1 Rcpp_1.1.2
#> [5] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [7] jquerylib_0.1.4 yaml_2.3.12
#> [9] fastmap_1.2.0 lattice_0.23-1
#> [11] R6_2.6.1 XVector_0.53.0
#> [13] S4Arrays_1.13.2 knitr_1.52
#> [15] DelayedArray_0.39.8 MatrixGenerics_1.25.0
#> [17] maketools_1.3.2 bslib_0.12.0
#> [19] rlang_1.3.0 cachem_1.1.0
#> [21] xfun_0.61 sass_0.4.10
#> [23] sys_3.4.3 otel_0.2.0
#> [25] SparseArray_1.13.4 cli_3.6.6
#> [27] digest_0.6.39 grid_4.6.1
#> [29] lifecycle_1.0.5 evaluate_1.0.5
#> [31] buildtools_1.0.0 abind_1.4-8
#> [33] matrixStats_1.5.0 tools_4.6.1
#> [35] htmltools_0.5.9