bio-genome-intervals-proximity-operations

What this skill does

## Version Compatibility

Reference examples tested with: bedtools 2.31+, pybedtools 0.10+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `bedtools --version` then `bedtools <subcommand> --help` to confirm flags
- Python: `pip show pybedtools` then `help(pybedtools.BedTool.closest)` to check signatures

`flank` and `slop` REQUIRE a chrom-sizes (`genome.txt`, two columns: chrom<TAB>length) file via `-g`; `closest` requires both inputs coordinate-sorted (`sort -k1,1 -k2,2n`). If code throws an error, introspect the installed tool and adapt rather than retrying.

# Proximity Operations

**"Which gene is nearest to each peak, and is that the gene it regulates?"** -> Compute interval geometry (nearest feature, signed distance, window membership, strand-aware promoters) with bedtools, then decide honestly whether geometry answers the biological question.
- CLI: `bedtools closest -D b -t first -a peaks.bed -b genes.bed`, `bedtools window -w 50000`, `bedtools slop -s -l 2000 -r 200 -g genome.txt`
- Python: `peaks.closest(genes.sort(), D='b', t='first')`, `peaks.window(genes, w=50000)`, `tss.slop(g='genome.txt', s=True, l=2000, r=200)` (pybedtools)

## The Single Most Important Modern Insight -- closest Answers a GEOMETRY Question Misread as a BIOLOGY Question

`bedtools closest` answers "what is the nearest annotated TSS?" - a coordinate fact. The user almost always wants "which gene does this element regulate?" - a biology claim. For distal regulatory elements these disagree the **majority of the time**. In the CRISPRi-FlowFISH gold standard (Fulco 2019 *Nat Genet* 51:1664), assigning each tested distal element to the **closest expressed gene gave only ~47% precision and ~37% recall** - the nearest gene was the wrong target most of the time, and the method missed nearly two-thirds of real links. Enhancers routinely skip intervening genes: the canonical case is the obesity-associated *FTO* intron regulating **IRX3 ~500 kb away**, not *FTO* (Smemo 2014 *Nature* 507:371). Do the bedtools arithmetic flawlessly here, then route real enhancer->gene linking to activity/contact/QTL methods (ABC: Fulco 2019, Nasser 2021; PCHi-C; eQTL-coloc) at atac-seq/enhancer-gene-linking - never present "nearest gene" as a regulatory call for a distal element.

The deeper twist - **two regimes, opposite advice, identical command:**
- **Enhancer -> target (closest is a TRAP).** Distal ATAC/H3K27ac peaks, enhancer GWAS variants: nearest gene is wrong most of the time. Use as a candidate generator, validate with ABC/PCHi-C/eQTL.
- **GWAS locus -> gene (closest is a fair PRIOR).** For a fine-mapped, colocalized credible-set SNP, the nearest **protein-coding** gene is right ~50-65% of the time - a strong, hard-to-beat baseline for which gene a locus implicates. Route to causal-genomics for the rigorous version, but nearest-coding-gene is a defensible first pass.

Conflating the two regimes is the real error. The discriminator: is the question the *target of an enhancer* (distrust nearest) or *the gene under a GWAS peak* (nearest is a fine first pass)?

## Operation Taxonomy

| Operation | What it computes | Strand-aware? | Needs genome file? |
|-----------|------------------|---------------|--------------------|
| closest | For each A, the nearest B (+ optional signed distance) | optional (`-s`/`-S`, `-D a`/`-D b`) | no |
| window | For each A, all B within +-W bp (fuzzy intersect) | optional (`-sw`/`-sm`/`-Sm`) | no |
| slop | Grow each interval by N bp, keeping it one feature | optional (`-s`) | yes (`-g`) |
| flank | Emit the regions BESIDE each interval, dropping the body | optional (`-s`) | yes (`-g`) |

`closest`/`window` are queries (A vs B); `slop`/`flank` are transforms (A only, + genome file). The slop-vs-flank distinction trips people: `slop -b 1000` makes a peak 2 kb wider (one feature); `flank -b 1000` returns only the left/right neighboring 1 kb regions and discards the peak itself (two features). `window -w 0` is approximately `intersect`.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Nearest gene to a promoter-proximal mark (H3K4me3, Pol II, CAGE) | `closest -D b -io -t first` | the peak really is at the gene it marks; closest is honest here |
| Distal enhancer / ATAC peak -> which gene? | `closest`/`window` as candidates, then -> atac-seq/enhancer-gene-linking | nearest is wrong the majority of the time (ABC/PCHi-C/eQTL link it) |
| GWAS credible-set SNP -> implicated gene | `closest` to nearest **protein-coding** gene, then -> causal-genomics/colocalization-analysis | nearest-coding-gene is a ~50-65% prior; a fair first pass |
| All candidate genes near an element | `window -w 50000` (or TAD-scale) | honest "candidate set", not a single call |
| Build promoters from a gene model | collapse to TSS, then `slop -s -l UP -r DOWN -g` | a promoter is an imposed definition, strand-aware, from the TSS |
| Distance-to-TSS distribution | `closest -D b -d` then plot signed distance | a distribution beats a binary "promoter vs distal" threshold |
| Upstream-only / downstream-only nearest | `closest -D b -iu` / `-id` | direction must be strand-relative (`-D b`), never `-D ref` |
| Peak-set GO enrichment from proximity | -> GREAT/rGREAT (regulatory-domain model) | avoids the `-t all` double-counting and distal mis-assignment |
| Regions flanking a feature (splice/boundary context) | `flank -s -b N -g` | the regions outside the feature, strand-aware |
| Peaks not yet called | -> chip-seq/peak-calling, atac-seq/atac-peak-calling | this skill operates on existing intervals |

## closest - Nearest Feature with Signed, Strand-Aware Distance

Default: for each A, report the single nearest B; **on ties, report ALL tied B** (`-t all` is the default - the double-counting trap below). Both inputs must be sorted. When A's chromosome has no B feature, bedtools prints `none` for B columns and **`-1`** for distance - filter this sentinel before any numeric summary.

```bash
# Nearest gene, signed distance by the GENE's strand, ignore overlaps, one row per peak
bedtools sort -i peaks.bed > peaks.sorted.bed
bedtools sort -i genes.bed  > genes.sorted.bed
bedtools closest -a peaks.sorted.bed -b genes.sorted.bed -D b -io -t first > nearest.bed
#                                                         ^^^^ sign by gene strand (biology, not coordinates)
#                                                              ^^^ closest non-overlapping gene
#                                                                  ^^^^^^^^ resolve ties deterministically (document this)

# k=3 nearest with unsigned distance (k>1 intentionally multiplies rows)
bedtools closest -a peaks.sorted.bed -b genes.sorted.bed -k 3 -d > top3.bed
```

```python
import pybedtools

peaks = pybedtools.BedTool('peaks.bed').sort()
genes = pybedtools.BedTool('genes.bed').sort()
near = peaks.closest(genes, D='b', io=True, t='first')      # -D b -io -t first
near = near.filter(lambda x: int(x.fields[-1]) != -1)        # drop the no-feature sentinel
near.saveas('nearest.bed')
```

Key flags: `-d` unsigned distance (overlaps = 0); `-D ref` signed by coordinate only (strand-agnostic - see Failure Modes); `-D a`/`-D b` signed by A's / B's strand; `-t all|first|last`; `-k N` k-nearest; `-io` ignore overlapping B; `-iu`/`-id` ignore upstream/downstream (require `-D`); `-fu`/`-fd` first upstream/downstream; `-s`/`-S` same/opposite strand; `-N` require different names; `-mdb each|all` and `-names`/`-filenames` for multiple `-b` files.

## window - Features Within a Search Radius

`window` reports all B within a window around each A (default 1000 bp each side). Use it for the honest "candidate genes near this element" framing.

```bash
# All genes within 50 kb of each peak, counted per peak
bedtools window -a peaks.bed -b genes.bed -w 50000 -c > peak_gene_counts.bed
```

Flags: `-w N` symmetric (default 1000); `-l N`/`-r N` asymmetric (coordinate left/