---
title: "SuperIntervals with Bioconductor: genomic overlap workflows"
output: rmarkdown::html_vignette
vignette: >
    %\VignetteIndexEntry{SuperIntervals with Bioconductor: genomic overlap workflows}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

```{r setup, message=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(GenomicRanges)
library(superintervalsr)
```

## Annotating peaks with genes

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.

```{r input}
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
```

## Retrieving source rows after indexing

`build()` 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:

```{r row_ids}
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]
```

To return new genomic ranges containing the selected metadata column:

```{r conversion}
hits <- search_items(gene_map, 5, 22)
as_granges(hits, seqname = "chr1", valuename = "gene_id")
```

The conversion assigns the requested sequence name. It does not restore original
strand, sequence lengths, or metadata columns omitted when constructing the map.

## Set operations for covered regions and nearby features

```{r operations}
peak_map <- IntervalMap(peaks, value_column = "gene_count")
shared <- superintervalsr::intersection(gene_map, peak_map)
as_granges(search_items(shared, 1, 50), seqname = "chr1")

# 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")

# Clamped search windows; results are already indexed.
windows <- superintervalsr::expand(peak_map, 5, 5, lo = 1, hi = 50)
search_items(windows, 1, 50)
```

Intersection 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:

```{r coverage}
superintervalsr::coverage(gene_map, 5, 10)
stopifnot(superintervalsr::coverage(gene_map, 5, 10)$total_coverage == 9)
```

## Repeated queries

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:

```{r batch}
individual <- vapply(seq_along(peaks), function(i) {
        count(gene_map, start(peaks)[i], end(peaks)[i])
}, integer(1))
stopifnot(identical(individual, counts))
```


## A real RNA-seq example: promoter overlaps in airway cells

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.

```{r airway-data, message=FALSE}
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.

```{r airway-query}
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)
```

The 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.

```{r airway-summary}
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, ])
```

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