bio-genome-intervals-overlap-significance

What this skill does

## Version Compatibility

Reference examples tested with: bedtools 2.31+, pybedtools 0.10+, regioneR 1.36+ (Bioconductor 3.18+), GAT 1.3+, LOLA 1.30+, rGREAT 2.4+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('regioneR')` then `?permTest` to verify parameters

Bioconductor packages (regioneR, LOLA, rGREAT) are version-pinned to the Bioconductor release, not just the package version - record the Bioconductor release with results. rGREAT 2.x runs GREAT locally with general background handling; rGREAT 1.x only proxied the (whole-genome-default) web server. If code throws an error, introspect the installed package and adapt the example to match the actual API rather than retrying.

# Overlap Significance

**"My peaks overlap enhancers a lot - is that more than chance?"** -> Compare the observed overlap to a null distribution from a structured-genome model, not to a uniform-random expectation, and report a permutation p-value/z-score, not a raw count.
- CLI: `bedtools fisher -a A -b B -g genome.txt` (fast screen); `gat-run.py --segments=A --annotations=B --workspace=accessible.bed --isochores=gc.bed --num-samples=10000`
- Python: `BedTool(A).shuffle(g='genome.txt', excl='blacklist.bed').jaccard(B)` looped for a null (pybedtools)
- R: `permTest(A=peaks, B=enhancers, randomize.function=circularRandomizeRegions, evaluate.function=numOverlaps, genome='hg38', mask=blacklist, ntimes=1000)` (regioneR)

## The Single Most Important Modern Insight -- A Raw Overlap Count Is an Observation in Search of a Null, and the Universe Choice Dominates

"847 of 1,000 peaks overlap enhancers" means nothing until the analysis can say what that number would have been *by chance* - and "by chance" is almost never "place the regions uniformly at random on the genome." The genome is structured: genes cluster, GC varies in megabase isochores, mappability is uneven, half the genome is repeat/gap, and query regions are drawn from a *biased universe* (open chromatin, callable space, exons). Two tracks that share nothing but a gene-rich, high-GC, high-mappability habitat overlap far more than uniform-random expectation, and a naive test returns p < 1e-300. The co-localization is real; the *interpretation* ("functional association") is false. Three load-bearing moves:

1. **The universe/background is the lever that moves the whole answer - bigger than the test choice.** Across LOLA (`userUniverse`), GREAT (background regions), GAT (`--workspace`), and regioneR (the `mask`/`resampleRegions` universe), the most consequential choice is *the set of regions the query could have come from*. ATAC/ChIP peaks can only be called in accessible chromatin; their honest universe is "all accessible regions," not the genome. Testing against the whole genome merely rediscovers that open chromatin is gene-rich - every gene-associated annotation lights up, none of it specific. The difference between LOLA/GAT/regioneR on the *same correct universe* is second-order; the difference between a correct universe and a whole-genome universe on the *same tool* is often p~1 vs p~1e-200. Care belongs on the background, not on tool selection.

2. **A correct null preserves three things, or it manufactures significance:** (a) the regions' **size distribution** (a 50 kb domain hits anything; a 200 bp peak rarely does - relocate intervals of the observed sizes, do not sprinkle points); (b) an **accessible workspace** excluding assembly gaps, centromeres, the ENCODE blacklist (Amemiya 2019), and unmappable bins; (c) **local structure** - GC/isochore, gene density, and clustering (GAT `--isochores`; regioneR `circularRandomizeRegions` for autocorrelated regions). The "right" answer is usually *less* significant than the naive test - that deflation is the methodology working.

3. **One ten-minute sanity move catches most false claims:** shuffle the query within the same workspace and re-run the same overlap pipeline. If shuffled regions also overlap the annotation a lot, the "enrichment" is workspace geography, not biology.

## Method Taxonomy

| Tool | Citation | Null model | When |
|------|----------|-----------|------|
| bedtools fisher | Quinlan 2010 *Bioinformatics* | analytic 2x2; estimates the unobserved "in-neither" cell from mean interval size + genome size; ignores genome structure | fast triage screen only - never the reported result |
| bedtools shuffle + jaccard/numOverlaps | Quinlan 2010 *Bioinformatics* | DIY size-preserving permutation; `-incl`/`-excl` make it matched | entry-level permutation in a shell/Python pipeline; full control, more code |
| GAT | Heger 2013 *Bioinformatics* | per-isochore size-preserving simulation in a workspace; GC/composition conditioned; built-in FDR across annotations | composition-aware enrichment vs many tracks with multiple-testing control |
| regioneR | Gel 2016 *Bioinformatics* | flexible permutation; `randomizeRegions` (mask) vs `circularRandomizeRegions` (preserves clustering) vs `resampleRegions` (real universe); any evaluator; `localZScore` | publication-grade, R-native, autocorrelated regions, where-in-the-region probing |
| LOLA | Sheffield & Bock 2016 *Bioinformatics* | NOT a shuffle - Fisher's exact of query vs a region database, relative to a `userUniverse` | ranking enrichment of a peak set against ENCODE/Roadmap region collections |
| GREAT / rGREAT | McLean 2010 *Nat Biotechnol*; Gu & Hubschmann 2023 *Bioinformatics* | gene regulatory-domain (basal 5 kb up/1 kb down, extend to 1 Mb); binomial-over-regions AND hypergeometric-over-genes | GO/ontology enrichment *from regions* (cis-regulatory), not from a gene list |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Quick "is this even worth permuting?" | `bedtools fisher` | one command, analytic; p~1 means stop; low means permute |
| Publication-grade colocalization of two region sets | regioneR `permTest` with a `mask` | flexible, reports p + z-score; circular randomization for clustered regions |
| Need GC/isochore/composition control + FDR over many tracks | GAT with `--workspace` + `--isochores` | per-isochore sampling removes the GC confounder; built-in qvalue |
| Enrich a peak set against a region database (TFBS/chromatin) | LOLA `runLOLA` with a correct `userUniverse` | universe-relative Fisher; ranked odds ratios + q-values |
| GO/ontology terms FROM regions (cis-regulatory) | rGREAT with a real background | regulatory-domain model; require binomial AND hypergeometric |
| GWAS/eQTL statistical colocalization (shared causal variant) | -> causal-genomics/colocalization-analysis | DISTINCT problem: coloc/SuSiE on summary stats, not interval overlap |
| Gene-list (not region) ontology enrichment | -> pathway-analysis/go-enrichment | start from genes; GREAT is the region-based analog |
| Query regions are autocorrelated/clustered | regioneR `circularRandomizeRegions` | rotating the set preserves inter-region spacing; uniform randomization inflates significance |
| CNV/SV concordance between call sets | `bedtools intersect -f 0.5 -r` (50% reciprocal) | "same event" convention; one-sided fractions let a giant call swallow a tiny one |
| Peaks not yet called | -> chip-seq/peak-calling, atac-seq/atac-peak-calling | this skill operates on existing interval sets |

## bedtools fisher - The Analytic Screen (weakest null)

```bash
bedtools fisher -a peaks.bed -b enhancers.bed -g genome.txt   # both inputs sorted; genome file required
```

`fisher` builds a 2x2 table (in-A/not x in-B/not) and runs Fisher's exact test. The trap: it cannot observe the "in-neither" cell (there is no negative class of intervals that do not exist), so it **estimates the table totals from a heuristic on mean interval size and genome size**, assuming intervals are independent points uniformly