Analysis using BEDTools

The BED format

BED (Browser Extensible Data) format is a simple text format for location-oriented data (genomic regions) developed to support UCSC Genome Browser tracks. Standard BED files have 3 to 6 Tab-separated columns, although up to 12 columns are defined. (Read more about the UCSC Genome Browser's official BED format.)

Memorize the 6 main BED fields

These 6 BED fields are so important that you should memorize them. Keep repeating "chrom, start, end, name, score, strand" until the words trip off your tongue

chrom (required) – string naming the chromosome or othre contig
start (required) – the 0-based start position of the region
end (required) – the 1-based end position of the region
name (optional) – an arbitrary string describing the region
- for BED files loaded as UCSC Genome Browser tracks, this text is displayed above the region
score (optional) – an integer score for the region
- for BED files to be loaded as UCSC Genome Browser tracks, this should be a number between 0 and 1000, higher = "better"
- for non-GenBrowse BED files, this can be any integer value (e.g. the length of the region)
strand (optional) - a single character describing the region's strand
- + – plus (Watson) strand region
- - – minus (Crick) strand region
- . – region is not associated with a strand (e.g. a transcription factor binding region)

Important rules for BED format:

The number of fields per line must be consistent throughout any single BED file
- e.g. they must all have 3 fields or all have 6 fields
The first base on a contig is numbered 0
- versus 1 for BAM file positions
- so the a BED start of 99 is actually the 100th base on the contig
- but end positions are 1-based
  - so a BED end of 200 is the 200th base on the contig
- the length of a BED region is end - start
  - not end - start + 1, as it would be if both coordinates with 0-based or both 1-based
- this difference is the single greatest source of errors dealing with BED files!

Note that the UCSC Genome Browser also defines many BED-like data formats (e.g. bedGraph, narrowPeak, tagAlign and RNA elements formats). See supported UCSC Genome Browser data formats for more information and examples.

In addition to standard-format BED files, one can create custom BED files that have at least 3 of the standard fields, followed by any number of custom fields. For example:

A BED3+ file contains the 3 required BED fields, followed by some number of user-defined columns (all records with the same number)
A BED6+ file contains the 3 required BED fields, followed by some number of user-defined columns

As we will see, BEDTools functions require BED3+ input files, or BED6+ if strand-specific operations are requested.

BEDTools overview

The BEDTools suite is a set of utilities for manipulating BED and BAM files. We call it the "Swiss army knife" for genomic region analyses because its sub-commands are so numerous and versatile. Some of the most common bedtools operations perform set-theory functions on regions: intersection (intersect), union (merge), set difference (subtract) – but there are many others. The table below lists some of the most useful sub-commands along with applicable use cases.

Sub-command	Description	Use case(s)
bamtobed	Convert BAM files to BED format.	You want to have the contig, start, end, and strand information for each mapped alignment record in separate fields. Recall that the strand is encoded in a BAM flag (0x10) and the exact end coordinate requires parsing the CIGAR string.
bamtofastq	Extract FASTQ sequences from BAM alignment records.	You have downloaded a BAM file from a public database, but it was not aligned against the reference version you want to use (e.g. it is hg18 and you want an hg38 alignment). To re-process, you need to start with the original FASTQ sequences.
getfasta	Get FASTA entries corresponding to regions.	Preparation to run motif analysis, which requires the original sequences, on a set of regions of interest. You must provide FASTA file(s) for the genome/reference used for alignment (e.g. the FASTA file used to build the aligner index).
coverage	Compute genome-wide coverage of your regsions, or generate per-base genome-wide signal trace.	You have performed a WGS (whole genome sequencing) experiment and want to know if has resulted in the desired coverage depth. Calculate what proportion of the (known) transcriptome is covered by your RNA-seq alignments. Provide the transcript regions as a BED or GFF/GTF file. Produce a per-base genome-wide signal (in bedGraph format) for a ChIP-seq experiment. After conversion to binary bigWig format, such tracks can be configured in the UCSC Genome Browser as custom tracks.
multicov	Count overlaps between one or more BAM or BED files and a set of regions of interest.	Count RNA-seq alignments that overlap a set of genes of interest. While this task is usually done with a specialized RNA-seq quantification tool (e.g. DESeq2 or EdgeR), bedtools multicov can provide a quick estimate, e.g. for QC purposes.
merge	Combine a set of possibly-overlapping regions into a single set of non-overlapping regions.	Collapse overlapping gene annotations into per-strand non-overlapping regions before counting. If this is not done, the source regions will potentially be counted multiple times, once for each (overlapping) target region it intersects.
subtract	Remove unwanted regions.	Remove rRNA gene regions from a merged gene annotations file before counting.
intersect	Determine the overlap between two sets of regions.	Similar to multicov, but can also report (not just count) the overlapping regions.
closest	Find the genomic features nearest to a set of regions.	For a set of significant ChIP-seq transcription factor (TF) binding regions ("peaks") that have been identified, determine nearby genes that may be targets of TF regulation.

We will explore a few of these functions in our exercises.

BEDTools versions

BEDTools is under active development and is always being refined and extended. Unfortunately, sometimes changes are made that are incompatible with previous BEDTools versions. For example, a major change to the way bedtool merge functions was made after bedtools v2.17.0.

So it is important to know which version of BEDTools you are using, and read the documentation carefully to see if changes have been made since your version. When you run the code below, you should see that the version on TACC is bedtools v2.25.0.

copy and paste exercise files

module load bedtools
bedtools --version

Input format considerations

Most BEDTools functions now accept either BAM or BED files as input.
- BED format files must be BED3+, or BED6+ if strand-specific operations are requested.
When comparing against a set of regions, those regions are usually supplied in either BED or GFF/GTF.
All text-format input files (BED, GFF, VCF) must use Unix line endings (linefeed only).

The most important thing to remember about comparing regions using BEDTools, is that all region files must share the same set contig names and be based on the same reference! For example, if an alignment was performed against a human GRCh38 reference genome from Gencode, annotations from the corresponding GTF annotations.

Exercises

Setup

copy and paste exercise files

cds
mkdir bedtools_practice
cd bedtools_practice
cp /corral-repl/utexas/BioITeam/core_ngs_tools/yeast_stuff/saccharomyces_cerevisiae_R64-1-1_20110208.gff .
cp /corral-repl/utexas/BioITeam/core_ngs_tools/yeast_stuff/*sort.bam .
cp /corral-repl/utexas/BioITeam/core_ngs_tools/yeast_stuff/*.bai .
cp /corral-repl/utexas/BioITeam/core_ngs_tools/human_stuff/bedtools_sample/*bed .

module swap intel gcc/4.7.1
module load bedtools

Q: Convert saccharomyces_cerevisiae_R64-1-1_20110208.gff into a 3 column bed file that includes 'gene' feature.

The gff file is under /corral-repl/utexas/BioITeam/core_ngs_tools/yeast_stuff

Answer

conversion with awk

cat saccharomyces_cerevisiae_R64-1-1_20110208.gff | awk 'BEGIN{FS="\t";OFS="\t"} $3=="gene" {print $1,$4-1,$5}' > sc_gene_3c.bed

Q: Convert saccharomyces_cerevisiae_R64-1-1_20110208.gff into a 4 column bed file that includes 'gene' feature (4th column has gene IDs)

Answer

conversion with awk

cat /corral-repl/utexas/BioITeam/core_ngs_tools/yeast_stuff/saccharomyces_cerevisiae_R64-1-1_20110208.gff | \
awk 'BEGIN{FS="\t";OFS="\t"} $3=="gene" {split($9,x,";") ; $9=substr(x[1],4) ; print $1,$4-1,$5,$9}' > sc_gene.bed

BEDTools merge

BEDTools multicov

BEDTools is a great utility for comparing genomic features in BAM, BED, VCF, and GFF formats. The documentation is well written in great detail. Here, we will practice three commonly used sub-commands: multicov, merge, and intersect.

We are going to use it to count the number of reads that map to each gene in the genome. Load the module and check out the help for bedtools and the multicov specific command that we are going to use:

bedtools
bedtools multicov

The multicov command takes a feature file (GFF/BED/VCF) and counts how many reads are in certain regions from many input files. By default it counts how many reads overlap the feature on either strand, but it can be made specific with the -s option.

Note: Remember that the chromosome names in your gff file should match the way the chromosomes are named in the reference fasta file used in the mapping step. For example, if BAM file used for mapping contains chr1, chrX etc, the GFF file must also call the chromosomes as chr1, chrX and so on.

In order to use the bedtools command on our data, do the following:

bedtools multicov -bams yeast_chip_sort.bam yeast_chip2_sort.bam -bed sc_gene.bed > gene_counts.txt

Then take a peek at the data...

head gene_counts.txt

Page tree