bio-genome-intervals-interval-arithmetic

What this skill does

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

| Tool | Role | Mechanism | When |
|------|------|-----------|------|
| bedtools | CLI interval algebra (reference implementation) | streaming sweep on sorted input; in-memory tree otherwise | shell pipelines, large files, reproducible one-liners |
| pybedtools | Python wrapper over bedtools | shells out to the bedtools binary; BedTool objects, flags as kwargs | inside a Python script; need exact bedtools parity; chaining with pandas |
| pyranges | pure-Python vectorized engine | native NumPy/pandas PyRanges; no bedtools dependency | large in-memory joins, no bedtools install, dataframe-native; mind the v0/v1 split |
| bioframe | functions on a plain pandas DataFrame | vectorized pandas merges; columns `chrom/start/end` | data already in pandas / the cooler-Hi-C ecosystem |

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

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Quick overlap on the command line | `bedtools intersect -u` | no Python overhead; reproducible one-liner |
| Inside a pandas/Python pipeline | pybedtools or pyranges/bioframe | stays in-process; pyranges/bioframe need no bedtools binary |
| Whole-genome-scale intersect | `intersect -sorted -g genome.txt` | low-memory sweep; modern bedtools errors on a sort/order mismatch, `-g` pins the expected chromosome order |
| Spliced reads / BED12 vs exons | add `-split` | otherwise the intron-spanning envelope is intersected (RNA-seq inflation) |
| Are two SV/CNV calls the same event? | `-f 0.5 -r` (50% reciprocal) | one-sided fractions let a giant interval swallow a tiny one |
| Transfer/aggregate B values onto A | `bedtools map -c COL -o OP` | columnar alternative to `intersect -wo \| groupby` |
| Build consensus peakset from replicates | `cat \| sort \| merge -d N` | collapses replicate peaks within N bp |
| Multi-sample shared-region map | `multiinter` / `unionbedg` | presence/absence (intervals) or stacked signal matrix |
| Is the overlap more than chance? | -> overlap-significance | raw count is length/coverage-confounded; needs a permutation null |
| Peaks not yet called | -> chip-seq/peak-calling or atac-seq/atac-peak-calling | this category operates on existing intervals |

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

```bash
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)
```

```python
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

```bash
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