Bedtools Tutorial: Genomic Interval Arithmetic You Will Use Every Week

Quick answer: Bedtools is the Swiss-army knife for genomic intervals. You hand it two or more files of genomic coordinates — BED, GFF, VCF, even BAM — and it answers set-style questions: which reads fall inside these peaks, what is the nearest gene to this variant, how much of the genome do these exons cover. The whole toolchain rests on one idea: intervals are rows in a text file, and once you can sort, intersect, merge, and subtract them, most daily bioinformatics data-wrangling problems collapse into a one-line command.

Companion to the file formats guide and the genome annotation article on this blog — bedtools is where those formats start talking to each other.

The one rule that prevents 90% of mistakes: coordinates

Every interval format answers “where?” slightly differently, and mixing them silently shifts your results:

• BED is 0-based, half-open: start is the first base (counting from 0), end is one past the last base. A feature covering the first 10 bases of a chromosome is start 0, end 10. Length = end - start.
• GFF/GTF is 1-based, closed: the same feature is start 1, end 10. Length = end - start + 1.
• VCF positions are 1-based; SAM/BAM positions are 1-based (despite the binary format).

Bedtools assumes the coordinates you feed it are already in the format the file claims. It does not check. The classic symptom is every intersect coming back off by one base — or off by one direction (one-base-short, one-base-long), which looks like a rounding problem but is a base-system problem. Convert deliberately:

# GFF3 (1-based) -> BED (0-based): subtract 1 from the start column
awk 'BEGIN{OFS="\t"} !/^#/ {print $1, $4-1, $5, $9, 0, $7}' genes.gff3 > genes.bed

Two more prerequisites before any command runs:

• Sort order must match. Either sort both files with bedtools sort -i in.bed > sorted.bed, or — for large files — sort with sort -k1,1 -k2,2n and pass the -sorted flag so bedtools streams instead of loading everything into memory.
• The genome file must be real. Several tools (genomecov, slop, complement) need a genome file: two tab-separated columns, chromosome and length (samtools faidx ref.fa produces one as .fai). Using the wrong assembly’s lengths gives plausible-looking nonsense.

The core verbs

Install (conda install -c bioconda bedtools) and the daily drivers, in order of how often you will type them:

# 1. intersect: which intervals from A overlap B?
bedtools intersect -a reads.bed -b peaks.bed                 # only overlapping A rows
bedtools intersect -a reads.bed -b peaks.bed -wa -wb         # report both sides
bedtools intersect -a reads.bed -b peaks.bed -wo             # ...plus overlap length
bedtools intersect -a reads.bed -b exons.bed -v              # A rows with NO overlap

# 2. subtract / closest
bedtools subtract -a genes.bed -b repeats.bed               # genes minus repeat regions
bedtools closest -a variants.bed -b genes.bed -d            # nearest gene + distance

# 3. merge / cluster
bedtools merge -i sorted.bed                                # fuse overlapping/abutting intervals
bedtools merge -i sorted.bed -c 4 -o count                  # ...counting original intervals

# 4. coverage and depth
bedtools genomecov -i reads.bed -g genome.fa.fai -bg        # per-base depth, BEDGRAPH
bedtools coverage -a regions.bed -b reads.bed               # count + fraction covered

# 5. convenience
bedtools flank -i genes.bed -g genome.fa.fai -b 1000        # 1 kb windows either side
bedtools slop -i peaks.bed -g genome.fa.fai -b 200          # extend peaks by 200 bp
bedtools getfasta -fi genome.fa -bed orfs.bed               # extract the actual sequence
bedtools jaccard -a a.bed -b b.bed                           # overlap statistic between sets

A worked question most analysts meet early — how many of my ChIP-seq peaks overlap annotated exons? — is exactly two commands: bedtools sort both files, then bedtools intersect -a peaks.bed -b exons.bed -u | wc -l.

Where it quietly goes wrong

• Strand is ignored unless you ask. intersect does not care about +/- by default; add -s (same strand) or -S (opposite strand) when direction matters, e.g. for paired reads or antisense overlap.
• Double counting. -wa reports each A row once per overlapping B row, so a read hitting three exons appears three times. Use -u (once) or -c (count) when you want one row per input.
• Unsorted input with -sorted. The flag trades memory for correctness only if both files are genuinely in the same chromosome order; otherwise results are silently truncated.
• BED12/block tricks. Spliced alignments and isoforms live in BED12; plain three-column BED treats an entire intron-containing span as one interval, which inflates overlap with intronic features.
• Huge files. Everything in BED is text. Compress (bgzip) and index (tabix) interval files when they exceed a few hundred megabytes, and prefer -sorted streaming.

What a defensible interval analysis hands over

The deliverable is not just the resulting BED: it is the BED plus the genome build it was computed against, the exact bedtools version, and the sorting/filtering steps that produced the inputs. Interval arithmetic is cheap and fast, which makes it dangerously easy to re-run in five years against a different assembly and get confident, incompatible answers.

Key takeaways

• BED is 0-based half-open, GFF is 1-based closed — mixing them is the classic off-by-one, and bedtools will not warn you.
• Sort your inputs (and match their chromosome order) before any set operation; use -sorted for large files with a valid genome file.
• Learn five verbs cold: intersect, subtract, closest, merge, coverage — they cover most daily questions.
• Use -u/-c instead of plain -wa when you need one row per input interval rather than one row per overlap.
• Strand and multi-overlap semantics are opt-in: decide them explicitly for every analysis.
• Version, genome build, and input provenance belong in the methods text of any report built on interval arithmetic.

Challenge

Take a BED file of promoters (flank the TSS of every gene with bedtools flank -b 500) and a BED of variants from a public VCF (grep -v '^#' in.vcf | awk 'BEGIN{OFS="\t"}{print $1, $2-1, $2}'). Compute: how many variants fall in promoters, how many of those are in repeats (subtract), and what the nearest gene is for the rest (closest -d). Every one of those answers is a single command — if you needed more than one, the interval probably was not in the form you assumed.

Want guided practice turning coordinate files into defensible results? Ampersand Academy offers one-to-one mentoring in bioinformatics — part of Mahadhi, a development and digital marketing company.

FAQ