Genomic Intervals

Genomic intervals are a fundamental abstraction in bioinformatics, representing regions on a chromosome or contig. BioLang provides the Interval type as a first-class value with set-theoretic operations, an implicit interval tree for efficient overlap queries, and windowing functions for sliding-window analyses.

Creating Intervals

Use the interval() constructor to create individual intervals. Coordinates are zero-based, half-open (matching BED convention): 

# Basic interval: chromosome, start, end
let region = interval("chr1", 1000, 2000)

# With strand information (4th positional argument)
let exon = interval("chr1", 1000, 2000, "+")
let antisense = interval("chr1", 1000, 2000, "-")

# Access fields
print(region.chrom)    # => "chr1"
print(region.start)    # => 1000
print(region.end)      # => 2000
print(exon.strand)     # => "+"

# Width/length
region |> len()        # => 1000

Intersect

intersect() returns the overlapping portion of two intervals. If they do not overlap, it returns nil: 

# Intersect two BED tables (tables with chrom, start, end columns)
let exons = read_bed("data/regions.bed")
let peaks = read_bed("data/regions.bed")

# Find all pairwise intersections
let shared = intersect(exons, peaks)
print(nrow(shared), "overlapping regions")

Merge

merge() combines overlapping or book-ended intervals into a single contiguous interval. On collections, it merges all overlapping intervals: 

# Merge overlapping intervals in a table (collapses overlapping rows)
# The table must have chrom, start, end columns
let regions = read_bed("data/regions.bed")
let merged = merge_intervals(regions)
print(nrow(merged), "merged regions")

Subtract

subtract() removes the portion of one interval that overlaps with another. The result may be zero, one, or two intervals: 

# Subtract one table's intervals from another
let exons = read_bed("data/regions.bed")
let repeats = read_bed("data/regions.bed")

# Remove repeat-overlapping portions from exons
let clean_exons = subtract(exons, repeats)
print(nrow(clean_exons), "clean exon regions")

Closest

closest() finds the nearest interval in a collection to a query interval. Returns both the interval and the distance: 

# Find closest intervals between two tables
# Returns a table with the closest match for each row in the first table
let queries = read_bed("data/regions.bed")
let genes = read_bed("data/regions.bed")

let result = closest(queries, genes)
print(result)

Interval Tree

For repeated overlap queries against a large collection, build an interval tree for O(log n + k) query time instead of O(n): 

# Build an interval tree from a BED file
let tree = read_bed("data/regions.bed") |> interval_tree()

# Query the tree for overlapping intervals
let hits = query_overlaps(tree, "chr1", 50000, 51000)
print(nrow(hits), "annotations overlap the query region")

# Find nearest interval to a region
let nearest = query_nearest(tree, "chr1", 50000, 51000)
print(nearest)

Sliding Window Analysis

Use windows() to extract typed sliding windows of a given size for genome-wide analyses: 

# Sliding-window GC content
let seq = dna"ATCGATCGATCGATCGATCGATCG"
windows(seq, 6, 1)
  |> map(|window| {
    start: window.pos,
    gc: gc_content(window.seq)
  })
  |> to_table()
  |> print()

Combining Operations

Interval operations compose naturally with BioLang's pipe syntax: 

# Find promoter regions that overlap with ChIP-seq peaks
# but not with known repeat elements
let genes = read_bed("data/regions.bed")
let peaks = read_bed("data/regions.bed")
let repeats = read_bed("data/regions.bed")

# Use flank() to get upstream regions (2kb)
let promoters = flank(genes, 2000)

# Intersect with peaks, then subtract repeats, then merge
let result = intersect(promoters, peaks)
result = subtract(result, repeats)
result = merge_intervals(result)
write_bed(result, "clean_promoter_peaks.bed")

Table-Based BED Operations

The following functions operate on Tables that have chrom, start, and end columns — matching standard BED format. They implement the most-used bedtools operations natively, without shelling out to external tools. All coordinate arithmetic follows the zero-based, half-open convention.

bed_intersect(a, b)

For each row in table a, finds every overlapping row in table b (same chrom, overlapping coordinates). Returns a table containing all columns from a plus three derived columns:

  • overlap_start — start of the shared region
  • overlap_end — end of the shared region
  • overlap_length — number of overlapping bases

Rows in a with no overlap in b are dropped. To keep them, use bed_closest() with a distance filter instead. 

# Find CDS intervals that overlap a variant set
let cds    = read_bed("annotation/cds.bed")
let vcf    = read_bed("variants/snps.bed")   # bed3+ export from VCF

let coding_variants = bed_intersect(cds, vcf)

# coding_variants now has all cds columns plus overlap_start, overlap_end, overlap_length
print(nrow(coding_variants), "variants fall inside coding sequence")

# Filter to only variants with at least 1 bp overlap (always true for bed_intersect,
# but useful after a join that may produce zero-length artefacts)
let confirmed = coding_variants |> filter(|r| r.overlap_length > 0)

# Summarise per gene
confirmed
  |> group_by("gene_id")
  |> summarize(|gene_id, rows| {
       { gene_id: gene_id, variant_count: len(rows) }
     })
  |> sort_by(|row| -row.variant_count)
  |> print()

bed_subtract(a, b)

Removes all bases in b from the intervals in a. An interval in a that is only partially covered by b is split into one or two shorter intervals. The returned table has the same columns as a; extra columns from b are not included. 

# Mask repeat regions from a gene annotation before peak analysis
let genes   = read_bed("annotation/genes.bed")
let repeats = read_bed("annotation/repeatmasker.bed")

let clean_genes = bed_subtract(genes, repeats)

print(nrow(genes),       "input gene intervals")
print(nrow(clean_genes), "intervals after repeat masking")

# Compare total coverage before and after
let bp_before = genes       |> map(|r| r.end - r.start) |> sum()
let bp_after  = clean_genes |> map(|r| r.end - r.start) |> sum()

print("Repeat content removed:", bp_before - bp_after, "bp")

# Write masked annotation for downstream analysis
write_bed(clean_genes, "annotation/genes_masked.bed")

bed_merge(a, gap?)

Merges overlapping or adjacent intervals within a. The optional gap parameter (default 0) allows intervals separated by at most gap bases to be merged as well — useful for collapsing nearly-adjacent peaks into a single region. The algorithm sorts by (chrom, start) and then greedily extends intervals. The returned table has the same columns as a; non-coordinate columns are taken from the first row of each merged group. 

# Merge overlapping exons produced by an RNA-seq assembler
let raw_exons = read_bed("rnaseq/assembled_exons.bed")

# Strict merge: only overlapping / book-ended intervals
let merged_strict = bed_merge(raw_exons)
print(nrow(merged_strict), "exons after strict merge")

# Permissive merge: also join exons separated by <= 50 bp (e.g. small introns)
let merged_50 = bed_merge(raw_exons, 50)
print(nrow(merged_50), "exons after 50 bp gap merge")

# Typical ChIP-seq workflow: load peaks, merge, write
read_bed("chipseq/peaks.bed")
  |> bed_merge(0)
  |> write_bed("chipseq/peaks_merged.bed")

# Count how many raw peaks collapsed into each merged region
let with_counts =
  bed_intersect(merged_strict, raw_exons)
    |> group_by("chrom")
    |> summarize(|chrom, rows| {
         { chrom: chrom, source_count: len(rows) }
       })
print(with_counts)

bed_closest(a, b)

For each row in a, finds the single nearest row in b on the same chromosome. The distance column is 0 when the intervals overlap; otherwise it is the number of bases between the nearest endpoints. Columns from b are prefixed with b_ to avoid name collisions. 

# Assign each ChIP-seq peak to its nearest annotated gene
let peaks = read_bed("chipseq/peaks_merged.bed")
let genes = read_bed("annotation/genes.bed")   # must have chrom, start, end, gene_name

let peak_gene = bed_closest(peaks, genes)

# peak_gene columns: all peak cols + b_chrom, b_start, b_end, b_gene_name, distance

# Keep only peaks within 50 kb of a gene
let proximal = peak_gene |> filter(|r| r.distance <= 50000)
print(nrow(proximal), "peaks within 50 kb of a gene")

# Summarise: how many peaks per gene, and what is the median distance?
proximal
  |> group_by("b_gene_name")
  |> summarize(|gene_name, rows| {
       {
         gene_name: gene_name,
         peak_count: len(rows),
         median_distance: rows |> col("distance") |> median()
       }
     })
  |> sort_by(|row| -row.peak_count)
  |> head(20)
  |> print()

# Full pipeline: load → merge peaks → assign genes → filter → export
read_bed("chipseq/peaks.bed")
  |> bed_merge(100)
  |> bed_closest(genes)
  |> filter(|r| r.distance <= 10000)
  |> select("chrom", "start", "end", "b_gene_name", "distance")
  |> write_bed("chipseq/peak_gene_assignments.bed")