Skip to main content

1.6 BED format and BEDtools

BED format

Bed files are a common format for representing genomic intervals. They are tab-delimited text files that contain information about the location of features on a genome, such as genes, exons, or regulatory elements. Each line in a bed file represents a single feature and contains at least three required fields: chromosome, start position, and end position. Additional optional fields can include name, score, strand, and more.

For example, a simple bed file might look like this:

chr1    1000    5000
chr1    6000    9000
chr2    2000    4000

or like this as a 6-column bed file with strand information:

chr1    1000    5000    .   0   +
chr1    6000    9000    .   0   -
chr2    2000    4000    .   0   +

Or like this with names:

chr1    1000    5000    geneA   0   +
chr1    6000    9000    geneB   0   -
chr2    2000    4000    geneC   0   +

Using awk to view bed files

Awk is a powerful command-line tool for processing and analyzing text files, including bed files. You can use the awk command to view and manipulate bed files in the terminal.

As bedfiles are by definition tab-delimited the FS and OFS variables are set to tab characters to ensure proper tab-delimited output. However these can be changed to spaces or other delimiters (such as commas ‘,’ or lines ‘|’) as needed.

Columns are defined by $1, $2, $3 etc. where $1 is the first column, $2 the second column and so on. $0 refers to the entire line.

To extract specific columns from a bed file:

awk 'BEGIN{FS=OFS="\t"}{print $1, $2, $3}' file.bed
awk 'BEGIN{FS=OFS="\t"}{print $1, $2, $3}' file.bed > output.bed # Save to new file output.bed

This command will print the first three columns (chromosome, start, end) of the bedfile file.bed.

To sort a bed file by chromosome, start and end position:

awk 'BEGIN{FS=OFS="\t"}{print $0}' | "sort -k1,1 -k2,2n -k3,3n" file.bed

To filter bed entries based on a specific chromosome:

awk 'BEGIN{FS=OFS="\t"}($1=="chr1"){print $0}' file.bed

To remove duplicate entries from a bed file:

awk 'BEGIN{FS=OFS="\t"}!seen[$0]++' file.bed

To remove entries with a score below a certain threshold (e.g., 10):

awk 'BEGIN{FS=OFS="\t"}($5>=10){print $0}' file.bed

To skip the header line (if present) in a bed file:

awk 'NR>1 {print}' file.bed

To remove entries on a specific chromosome (e.g., chrM):

awk 'BEGIN{FS=OFS="\t"}($1!="chrM"){print $0}' file.bed

BEDtools

Bedtools is a powerful suite of command-line tools for manipulating and analyzing bed files and other genomic interval formats. It allows users to perform a wide range of operations on genomic intervals, such as intersecting, merging, subtracting, and shifting intervals.

It is loaded on the CCB cluster by:

module load bedtools

The documentation for bedtools is incredibily informative and gives the full breakdown of what each function can be used for

BEDtools Documentation Pages.

Use case to make a genome wide promoter list from a gene annotation file

Step 1 - Download the necessary files

The UCSC server hosts a variety of genome annotation files that can be downloaded via wget. This includes annotations of genes, transcripts, and other genomic features for a wide range of organisms, including human and mouse.

Here we download the RefGene annotation file for mouse genome version mm9, as well as the chromosome sizes file.

wget -c -O RefGene.mm9.txt.gz http://hgdownload.soe.ucsc.edu/goldenPath/mm9/database/refGene.txt.gz
wget http://hgdownload.cse.ucsc.edu/goldenPath/mm9/database/chromInfo.txt.gz
gunzip chromInfo.txt.gz
cat chromInfo.txt

Step 2 - Create a bed file of promoters from the refGene annotation file

The refGene file contains a variety of information about each gene, including its chromosome, start and end positions, strand, and more. To look at the first few lines of the refGene file:

gunzip RefGene.mm9.txt.gz
head -2 RefGene.mm9.txt
#bin    name    chrom   strand  txStart txEnd   cdsStart        cdsEnd  exonCount       exonStarts      exonEnds        score   name2   cdsStartStat    cdsEndStat      exonFrames
0       NM_028778       chr1    +       134212701       134230065       134212806       134228958       7       134212701,134221529,134224273,134224707,134226534,134227135,134227897,  134213049,134221650,134224425,134224773,134226654,134227268,134230065,  0       Nuak2   cmpl    cmpl    0,0,1,0,0,0,1,

We can use awk to extract the relevant columns and create a bed file, in this case we extract chromosome (column 3), txStart (column 5), txEnd (column 6), name (column 2 and 13 combined), score (set to 0) and strand (column 4):

gunzip RefGene.mm9.txt.gz
awk 'BEGIN{FS=OFS="\t"}{print $3,$5,$6,$2"|"$13,0,$4}' RefGene.mm9.txt > RefGene.mm9.bed
head -3 RefGene.mm9.bed
#chrom  txStart txEnd   name|name2  0   strand
chr1    134212701   134230065   NM_028778|Nuak2 0   +
chr1    134212701   134230065   NM_001195025|Nuak2  0   +
awk 'BEGIN{FS=OFS="\t"}{print $0}' RefGene.mm9.bed | sort -k1,1 -k2,2n -k3,3n > RefGene.mm9.sorted.bed

Step 3 - Extract TSS positions

The transcription start site (TSS) is the position where transcription of a gene begins. For genes on the positive strand, the TSS is at the start position, while for genes on the negative strand, the TSS is at the end position. We can use awk to extract the TSS positions from the bed file:

awk 'BEGIN{FS=OFS="\t"}($6=="+"){print $1,$2,$2,$4,0,$6} ($6=="-"){print $1,$3,$3,$4,0,$6}' RefGene.mm9.sorted.bed > RefGene.mm9.TSS.bed
head -2 RefGene.mm9.TSS.bed
chr1    134212701       134212701       NM_028778|Nuak2 0       +
chr1    134212701       134212701       NM_001195025|Nuak2      0       +

Step 4 - Sort and remove duplicates

Finally, we sort the promoter bed file and remove any duplicate entries:

awk 'BEGIN{FS=OFS="\t"}{print $0}' RefGene.mm9.TSS.bed | sort -k1,1 -k2,2n -k3,3n > RefGene.mm9.TSS.sorted.bed
awk 'BEGIN{FS=OFS="\t"}!seen[$0]++' RefGene.mm9.TSS.sorted.bed > RefGene.mm9.TSS.sorted.unique.bed
head -4 RefGene.mm9.TSS.sorted.unique.bed
chr1    3661579 3661579 NM_001011874|Xkr4       0       -
chr1    4350384 4350384 NM_001370921|Rp1        0       -
chr1    4350395 4350395 NM_011283|Rp1   0       -
chr1    4399322 4399322 NM_001195662|Rp1        0       -
  • This is a limitation of this approach - some genes have multiple TSS positions due to alternative promoters or transcription start sites.

A way of removing the duplicates while retaining only one TSS per gene would be to use bedtools groupby on the name field, taking the min or max of the start position depending on strand or by filtering for unique values in the name field. This could be explored as an excercise.

Step 5 - Create rough promoter regions

Promoter regions are typically defined as a certain number of base pairs upstream and downstream of the transcription start site (TSS). Here we define promoters as 2000 base pairs upstream and 500 base pairs downstream of the TSS, taking into account the strand information. We can use bedtools to create these promoter regions:

bedtools slop -i RefGene.mm9.TSS.sorted.unique.bed -g chromInfo.txt -l 2000 -r 500 -s > RefGene.mm9.2000_500_unique_promoters.bed

The final output file RefGene.mm9.2000_500_promoters.sorted.unique.bed contains a list of unique promoter regions for all genes in the refGene annotation file, defined as 2000 base pairs upstream and 500 base pairs downstream of the TSS.

chr1    3661079 3663579 NM_001011874|Xkr4   0   -
chr1    4349884 4352384 NM_001370921|Rp1    0   -
chr1    4486935 4489435 NM_011441|Sox17 0   -
chr1    4775307 4777807 NM_025300|Mrpl15    0   -
chr1    4795903 4798403 NM_001355712|Lypla1 0   +

Advanced BEDtools functions - Intersecting and merging intervals

We may want to find overlapping regions between two bed files, such as finding peaks from a ChIP-seq experiment that overlap with promoter regions. In this case this is a ChIP.peaks.bed file. Peak calling is covered in Chapter 4.

We can use the bedtools intersect command to achieve this:

bedtools intersect -a ChIP.peaks.bed -b RefGene.mm9.2000_500_promoters.sorted.unique.bed -wa -wb > mm9.Promoters_with_ChIP_signal.bed

This command will output a new bed file containing the overlapping regions between the peaks and promoter regions, along with the original entries from both files.

Peak Visualization

Converting Peaks to BigBed

To load peak files into genome browsers like UCSC or IGV, it is often useful to convert them to bigBed format for efficient viewing.

# Sort peaks
sort -k1,1 -k2,2n ChIP.peaks.bed > ChIP.peaks_sorted.bed

# Convert to bigBed
bedToBigBed ChIP.peaks_sorted chrom.sizes ChIP.peaks.bb
# You can get chrom.sizes file from UCSC genome browser downloads page using wget 
# or from the /databank/igenomes/ folder on the CCB

Additional Resources

I wuld again recommend the BEDtools Documentation Pages for further reading and exploration of the many functions available in BEDtools.