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:
module spider bedtools
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.
ssh user@login8.stampede.tacc.utexas.edu #if you are not already logged in! cds mkdir bedtools cd samtools cp yeast_pairedend_sort.mapped.q1.bam ../bedtools cd ../bedtools
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.
bedtools bamtofastq -i input.bam -fq output.fastq
Exercise 1: convert BAM to fastq and look at the quality scores
Bed file format: 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.
bedtools bamtobed -i input.bam > output.bed #output to a file bedtools bamtobed -i input.bam | more #output to more
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.
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.
bedtools coverage -a experiment.bed -b reference_file.bed
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:
cd bedtools #if you aren't already there cp /scratch/01786/awh394/core_ngs/day4_2015/sacCer3.chrom.sizes.bed .
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
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 [OPTIONS] -i experiment.bed > experiment.merge.bed
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!
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] -a <FILE> \ -b <FILE1, FILE2, ..., FILEN>
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:
cds mkdir intersect cd intersect cp /scratch/02423/nsabell/core_ngs/alignment/bam/human_rnaseq_bwa.bam . cp /scratch/02423/nsabell/core_ngs/alignment/bam/human_mirnaseq_hg19.bam .
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.
A little bit of filtering, using awk
As a final note, yesterday Nathan 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. 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:
cat file.bed | awk 'BEGIN{FS="\t";OFS="\t";}{if ($6 == '+'){print}}' > file.plusStrand.bed
- In the first section, we open the bed file of interest. Then we pipe that filestream to the awk program.
- 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.
- Generally, you can leave this section constant (if you are working with tab delimited files).
- The second section: {if ($6 == '+'){print}} is our selection and printing criteria.
- "$6" indicates column 6, and == indicates "equals" or "matches".
- The {print} command tells awk to print the whole line if the statement for column 6 evaluates to true.
- 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.
cat hg19_rnaseq_mirnaseq_intersect.bed | awk 'BEGIN{FS="\t";OFS="\t";}{if ($6 == "+"){print}}' | more
You could also insist on columns 6 and 12 both being the plus strand as such:
cat hg19_rnaseq_mirnaseq_intersect.bed | awk 'BEGIN{FS="\t";OFS="\t";}{if ($6 == "+" && $12 == "+"){print}}' | more