Bedtools: Analyzing your aligned experiment

Use bedtools coverage to create a signal track

A signal track is a bedGraph (BED3+) file with an entry for every base in a defined set of regions (see https://genome.ucsc.edu/goldenpath/help/bedgraph.html). bedGraph files can be visualized in the Broad's IGV (Integrative Genomics Viewer) application (https://software.broadinstitute.org/software/igv/download) or in the UCSC Genome Browser (https://genome.ucsc.edu/).

The bedtools coverage function (https://bedtools.readthedocs.io/en/latest/content/tools/coverage.html), with the -d (per-base depth) option produces output that can be made into a bedGraph. Here we'll analyze the per-base coverage of yeast RNAseq reads in our merged yeast gene regions.

Make sure you're in an idev session, then prepare a directory for this exercise.

Prepare for bedtools coverage

Then calling bedtools coverage is easy. The "A" file will be our gene regions, and the "B" file will be the yeast RNAseq reads. We also use the -d (per-base depth) and -s (force "strandedness") options.

It will complain a bit because our genes file includes the yeast plasmid "2-micron" but the RNAseq BAM doesn't include that contig. We'll ignore that warning.

The bedtools coverage output is a bit strange. It lists each region in the A file, followed by information from the B reads. Here the column order will be gene_chrom gene_start gene_end gene_name gene_score gene_strand offset_in_the_gene_region read_overlap count.

Let's look at coverage for gene YAL067C:

Will look like this:

A proper bedGraph file has only 4 columns: chrom start end value and does not need to include positions with 0 reads, so we'll convert the bedtools coverage output to bedGraph using awk. We re-sort the output so that plus and minus strand positions are adjacent.

While we probably could consider this file to have bedGraph format, it's preferable to combine adjacent per-base coordinates with the same count into larger regions, e.g.

Here's some awk to do this:

Make sure the total counts match!

Now our yeast_mrna.gene_coverage.bedGraph file is a proper bedGraph, whose first lines look like this:

x

A brief introduction to bedtools

Now that we have a BAM file with only the reads we want included, we can do some more sophisticated analysis using  bedtools.  Bedtools changes from version to version, and here we are using version 2.22, the newest version, and what is currently installed on stampede.  You can check what versions of bedtools are installed by using the following command on stampede:

First, log on to the login8 node on stampede and make a directory in scratch called bedtools in your scratch folder.  Then copy your filtered BAM file from the samtools section into this folder.

If you were unable to make the filtered and sorted BAM file from the previous section, you can copy it over from my scratch directory:

bedtools bamtofastq: converting a BAM file to a fastq file

Sometimes, especially when working with external data, we need to go from a BAM file back to a fastq file.  This can be useful for re-aligning reads using a different aligner, different settings on the original aligner used.  It can also be useful for extracting the sequence of interesting regions of the genome after you have manipulated your BAM file.

For this exercise, you'll be using bamtofastq.  This function takes an aligned BAM file as input and outputs a fastq format file.  You can use the options if you have paired end data to output R1 and R2 reads for your fastq file.  This type of function is especially useful if you need to to  analyze sequences after you've compared several BAM or bed files.

Exercise 1: convert BAM to fastq and look at the quality scores

bedtools bamtobed: converting a BAM file into a bed file

While it's useful to be able to look at the fastq file, many analyses will be easiest to perform in bed format.  Bed format is a simple tab delimited format that designates various properties about segments of the genome, defined by the chromosome, start coordinates and end coordinates.  Bedtools provides a simple utility to convert BAM files over into bed files, termed bamtobed.

Note that the output will be piped to standard out unless you redirect to a program (head, more, less) or a file (output.bed).  Now we'll put this example to use and convert our filtered BAM file from the samtools section into a bed file.

Exercise 2: Convert the filtered yeast paired end BAM to bed using bamtobed, look at your file in more, and find the number of lines in the file

Hint: direct the output to a file first, then use more to look at the converted file visually; use ctrl+c to quit more.

See also: bedtools bedtobam, if you need to get back to a bam file from a bed file (some programs take bam files as input).  Documentation here.

bedtools coverage: how much of the genome does my data cover?

One way of characterizing data is to understand what percentage of the genome your data covers.  What type of experiment you performed should affect the coverage of your data.  A ChIP-seq experiment will cover binding sites, and a RNA-seq experiment will cover expressed transcripts.  Bedtools coverage allows you to compare one bed file to another and compute the breadth and depth of coverage. 

The resulting output will contain several additional columns which summarize this information:

For this exercise, we'll use a bed file that summarizes the S. cerevisiae genome, version 3 (aka sacCer3).  For this class, I've made a bed file out of the genome, using the file sacCer3.chrom.sizes.  First go and copy the file over from my scratch directory:

Now use bedtools coverage to find the coverage of the file output.bed over the sacCer3 genome and examine the output coverage.

Exercise 3: Find the coverage of your bed file over the sacCer3 genome

It's worth noting that just computing coverage over the genome isn't the most useful thing, but you might compute coverage over a set of genes or regions of interest.  Coverage is really useful coupled with intersect or subtract as well.

bedtools merge: collapsing bookended elements (or elements within a certain distance)

When we originally examined the bed files produced from our BAM file, we can see many reads that overlap over the same interval.  While this level of detail is useful, for some analyses, we can collapse each read into a single line, and indicate how many reads occurred over that genomic interval.  We can accomplish this using bedtools merge.

Bedtools merge also directs the output to standard out, to make sure to point the output to a file or a program.  While we haven't discussed the options for each bedtools function in detail, here they are very important.  Many of the options define what to do with each column (-c) of the output (-o).  This defines what type of operation to perform on each column, and in what order to output the columns.  Standard bed6 format is chrom, start, stop, name, score, strand and controlling column operations allows you to control what to put into each column of output.  The valid operations defined by the -o operation are as follows:

For this exercise, we'll be summing the number of reads over a region to get a score column, using distinct to choose a name, and using distinct again to keep track of the strand.  For the -c options, define which columns to operate on, in the order you want the output.  In this case, to keep the standard bed format, we'll list as -c 5,4,6 and -o count_distinct,sum,distinct, to keep the proper order of name, score, strand.  Strandedness can also be controlled with merge, using the -s option. 

Exercise 4: Use bedtools merge to merge an experiment, look at the output and see how many lines there are in the file.

Hint: make sure to remove whitespace between lists for the -c and -o options!

Note the change in column order in the first set of commands.  We can use awk like this to change the column order, either piped in the original command or after the fact:

using awk for column reordering

bedtools intersect: identifying where two experiments overlap (or don't overlap)

One useful way to compare two experiments (especially biological replicates, or similar experiments in two yeast strains/cell lines/mouse strains) is to compare where reads in one experiment overlap with reads in another experiment.  Bedtools offers a simple way to do this using the intersect function.

bedtools intersect options

The intersect function has many options that control how to report the intersection.  We'll be focusing on just a few of these options, listed below.

-a and -b indicate what files to intersect.  in -b, you can specify one, or several files to intersect with the file specified in -a.

In this section, we'll intersect two human experiments - one from sequencing RNA, and one from sequencing micro RNA.  Copy these files over to your directory:

copy some files over to intersect

Before we can intersect these files, we need to perform the pipeline we used in samtools to index, sort and filter the files, and bedtools to convert from BAM over to bed, then collapse down the files using merge.  Below is a little workflow to help you through it on the files you just copied above.

If we run low on time, you can copy the merged bed files over from my directory on scratch:

Exercise 5: Intersect two experiments using intersect and examine the output

Using the options we've specified (-wo) the resulting file will have entries for file A, file B and the number of base pairs overlap between the feature in A and the features in B, but we'll only retain lines where there is an overlap between A and B.  We could also use the -v option to only contain areas with NO intersection, or control the intersections with -f and -r options.  Bedtools intersect is a powerful tool, and it's always a good idea to ask "what is this code going to do?" while you're testing analysis workflows.  It can be very useful to pipe your output to more when you are unsure of the output of a command, as such:

pipe-ing output to more

bedtools closest: when you want to know how far your regions are from a test set

The manual page for bedtools closest has a really nice image of how closest behaves with overlapping options.  Bedtools closest first looks for any overlaps of B with A, if it finds an overlap, the overlap in B with the highest proportional overlap with A is reported.  If there are no overlaps, then it looks for the closest genomic feature proximal to A (using distance from the start or end of A to do this).

bedtools intersect options

Much like bedtools intersect, bedtools closest takes an A file and a series of B files.  So if you wanted to determine the distance of your regions of interest to several different classes of genes, bedtools closest would be a useful tool for that analysis. 

In this section, we'll intersect the human_rnaseq_bwa_sort.mapped.q1.merge.bed file with some protein coding genes from Gencode (hg19).  First go copy a couple files from my scratch directory:

copy some gencode files over

Exercise 6: Identify the closest protein coding genes (on the same strand) for the human rnaseq file using closest, then sort by the distance column (largest, then smallest distance first).

This is a nice way to examine your reads over annotated protein-coding genes.  Note the strand specificity - only reads on the correct strand will be reported when there is a + strand gene and a - strand gene over the same location.

bedtools subtract:  removing features from your bed file

Bedtools subtract takes an A file and a B file, then searches for features in B that overlap A.  When/if these features are identified, the overlapping portion is removed from A and the remaining portion of A is reported.  If a feature in B overlaps all of a feature in A, that feature will not be reported.

bedtools subtract options

Note that bedtools subtract is performed on two files, and unlike some of the other utilities we've used, you can't use multiple B features here.  However, you can use cat to join together features you'd like to subtract from your A file.

Let's do a little set-up for the next exercise:

copy some gencode files over

Exercise 7:  remove the protein-coding genes from a gencode list of genes using subtract, then give a count of the non-protein-coding gene entries

This allows you to identify which gene regions are not protein coding, and are likely pseudogenes, but could also be miRNAs, snRNAs or other genes that aren't translated into a peptide sequence.

While the above example is not super useful in all cases, one might use the above workflow to remove genes that aren't of interest from a larger set.

A little bit of filtering, using awk

As a final note, yesterday we taught you about using a lot of unix utilities, including uniq, sort and cut.  One last utility I'd like to add, that is very useful for manipulating these types of tab delimited files, is awk.  Awk isn't a command, but rather a little text manipulation language in it's own right (which we briefly used above to rearrange the columns in a file).  While awk can be used to do many different things, here we'll primarily use it to sort tab delimited files based on the values present in those files.  That is useful to filter your files for entries on a given chromosome, or greater than/less than a given score.  If your dataset is large, this type of filtering can be invaluable!  Below is an example of a simple awk script:

a simple awk script

1. In the first section, we open the bed file of interest.  Then we pipe that filestream to the awk program. 

2. The section: BEGIN{FS="\t";OFS="\t";} tells awk to begin a filter, the input file is tab delimited, and the output file is also tab delimited. 

   1. Generally, you can leave this section constant (if you are working with tab delimited files). 

3. The second section: {if ($6 == '+'){print}} is our selection and printing criteria. 

   1. "$6" indicates column 6, and == indicates "equals" or "matches". 

4. The {print} command tells awk to print the whole line if the statement for column 6 evaluates to true. 

5. Thus, the output file only contains those lines which satisfy the criteria in the selection statement.

We can do this filtering on the hg19_rnaseq_mirnaseq_intersect.bed file we just created using bedtools intersect.

You could also insist on columns 6 and 12 both being the plus strand as such: