Bio Genome Intervals Interval Arithmetic

Version Compatibility

Reference examples tested with: bedtools 2.31+, pybedtools 0.10+, pyranges 0.1+ (or 1.0+ - see note), bioframe 0.7+.

Before using code patterns, verify installed versions match. If versions differ:

• CLI: bedtools --version then bedtools <subcommand> --help to confirm flags
• Python: pip show <package> then help(module.function) to check signatures

pyranges has a major-version API split: pyranges 0.x and the 1.0 rewrite (package pyranges1) differ in method names and return shapes. Verify with import pyranges; pyranges.__version__ before pasting v0 idioms. If code throws an error, introspect the installed package and adapt rather than retrying.

Interval Arithmetic

"Which of my peaks overlap promoters, and how do I combine/subtract/annotate interval sets?" -> Apply exact, deterministic set operations to sorted interval files, guarding the preconditions (prior sort, the -sorted chromosome-order contract, -split) that otherwise corrupt the answer.

• CLI: bedtools intersect -a a.bed -b b.bed -u, bedtools merge, bedtools subtract, bedtools map -c 4 -o mean
• Python: a.intersect(b, u=True), a.merge() (pybedtools); pr_a.overlap(pr_b) (pyranges); bf.overlap(df1, df2) (bioframe)

The Single Most Important Modern Insight -- The Arithmetic Is Exact; the Danger Is the Silent Preconditions

The set operations themselves are exact and deterministic - bedtools, pyranges, and bioframe compute identical geometry on the same 0-based half-open intervals. The bugs are never in the arithmetic; they hide in four preconditions that fail quietly, returning a plausible wrong answer with exit code 0:

1. merge, map, closest, groupby require prior sort. merge only collapses records that are adjacent in file order - on unsorted input, overlapping intervals survive un-merged and downstream counts are wrong, with no warning.
2. -sorted requires sorted input in a shared chromosome order. It swaps intersect's in-memory interval tree for a low-memory chromosome sweep. Modern bedtools (>=~2.25) detects unsorted or differently-ordered -sorted input and errors out (exit 1: ... is not sorted / chromomsome sort ordering ... is inconsistent); older versions silently swept past overlaps and under-reported. Pass -g genome.txt to pin the expected chromosome order (reproducible, and it catches the subtler missing-chromosome cases). The mismatch that stays SILENT on every version is a chromosome-NAME difference (chr1 vs 1), which returns an empty result with no error.
3. -split changes whether the count is exons or the spanning envelope. A BED12 record or spliced BAM read (CIGAR N) spans introns; without -split bedtools intersects the whole intron-spanning envelope, silently inflating RNA-seq overlaps. With -split it intersects only the blocks (exons).
4. A raw overlap count is not association. Long features, clustered features, and uneven coverage all inflate it; the number means nothing without a null. bedtools fisher is a weak analytic screen, not the answer - route rigorous significance to overlap-significance.

Tool Taxonomy

All three Python engines compute the same overlaps; the porting bugs are about default strand handling and return shape (pyranges overlap vs join vs intersect; bioframe overlap with how=), not geometry.

Decision Tree by Scenario

Intersect - the Workhorse

The output-mode flags do not change what overlaps; they change what gets printed (the #1 source of confusion). Full flag semantics are in usage-guide.md.

bedtools intersect -a peaks.bed -b genes.bed -u            # whole A, once, if it overlaps >=1 B
bedtools intersect -a peaks.bed -b genes.bed -v            # A features with NO overlap (set difference)
bedtools intersect -a peaks.bed -b genes.bed -c            # per-A count of B hits (0 if none)
bedtools intersect -a peaks.bed -b genes.bed -wa -wb       # whole A + whole B, one line per pair ("join")
bedtools intersect -a peaks.bed -b genes.bed -loj          # left outer join: every A, NULL B if none
bedtools intersect -a peaks.bed -b genes.bed -wo           # A+B+bp-of-overlap, only A with overlap
bedtools intersect -a peaks.bed -b genes.bed -wao          # like -wo but A-with-no-overlap kept (B=., overlap=0)

import pybedtools

a = pybedtools.BedTool('peaks.bed')
b = pybedtools.BedTool('genes.bed')
a.intersect(b, u=True)            # flags become kwargs
a.intersect(b, wa=True, wb=True)
a.intersect(b, c=True)

Subtract, Merge, Complement, Cluster

bedtools subtract -a a.bed -b b.bed             # clip the overlapping portions out of A (A can fragment)
bedtools subtract -a a.bed -b b.bed -A          # drop the ENTIRE A feature if any part overlaps B
bedtools sort -i a.bed | bedtools merge -d 0    # collapse overlapping + book-ended; -d 0 is the default
bedtools sort -i a.bed | bedtools merge -c 4,5 -o distinct,sum   # summarize columns while merging
bedtools complement -i a.bed -g genome.txt      # the gaps: genome NOT covered by A (genome file required)
bedtools sort -i a.bed | bedtools cluster -d 0  # assign a cluster id to overlapping/adjacent features

Valid -o operations: sum, min, max, absmin, absmax, mean, median, mode, antimode, stdev, sstdev, collapse, distinct, count, count_distinct, first, last. merge -d 0 merges overlapping and book-ended (touching) features but NOT a 1 bp gap; -d 1 does.

Map - Transfer Values, and Groupby - Aggregate

Goal: Summarize a column of overlapping B features onto each A interval (e.g. mean signal per gene).

Approach: For each sorted A interval, map collects overlapping B features and applies an aggregation -o to a B column -c; groupby is the single-file SQL-style aggregator after an intersect -wo.

bedtools map -a genes.bed -b scores.bedgraph -c 4 -o mean      # both inputs MUST be sorted
bedtools intersect -a genes.bed -b peaks.bed -wo \
  | bedtools groupby -g 1,2,3,4 -c 13 -o sum                   # group on A cols, sum the overlap-bp col

import pybedtools

genes = pybedtools.BedTool('genes.bed').sort()
scores = pybedtools.BedTool('scores.bedgraph').sort()
genes.map(scores, c=4, o='mean')

Multi-Sample: Multiinter and Unionbedg

bedtools multiinter -header -names s1 s2 s3 -i s1.bed s2.bed s3.bed   # which files cover each sub-interval
bedtools unionbedg -header -names s1 s2 s3 -i s1.bg s2.bg s3.bg       # stack bedGraph signal into a matrix

multiinter is the interval presence/absence map (build a consensus by filtering its num/list columns); unionbedg is its signal-track analog.

Jaccard and Fisher - Mechanics Only

jaccard is a single similarity scalar |A n B| / |A u B| in [0,1], useful for all-vs-all dataset clustering - it is NOT a significance test (no p-value). fisher builds a 2x2 table and returns a Fisher p, but it estimates the in-neither cell from a mean-interval-size/genome-size heuristic, ignores genome structure, and is prone to inflation - treat it as a fast triage screen only.

bedtools jaccard -a a.bed -b b.bed -g genome.txt        # both sorted; reports jaccard + n_intersections
bedtools fisher  -a a.bed -b b.bed -g genome.txt         # weak analytic null; validate any low p by simulation

For a defensible enrichment p-value (size-preserving permutation in an accessible workspace, GAT/regioneR/LOLA/GREAT), route to overlap-significance.