Bio Genome Intervals Proximity Operations

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

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

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.

# 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

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.

# 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/right); -sw define -l/-r BY STRAND; -sm/-Sm keep only same/opposite-strand B; -u boolean (A once if any B); -c count of B per A; -v A with no B in window. -sw controls where the window is; -sm/-Sm control which B count - distinct concerns.

slop / flank - Extend or Find Adjacent Regions (genome file REQUIRED)

Both need -g genome.txt precisely so they can clip at chromosome boundaries - extension past coordinate 0 or past chrom length is silently truncated (start floored at 0, end capped). Flags: -b N both sides; -l N/-r N per side (coordinate unless -s); -s strand-aware (on a --strand feature -l adds to the END, so -l always means "upstream of the feature"); -pct treat N as a fraction of feature length; -header echo input header.

Build Strand-Aware Promoters from a Gene Model

Goal: Produce a promoter BED (TSS -2000 / +200 bp, strand-aware) that is correct for both strands - the right way to define "promoter", which is a choice imposed on a TSS, not an annotated feature.

Approach: Collapse genes to their TSS first (start for +, end-1 for -), THEN slop -s so "upstream" tracks strand. Running slop -b 2000 on a gene BODY is the wrong promoter (it grows the whole gene, ignores strand).

# 1) TSS BED from a BED6 gene model (strand-aware single base)
awk -v OFS='\t' '{ if ($6=="+") print $1,$2,$2+1,$4,$5,$6; else print $1,$3-1,$3,$4,$5,$6 }' genes.bed > tss.bed

# 2) Promoter = TSS -2000 / +200, strand-aware (-l is always the upstream side under -s)
bedtools slop -i tss.bed -g genome.txt -s -l 2000 -r 200 > promoters.bed

import pybedtools

UP = 2000   # bp upstream of TSS; common core-promoter convention, NOT a fact -- report it and tune per assay
DOWN = 200  # bp downstream of TSS; asymmetric on purpose (+1 nucleosome / 5'UTR sit downstream)

genes = pybedtools.BedTool('genes.bed')
tss = genes.each(lambda f: pybedtools.create_interval_from_list([f[0], str(f.start) if f.strand == '+' else str(f.end - 1), str(f.start + 1) if f.strand == '+' else str(f.end), f.name, f.score, f.strand])).saveas()
promoters = tss.slop(g='genome.txt', s=True, l=UP, r=DOWN).saveas('promoters.bed')

flank shares the flag vocabulary but emits the regions BESIDE each feature and drops the original (two intervals per input, used for splice/boundary context):

bedtools flank -i exons.bed -g genome.txt -s -b 1000 > exon_flanks.bed   # 1 kb each side, strand-aware

Per-Method Failure Modes

-D ref silently mis-signs minus-strand genes

Trigger: using closest -D ref and interpreting the sign as upstream/downstream. Mechanism: -D ref signs by genomic coordinate only (lower = negative); for a --strand gene the TSS is at the HIGHER coordinate, so "upstream" runs to higher coordinates and the coordinate sign is inverted relative to biology. Symptom: half the genes (the --strand ones) are folded the wrong way; symmetric QC (TSS-enrichment plot) still looks fine, but any "enhancers preferentially upstream" claim washes out or inverts. Fix: use -D b (sign by the gene's strand) for any upstream/downstream biology; reserve -D ref for pure left/right genomic distance.