BEDTools Intersect Tutorial: Step-by-Step Guide

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-10

This BEDTools intersect tutorial takes you from a clean terminal to a verified overlap analysis you can defend in a paper or a lab meeting. By the end you will know how to run bedtools intersect on BED, GFF, VCF and BAM inputs, how to read its default output (which surprises almost everyone the first time), how to control overlap thresholds, how to make it fast on multi-gigabyte files, and how to prove the answer is correct before you build anything on top of it.

You need four things on hand. First, BEDTools installed and reporting a version of at least 2.30 on bedtools, version. Second, two interval files, one used as -a and one as -b. Third, a genome file (chromosome names and sizes) if you plan to use the sorted-input mode. Fourth, samtools if any of your inputs are BAM. Everything else in this guide is optional.

BEDTools is open-source software maintained by the Quinlan laboratory and distributed through Bioconda, Homebrew and direct compilation. The canonical reference is Quinlan AR and Hall IM, "BEDTools: a flexible suite of utilities for comparing genomic features," Bioinformatics 26(6):841-842, 2010. The command-level behavior described here follows the official BEDTools documentation at bedtools.readthedocs.io.

What bedtools intersect Actually Does

bedtools intersect compares two sets of genomic intervals and reports where they overlap. It reads intervals from the -a file and tests each one against the intervals in the -b file. Anything in -a that shares at least one base with anything in -b is an overlap.

The term "interval" here means a chromosome name plus a start and end position. That is the minimum requirement. BEDTools also understands strand, name, score and block structure when those columns are present.

The tool answers four distinct questions depending on the flags you pass:

1. Which regions of A overlap B, and where exactly inside A does the overlap sit?
2. Which full A features overlap B (keeping A's original coordinates)?
3. Which A features do not overlap B at all?
4. How many bases of overlap exist for each A-B pair?

Most confusion in this tutorial traces back to running question 1 when the user wanted question 2. That distinction drives the first half of the walkthrough below.

The Bioinformatics File Formats You Need First

bedtools intersect consumes interval files, not sequence files. Knowing which file is which prevents the single most common wasted hour in this workflow.

BED

BED is a tab-delimited text format with no header. The first three columns are mandatory: chromosome, start, end. Columns four through twelve are optional and include name, score, strand, thickStart, thickEnd, itemRgb, blockCount, blockSizes and blockStarts.

BED coordinates are zero-based for the start and half-open for the end. A record chr1 100 200 covers the 101st through the 200th base in one-based counting and spans exactly 100 bases. That arithmetic matters when you compare BEDTools output to a genome browser.

GFF and GTF

GFF3 and GTF are one-based, inclusive, tab-delimited annotation formats. They carry feature type, source, strand, phase and attributes. BEDTools reads them natively and detects the format from the file extension and the column structure. Because GFF is one-based and BED is zero-based, the same physical feature looks one base different in the two files. BEDTools normalizes internally, so an overlap is still an overlap. The output is written in the format of the -a file.

VCF

Variant Call Format is one-based and describes positions plus reference and alternate alleles. BEDTools treats each VCF record as a one-base interval at the POS field unless you convert it to BED first with a wider window. Use -header to preserve the VCF header lines in the output.

BAM

BAM is the binary, compressed version of SAM alignment format. It stores reads with their aligned coordinates, CIGAR strings, mapping quality and flags. BEDTools reads BAM directly and can write BAM back out. This is the format you intersect when you care about reads rather than annotations.

FASTA versus FASTQ

FASTA and FASTQ are sequence formats, and neither has genomic coordinates. That single fact excludes both from bedtools intersect input.

FASTA stores a sequence under a header line beginning with >. It carries nucleotides or amino acids and nothing else. Reference genomes, transcript sequences and protein databases are distributed as FASTA.

FASTQ stores a sequence under a header line beginning with @, followed by the sequence, a + separator line, and a quality string with one ASCII character per base. The quality characters encode Phred scores. Raw sequencing reads arrive as FASTQ because the instrument produces both bases and per-base confidence.

The practical difference for this tutorial: FASTQ is what comes off the sequencer, FASTA is what you align against, and BAM is what the aligner produces. Only the BAM (or a BED derived from it) is an interval file. If someone hands you a FASTQ and asks you to intersect it with a BED file, the correct first step is alignment with BWA, Bowtie2, STAR or HISAT2, not bedtools intersect.

Coordinate systems in one sentence

BEDTools converts every input to zero-based, half-open coordinates internally, so mixing BED with GFF or VCF produces correct overlaps even though the raw files disagree by one base.

Choosing A and B

The -a file is the query. The -b file is the reference set. Order matters for output format, for computational cost and for what the columns mean.

A useful rule: put the file whose rows you want to keep in the -a slot. If you want one row per gene, genes go in -a. If you want one row per peak, peaks go in -a.

BEDTools loads the -b file into memory when the sorted mode is off. That makes -b the expensive slot at scale. A 5-million-interval BED file as -b can consume a large amount of RAM. Swap the files and the memory profile changes completely.

Step-by-Step Walkthrough

Step 1: Install BEDTools and confirm the version

mamba install -c bioconda -c conda-forge bedtools
bedtools --version

Expected output:

bedtools v2.31.1

Any 2.30 or later release supports every flag in this tutorial. If you see "command not found" after installation, your conda environment is not activated or the binary directory is not on your PATH.

Step 2: Build a small test dataset

Never debug a pipeline on real data. Build a toy set first.

printf 'chr1\t100\t200\tgeneA\t0\t+\nchr1\t300\t400\tgeneB\t0\t-\nchr1\t500\t600\tgeneC\t0\t+\nchr2\t100\t250\tgeneD\t0\t+\n' > genes.bed

printf 'chr1\t150\t250\tpeak1\nchr1\t550\t580\tpeak2\nchr2\t50\t120\tpeak3\n' > peaks.bed

Note the \t separators. BED files are tab-delimited, and spaces will break the parser later. Confirm the files look right:

cat -A genes.bed | head -2

You should see ^I between fields, which is how cat -A renders a tab.

Step 3: Run the simplest possible intersect

bedtools intersect -a genes.bed -b peaks.bed

Expected output:

chr1	150	200
chr1	550	580
chr2	100	120

Three lines, one per overlapping gene. Now look carefully at the coordinates. geneA runs from 100 to 200, but the output says 150 to 200. geneC runs 500 to 600, but the output says 550 to 580. BEDTools has clipped each gene down to the portion that actually overlaps a peak.

This is the default behavior, and it is the number one source of "my output is wrong" messages on bioinformatics forums. The default reports the intersection, not the original feature.

Step 4: Keep the original A intervals with -wa

bedtools intersect -a genes.bed -b peaks.bed -wa

Expected output:

chr1	100	200	geneA	0	+
chr1	500	600	geneC	0	+
chr2	100	250	geneD	0	+

Now the coordinates are the genes as you wrote them. All six columns survive. Use -wa whenever downstream analysis needs the original feature boundaries, which is almost always.

Step 5: See which B feature caused each hit with -wb

bedtools intersect -a genes.bed -b peaks.bed -wa -wb

Expected output:

chr1	100	200	geneA	0	+	chr1	150	250	peak1
chr1	500	600	geneC	0	+	chr1	550	580	peak2
chr2	100	250	geneD	0	+	chr2	50	120	peak3

The A columns come first, then the B columns. -wb on its own (without -wa) appends B to the clipped A interval, which is rarely what you want.

Step 6: Find genes with no peak at all using -v

bedtools intersect -a genes.bed -b peaks.bed -v

Expected output:

chr1	300	400	geneB	0	-

-v inverts the logic and reports only A features with zero overlap. Note that -v returns the full original A record, not a clipped one, because a clipped interval of length zero would be meaningless.

Step 7: Use -u when you want each A once

If a gene overlaps five peaks, -wa prints it five times. -u prints it once.

bedtools intersect -a genes.bed -b peaks.bed -u

The output is identical to Step 4 here because each gene happens to hit exactly one peak. On real data with clustered peaks, -u is the flag that prevents inflated counts.

Step 8: Quantify overlap with -wo

bedtools intersect -a genes.bed -b peaks.bed -wo

Expected output:

chr1	100	200	geneA	0	+	chr1	150	250	peak1	50
chr1	500	600	geneC	0	+	chr1	550	580	peak2	30
chr2	100	250	geneD	0	+	chr2	50	120	peak3	20

The final column is the number of overlapping bases. Check the arithmetic by hand: geneA is 100 bp, peak1 starts at 150, so the shared region is 150 to 200, which is 50 bases. Correct.

Step 9: Keep non-overlapping A with -wao

-wo drops genes with no overlap. -wao keeps them with a zero in the overlap column and placeholder values where B fields would go.

bedtools intersect -a genes.bed -b peaks.bed -wao

Expected output:

chr1	100	200	geneA	0	+	chr1	150	250	peak1	50
chr1	300	400	geneB	0	-	.	-1	-1	-1	0
chr1	500	600	geneC	0	+	chr1	550	580	peak2	30
chr2	100	250	geneD	0	+	chr2	50	120	peak3	20

The placeholder pattern is a dot for the first B column and -1 for the remaining B columns. -loj gives you the same left-outer-join shape without the overlap count. Row counts after -wao are always greater than or equal to the number of A features, because a single A can produce multiple rows when it hits multiple B features.

Step 10: Apply overlap thresholds with -f, -F, -r and -e

The -f flag sets the minimum fraction of A that must be covered. The default is 1E-9, which means "any single base."

bedtools intersect -a genes.bed -b peaks.bed -wa -f 0.5

Expected output:

chr1	100	200	geneA	0	+

Only geneA passes. Its overlap is 50 of 100 bases, exactly 0.5. geneC covers 30 of 100 bases and fails. geneD covers 20 of 150 bases and fails.

The companion flags:

• -F 0.5 requires half of the B feature to be covered. Use it when B features are small and A features are large, for example asking which transcription factor peaks fall entirely inside a gene body.
• -r requires both -f and -F to be satisfied simultaneously.
• -e requires either -f or -F to be satisfied.

These four flags compose. -f 0.5 -r means "at least half of A and at least half of B are in the overlap."

Step 11: Add strand awareness with -s and -S

Both files need a strand column for this to work.

bedtools intersect -a genes.bed -b peaks.bed -wa -s

-s keeps only hits where A and B share the same strand. -S keeps only opposite-strand hits. Without a strand column in either file, BEDTools stops with a strandedness error, covered below.

Step 12: Intersect one A against many B files

bedtools intersect -a genes.bed -b peaks_rep1.bed peaks_rep2.bed -wa -wb -names rep1 rep2

-names adds a label column identifying which B file produced each hit. Without it, you cannot tell the replicates apart. -filenames does the same using file paths instead of labels you supply.

Step 13: Switch on the sorted algorithm for large files

sort -k1,1 -k2,2n genes.bed > genes.sorted.bed
sort -k1,1 -k2,2n peaks.bed > peaks.sorted.bed

bedtools intersect -a genes.sorted.bed -b peaks.sorted.bed -wa -sorted

The -sorted flag replaces the in-memory index with a line-sweep algorithm that streams both files. Peak memory drops dramatically and runtime improves on large inputs. Both files must be sorted by chromosome and then by start position, and the chromosome order must match between the two files. Sorting both with the same command guarantees that.

If your chromosome order is not plain lexicographic, supply a genome file with -g:

bedtools intersect -a genes.sorted.bed -b peaks.sorted.bed -wa -sorted -g hg38.genome

A genome file is a two-column tab-delimited file of chromosome name and length. The .fai index produced by samtools faidx works directly as a genome file.

Step 14: Work with BAM input

bedtools intersect -a reads.bam -b targets.bed -bed > hits.bed

When -a is BAM, the default output is BAM. -bed writes BED12 instead, which is easier to inspect. Use -ubam if you need uncompressed BAM for a downstream tool that cannot read BGZF.

For RNA-seq alignments with spliced reads, add -split. This tells BEDTools to treat the N CIGAR operations as separate blocks rather than one long interval spanning the intron.

bedtools intersect -a rnaseq.bam -b exons.bed -split -bed > exon_hits.bed

Without -split, a read spanning a 10 kb intron is treated as covering the whole 10 kb. That inflates overlap counts badly for spliced data.

Options at a glance

Common Errors and Their Fixes

These are the messages you will actually see. The wording may vary slightly by version.

"Error: Unable to open file"

Error: Unable to open file peaks.bed. Exiting.

The path is wrong, the file is a directory, or the extension is not recognized. Check with ls -l peaks.bed. BEDTools infers the file type partly from the extension, so a BED file named peaks.txt may be rejected or misparsed. Rename it to .bed.

"Error: Sorted input specified, but the file has the following out of order record"

Error: Sorted input specified, but the file peaks.bed has the following out of order record
chr1	500	600	...

You used -sorted on an unsorted file. The message prints the offending record so you can find the break point. Fix it with bedtools sort -i peaks.bed > peaks.sorted.bed or with sort -k1,1 -k2,2n.

"chromosome ... found in your BED file but not in your genome file"

Your genome file is incomplete or uses different chromosome names from your BED file. A common cause is a genome file listing 1, 2, 3 while the BED file uses chr1, chr2, chr3. Regenerate the genome file from the same assembly your BED file was built on.

"Error: Type checker found wrong number of fields"

A line has a different column count from the rest of the file. This usually comes from a stray tab, a trailing tab, or a header line that was not removed. Find it with:

awk -F'\t' '{print NF}' genes.bed | sort -n | uniq -c

Every line should report the same field count.

"Error: malformed BED entry at line N"

The start coordinate is greater than the end coordinate, or a coordinate is not an integer. BEDTools cannot represent an inverted interval. Check with:

awk -F'\t' '$3 <= $2 {print NR": "$0}' genes.bed

Empty output with no error

The command runs cleanly and prints nothing. Three causes dominate. Chromosome naming differs between the files (chr1 versus 1, or chr1 versus Chr1). The coordinate systems are mixed without conversion, shifting everything by one base. Or the two files genuinely do not overlap, which is common when one is on hg19 and the other on hg38.

Check the chromosome vocabularies directly:

cut -f1 genes.bed | sort -u
cut -f1 peaks.bed | sort -u

If those lists differ, fix the naming before anything else.

Strandedness errors

Passing -s or -S to a three-column BED file produces an error stating that strand information is unavailable. Add a strand column, or drop the flag.

-split rejected on BED input

-split only applies to BAM alignments with CIGAR strings. On a BED file it has no meaning and BEDTools will reject it or silently ignore it depending on version. Convert to BAM first if you truly need block-level splitting.

How to Verify the Result Is Correct

Run three checks before trusting any overlap count.

Check 1: The partition identity. Every A feature either overlaps B or does not. Those counts must add up.

total=$(wc -l < genes.bed)
nohit=$(bedtools intersect -a genes.bed -b peaks.bed -v | wc -l)
hit=$(bedtools intersect -a genes.bed -b peaks.bed -u | wc -l)
echo "$((nohit + hit)) vs $total"

The two numbers must match. If they do not, you have a duplicate-row problem, most often from overlapping B intervals.

Check 2: Hand-check one overlap. Pick a random output row with -wo and recompute the overlap length manually as min(A_end, B_end) - max(A_start, B_start). If the reported number matches, the parser read your file correctly.

Check 3: Confirm the coordinate convention. Run the intersect with -a and -b swapped. The set of overlapping pairs should be identical even though the output looks different. If the pair sets differ, one file is being misparsed.

For VCF and GFF inputs, add a fourth check: convert both files to BED with bedtools slop or a converter, run the intersect again, and confirm the same pairs come back. That catches off-by-one errors introduced before BEDTools ever sees the data.

Scaling Up: Performance Notes

The default mode builds an interval index of the -b file in memory. Memory scales with the number of B intervals and their span. For a few hundred thousand intervals this is trivial. For tens of millions, it can exhaust a workstation.

-sorted changes the algorithm to a sweep that walks both sorted files in lockstep. Memory use becomes roughly constant, and runtime scales close to linearly. The trade-off is the sorting requirement and the requirement that chromosome ordering match.

For very large jobs, bedtools intersect itself is single-threaded. Parallelize by splitting the A file by chromosome and running one process per chromosome, then concatenating. That keeps memory bounded per process and uses all cores.

for c in $(cut -f1 genes.bed | sort -u); do
  grep -P "^${c}\t" genes.sorted.bed > work/${c}.a.bed
  grep -P "^${c}\t" peaks.sorted.bed > work/${c}.b.bed
  bedtools intersect -a work/${c}.a.bed -b work/${c}.b.bed -wa -sorted > work/${c}.hits.bed &
done
wait
cat work/*.hits.bed > all_hits.bed

GNU parallel or a workflow manager such as Snakemake or Nextflow does the same thing with better bookkeeping. The -sorted flag stays valid inside each per-chromosome job because a single chromosome is trivially sorted.

Common Mistakes and Limitations

Mistaking the default output for the original A interval. This is the mistake that produces most silent errors. The default output is the clipped overlap. Add -wa unless you specifically want the intersection region.

Using spaces instead of tabs. BED is tab-delimited. Files edited in a spreadsheet and saved as "text" often contain spaces. BEDTools will either merge columns or throw a field-count error. Run cat -A to check.

Forgetting that overlapping B intervals multiply rows. If peak1 and peak2 overlap each other, an A feature hitting both region produces two rows with -wa. Counts inflate. Collapse B first with bedtools merge -i peaks.bed > peaks.merged.bed when you want a union.

Mixing genome builds. A BED file on hg19 intersected against a BED file on hg38 gives near-zero output because the coordinate systems refer to different assemblies. Liftover with liftOver or the UCSC tools before intersecting. This is the most common silent failure in published analyses.

Assuming chromosome naming is universal. Ensembl uses 1, UCSC uses chr1, and some pipelines strip the prefix. Normalize with sed 's/^chr//' or the reverse, applied consistently to both files.

Treating one-based and zero-based formats as interchangeable. BEDTools handles the conversion internally, but the moment you export to a tool that does not, a one-base shift appears. Verify the convention every time you cross a tool boundary.

Using -f 1.0 and expecting "fully contained." A feature is fully contained only if the overlap equals the full A length under the same coordinate convention. Because -f is computed on floats, boundary cases can behave unexpectedly. For strict containment, -F 1.0 with B as the smaller feature is often the cleaner formulation.

Running -sorted without matching sort orders. Two files sorted with different tools can produce different chromosome orders. Sort both with the same command, or use -g.

Applying bedtools intersect to FASTQ or FASTA files. These formats have no coordinates. Align the reads first, then intersect the resulting BAM.

Confusing intersect with coverage. bedtools intersect answers "do these overlap." bedtools coverage answers "what depth or breadth of coverage does B provide across A." If you need per-base depth, coverage is the correct tool.

Assuming the tool validates your biology. bedtools intersect checks syntax and coordinate consistency. It does not check that your BED file came from the assembly you think it did, that your peaks are from the right实验, or that your gene annotations match your variant calls. Those checks are yours.

Individual research questions and clinical interpretations need input from a qualified domain expert before you act on them.

Frequently Asked Questions

What is the difference between bedtools intersect and bedtools coverage?

bedtools intersect reports whether intervals overlap and where. bedtools coverage reports how much of each A feature is covered by B, including depth at each position. Use intersect for presence and coverage for depth or breadth.

Does bedtools intersect report the full A interval or just the overlap?

By default it reports only the portion of A that overlaps B. Add -wa to get the original A interval back.

Can bedtools intersect read VCF files directly?

Yes. BEDTools parses VCF natively and treats each record as a one-base interval at POS unless you widen it. Use -header to keep the VCF header in the output.

Why is my bedtools intersect output empty?

The most common causes are mismatched chromosome naming, mismatched genome builds, and genuinely non-overlapping regions. Compare the chromosome name lists from both files first.

What does -f 1.0 mean?

It requires 100 percent of the A interval to be covered by B. Any uncovered base fails the test. For strict containment checks, -F 1.0 with the smaller feature as B is often more intuitive.

Do I need to sort my files before running bedtools intersect?

No for the default mode, which builds an in-memory index. Yes for -sorted, which requires both files sorted by chromosome and start position with matching chromosome order.

Can bedtools intersect work on FASTQ files?

No. FASTQ and FASTA are sequence formats without coordinates. Align the reads to produce BAM, then intersect the BAM or a BED derived from it.

How do I get only the A features that do not overlap B?

Use -v. It returns the full original A record for every feature with zero overlap, which is the complement of the -u result.