Before you start
- You are comfortable with the common interval formats from File formats, and have basic command-line comfort (see Bash basics).
- Any term new? The Glossary has it.
Learning objectives
By the end of this lesson you will be able to: read a genomic interval as a row of (chrom, start, end), describe what the core bedtools subcommands do (intersect, merge, subtract, closest, complement, genomecov, slop and flank), and write an intersect and a sort-then-merge command correctly.
Intervals are just rows of numbers
A genomic interval is a region on a chromosome. In a BED file it is written as three columns: the chromosome, the start, and the end. Because BED is 0-based and half-open, the start is counted from 0 and the end is one past the last base, so the length of an interval is simply end - start with no "+1" to remember.
# peaks.bed (chrom start end) -- 0-based, half-open
chr1 100 200
chr1 400 450
chr2 50 90
Once your features are rows like this, bedtools treats them like math. "Where do two sets of intervals overlap?" is set intersection. "Combine touching intervals into one" is a kind of addition. "Remove these regions from those" is subtraction. The whole toolkit is genome arithmetic.
1The core subcommands
You will reach for a handful of subcommands again and again. Here is what each one answers.
| Subcommand | What it does | Handy flags |
|---|---|---|
| intersect | Find where intervals in A overlap intervals in B | -wa, -wb, -v |
| merge | Combine overlapping or touching intervals into one (needs sorted input) | -d |
| subtract | Remove the parts of A that overlap B | -A |
| closest | For each interval in A, find the nearest interval in B | -d |
| complement | Return the regions of the genome NOT covered by A | needs a genome file |
| genomecov | Report coverage depth across the genome | -bg, -d |
| slop / flank | Extend intervals (slop) or make new flanking regions (flank) | -b, -l, -r |
Why -wa and -wb matter
By default intersect reports only the overlapping piece. -wa ("write A") gives you the full original A interval, -wb ("write B") appends the full B interval, and -v inverts the question to report A intervals that do NOT overlap B at all.
2Picturing intersect
Intersect lines up two tracks and reports where they overlap. In the figure below, track A (peaks) and track B (genes) share one overlapping region; that shared stretch is what a plain intersect returns.
3Worked example: intersect with -wa -wb
Say you have ChIP-seq peaks and a set of gene intervals, and you want to know which peaks fall inside which genes, keeping the full record of both. That is one command:
# report each peak (A) alongside the gene (B) it overlaps bedtools intersect -a peaks.bed -b genes.bed -wa -wb # only the peaks that overlap NO gene at all bedtools intersect -a peaks.bed -b genes.bed -v
The first command prints the full peak followed by the full gene on each line where they overlap. The second flips the question and lists the peaks with no gene overlap.
4Worked example: sort then merge
Many subcommands, including merge, require the input to be sorted by chromosome and then by start position. The standard sort is sort -k1,1 -k2,2n: sort by column 1 as text, then by column 2 numerically.
# sort first (chrom as text, start numerically), then merge overlaps sort -k1,1 -k2,2n peaks.bed > peaks.sorted.bed bedtools merge -i peaks.sorted.bed # or pipe it in one go sort -k1,1 -k2,2n peaks.bed | bedtools merge -i -
Why sorting matters
merge walks through intervals in order and fuses each one with the next if they overlap or touch. If the file is not sorted, neighbouring intervals are not actually next to each other in the file, so the merge is wrong. Sort first, every time.
Check your understanding
bedtools intersect report A intervals that do NOT overlap B?bedtools merge, what must you usually do to the input?