intersectBed

Updated 2014 June 25th

The tool intersectBed is part of the BEDTools suite of tools and performs an intersection between two BED files. For example, given two BED files, you may be interested in finding the entries that overlap. To install the latest version of BEDTools, download the source code from GitHub and compile:

git clone https://github.com/arq5x/bedtools2.git
cd bedtools2
make clean; make all
bin/bedtools --version
#bedtools v2.20.1-4-gb877b35

BEDTools usage examples are available in the official documentation.

Getting started

When I first used BEDTools in 2010, intersectBed was a separate program that was part of the BEDTools suite. Since then, BEDTools is released as one main program called bedtools. If we look inside intersectBed, it is actually just a shell script calling bedtools:

cat intersectBed 
#!/bin/sh
${0%/*}/bedtools intersect "$@"

The $@ symbol stores all the input arguments supplied, for example:

cat blah.sh 
#!/bin/sh
echo "$@"

blah.sh -a a.bed -b b.bed
-a a.bed -b b.bed

The ${0%/*} refers to $0 using parameter expansion

cat blah.sh 
#!/bin/sh
echo ${0} "$@"

blah.sh -a a.bed -b b.bed
./blah.sh -a a.bed -b b.bed

and the %/* removes everything before the /

cat blah.sh 
#!/bin/sh
echo ${0%/*} "$@"

blah.sh -a a.bed -b b.bed
. -a a.bed -b b.bed

So in the end:

cat blah.sh 
#!/bin/sh
echo ${0%/*}/bedtools intersect "$@"

blah.sh -a a.bed -b b.bed
./bedtools intersect -a a.bed -b b.bed

I'm used to calling intersectBed, it's exactly the same as calling bedtools intersect. When we call intersectBed without any parameters, we get the usage on STDERR. We can redirect STDERR to STDOUT, so that we can use a terminal pager:

intersectBed 2>&1 | less

Here is the entire usage:

Tool:    bedtools intersect (aka intersectBed)
Version: v2.20.1-4-gb877b35
Summary: Report overlaps between two feature files.

Usage:   bedtools intersect [OPTIONS] -a  -b 

Options: 
	-abam	The A input file is in BAM format.  Output will be BAM as well.

	-ubam	Write uncompressed BAM output. Default writes compressed BAM.

	-bed	When using BAM input (-abam), write output as BED. The default
		is to write output in BAM when using -abam.

	-wa	Write the original entry in A for each overlap.

	-wb	Write the original entry in B for each overlap.
		- Useful for knowing _what_ A overlaps. Restricted by -f and -r.

	-loj	Perform a "left outer join". That is, for each feature in A
		report each overlap with B.  If no overlaps are found, 
		report a NULL feature for B.

	-wo	Write the original A and B entries plus the number of base
		pairs of overlap between the two features.
		- Overlaps restricted by -f and -r.
		  Only A features with overlap are reported.

	-wao	Write the original A and B entries plus the number of base
		pairs of overlap between the two features.
		- Overlapping features restricted by -f and -r.
		  However, A features w/o overlap are also reported
		  with a NULL B feature and overlap = 0.

	-u	Write the original A entry _once_ if _any_ overlaps found in B.
		- In other words, just report the fact >=1 hit was found.
		- Overlaps restricted by -f and -r.

	-c	For each entry in A, report the number of overlaps with B.
		- Reports 0 for A entries that have no overlap with B.
		- Overlaps restricted by -f and -r.

	-v	Only report those entries in A that have _no overlaps_ with B.
		- Similar to "grep -v" (an homage).

	-f	Minimum overlap required as a fraction of A.
		- Default is 1E-9 (i.e., 1bp).
		- FLOAT (e.g. 0.50)

	-r	Require that the fraction overlap be reciprocal for A and B.
		- In other words, if -f is 0.90 and -r is used, this requires
		  that B overlap 90% of A and A _also_ overlaps 90% of B.

	-s	Require same strandedness.  That is, only report hits in B
		that overlap A on the _same_ strand.
		- By default, overlaps are reported without respect to strand.

	-S	Require different strandedness.  That is, only report hits in B
		that overlap A on the _opposite_ strand.
		- By default, overlaps are reported without respect to strand.

	-split	Treat "split" BAM or BED12 entries as distinct BED intervals.

	-sorted	Use the "chromsweep" algorithm for sorted (-k1,1 -k2,2n) input.

	-g	Provide a genome file to enforce consistent chromosome sort order
		across input files. Only applies when used with -sorted option.

	-header	Print the header from the A file prior to results.

	-nobuf	Disable buffered output. Using this option will cause each line
		of output to be printed as it is generated, rather than saved
		in a buffer. This will make printing large output files 
		noticeably slower, but can be useful in conjunction with
		other software tools and scripts that need to process one
		line of bedtools output at a time.

Notes: 
	(1) When a BAM file is used for the A file, the alignment is retained if overlaps exist,
	and exlcuded if an overlap cannot be found.  If multiple overlaps exist, they are not
	reported, as we are only testing for one or more overlaps.

Something to remember is that:

1.3.6 File B is loaded into memory
Whenever a BEDTool compares two files of features, the “B” file is loaded into memory. By contrast, the “A” file is processed line by line and compared with the features from B. Therefore to minimize memory usage, one should set the smaller of the two files as the B file.

Now, let's do some sanity checks. Let's test the -f parameter:

cat a.bed
#chr1    500 1000

cat b.bed
#chr1    750 1000

intersectBed -a a.bed -b b.bed -f 0.50
#chr1    750 1000

#if I increase b.bed, it should still return a result
cat a.bed
#chr1    500 1000

cat b.bed
#chr1    750 2000

intersectBed -a a.bed -b b.bed -f 0.50
#chr1    750 1000

#if I make a.bed longer, the -f 0.50 condition is not fulfilled
cat a.bed
#chr1    499 1000

cat b.bed
#chr1    750 2000

intersectBed -a a.bed -b b.bed -f 0.50
#returns nothing

So if I want to check whether a feature of interest has at least 50% of its length overlapping another feature, the feature of interest should be the -a BED file i.e. the first BED file.

Testing reciprocity (-r):

cat a.bed
#chr1    500 1000

cat b.bed
#chr1    750 1250

intersectBed -a a.bed -b b.bed -f 0.50 -r
#chr1    750 1000

cat a.bed
#chr1    500 1000

cat b.bed
#chr1    750 1251

intersectBed -a a.bed -b b.bed -f 0.50 -r
#returns nothing

intersectBed with BED12

I made this BED12 file:

chr1 1 1000 test 0 + 1 1000 0 2 100,100 0,899

So there are two blocks/exons at positions 1-100 and at 899-1000.

I then made this BAM file:

read1 0 chr1 50 37 21M * 0 0 GTATGCAGTATATACATGACT gggfcfcfcff^ffdgfgg]g XT:A:U NM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:21
read2 0 chr1 500 37 21M * 0 0 GTATGCAGTATATACATGACT gggfcfcfcff^ffdgfgg]g XT:A:U NM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:21

Then I ran:


#Version: 2.15.0
intersectBed -abam haha.bam -b test.bed -bed -wo


and got:

chr1 49 70 read1 37 + chr1 1 1000 test 0 + 1 1000 0 2 100,100 0,899 21
chr1 499 520 read2 37 + chr1 1 1000 test 0 + 1 1000 0 2 100,100 0,899 21

So intersectBed does not take into account the different blocks. This was posted by the developer (Aaron Quinlan) on SEQanswers:

(3) Native support for "blocked" BED features (aka BED12). Note that each block is not considered separately. BEDTools merely allows one to use BED12 files and the last 6 fields are "passed through" the tools.

So in order to intersect with blocks/exons, they should be separate entities in the BED file.

Other tips

#standard in / stream
cat a.bed | intersectBed -a stdin -b b.bed

#this is useful when you want to perform two intersect steps
intersectBed -a first.bed -b second.bed | intersectBed -a stdin -b third.bed

Lastly, below is what Sean Eddy has to say about BEDTools:

#BOSC2013 Eddy: Slightly insane level of commitment required to build a world class tool, eg now using excellent bedtools from @aaronquinlan

— Peter Cock (@pjacock) July 20, 2013