/bedtools-genomic-intervals

name: "bedtools-genomic-intervals"
description: "Genomic interval ops on BED/BAM/GFF/VCF. Find overlaps, merge intervals, compute coverage, extract FASTA, find nearest features. Core for ChIP-seq peak annotation, region filtering, genome arithmetic. Use tabix for indexed single-region queries; use deeptools for normalized bigWig coverage."
license: "GPL-2.0"

bedtools — Genomic Interval Analysis Toolkit

Overview

bedtools is the standard toolkit for operating on genomic intervals in BED, BAM, GFF, and VCF formats. It solves the core problem of genome arithmetic: finding overlaps between feature sets, computing coverage, extracting sequences, merging adjacent regions, and annotating features with nearest neighbors. bedtools operates on sorted coordinate lists and runs at C speed, making it practical for whole-genome analyses.

When to Use

• Intersecting ChIP-seq peaks with gene annotations to find promoter-overlapping peaks
• Merging overlapping ATAC-seq peaks or called regions across replicates
• Computing read coverage depth over target capture regions
• Extracting FASTA sequences for motif discovery or primer design
• Finding the nearest gene to each regulatory element or variant
• Subtracting blacklist or repeat regions from peak calls
• Expanding genomic intervals by fixed distance (promoter regions)
• Use `tabix` instead for fast indexed queries of a single genomic region
• For normalized coverage bigWig tracks, use `deeptools bamCoverage` instead
• Use `mosdepth` instead for whole-genome per-base depth (10× faster)

Prerequisites

• **Python packages**: None required (command-line only)
• **Input requirements**: BED/BAM/GFF/VCF files; FASTA reference for `getfasta`; genome file (chromosome sizes) for `slop`/`flank`/`genomecov`
• **Sorting**: Most operations require coordinate-sorted input

> **Check before installing**: The tool may already be available in the current environment (e.g., inside a `pixi` / `conda` env). Run `command -v bedtools` first and skip the install commands below if it returns a path. When running inside a pixi project, invoke the tool via `pixi run bedtools` rather than bare `bedtools`.

# Bioconda (recommended)
conda install -c bioconda bedtools

# Homebrew (macOS)
brew install bedtools

# Verify
bedtools --version
# bedtools v2.31.0

# Create genome file from FASTA index
samtools faidx reference.fa
cut -f1,2 reference.fa.fai > genome.txt  # chr → size table

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: overlapRequirement
    kind: required
    source: user
    ask: "How much of a feature must overlap before the two count as intersecting - a single base, or a stated fraction?"
    default: "one base"

  - id: D2
    param: reciprocalOverlap
    kind: required
    source: user
    depends_on: [D1]
    ask: "Must the fraction hold for both features, or only for the query?"
    default: "query only"
    skip_if: "overlap requirement left at a single base"

  - id: D3
    param: strandedness
    kind: required
    source: data
    ask: "Do the features carry strand, and must two features share it to be considered overlapping?"
    default: "strand ignored"

  - id: D4
    param: mergeDistance
    kind: optional
    source: user
    ask: "How close may two intervals be before they are merged into one?"
    default: "touching or overlapping only"

  - id: D5
    param: coverageScaling
    kind: optional_conditional
    source: upstream
    ask: "Should coverage be scaled to a per-million factor so samples are comparable?"
    default: "raw coverage"

  - id: D6
    param: sortedInput
    kind: never_ask
    source: data
    reason: "The sweep algorithm needs sorted input and changes memory and speed, not the intervals returned"
    default: "used when inputs are sorted"

D3 is the quiet one. For stranded features - genes, stranded RNA coverage, motif hits - ignoring strand reports antisense overlaps as real, and the output is a perfectly ordinary BED file with too many rows.

Quick Start

# Find peaks overlapping genes, then merge overlapping peaks
bedtools intersect -a peaks.bed -b genes.bed -wa -wb > peaks_with_genes.bed
bedtools merge -i peaks.bed > merged_peaks.bed
bedtools coverage -a genes.bed -b reads.bam > gene_coverage.bed

Core API

Module 1: Interval Intersection and Overlap Analysis

Find regions that overlap between two feature sets.

# Basic intersection: output overlapping regions
bedtools intersect -a peaks.bed -b genes.bed

# Report original A and B features for each overlap
bedtools intersect -a peaks.bed -b genes.bed -wa -wb

# Count B overlaps per A feature (adds column)
bedtools intersect -a peaks.bed -b genes.bed -c
# Output: chr1  1000  2000  peak1  gene_count

# Peaks with ANY overlap (report each peak once)
bedtools intersect -a peaks.bed -b genes.bed -u

# Peaks with NO overlap in B (invert filter)
bedtools intersect -a peaks.bed -b blacklist.bed -v

# Require reciprocal 50% overlap both ways
bedtools intersect -a exp1.bed -b exp2.bed -f 0.5 -F 0.5 -r

# Same-strand intersections only
bedtools intersect -a peaks.bed -b genes.bed -s

# Multiple database files with overlap counts per file
bedtools intersect -a query.bed -b enhancers.bed promoters.bed \
    -names enh prom -C

# Memory-efficient mode for pre-sorted large files
bedtools intersect -a sorted_peaks.bed -b sorted_genes.bed -sorted

Module 2: Interval Merging and Arithmetic

Combine overlapping intervals and perform set operations.

# Merge overlapping and adjacent intervals
sort -k1,1 -k2,2n peaks.bed | bedtools merge -i stdin

# Merge intervals within 500 bp of each other
bedtools merge -i peaks.bed -d 500

# Merge and count original features
bedtools merge -i peaks.bed -c 1 -o count
# Output: chr1  1000  5000  3 (3 original peaks merged)

# Merge and collapse feature names
bedtools merge -i peaks.bed -c