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
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 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 occured 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 distinct,sum,distinct, to keep the proper order of name, score, strand.
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 miRNA. 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.bam .
As above, let's convert the high quality reads from these files over into a filtered bam file, then a bed file.
Exercise 5: Intersect two experiments using intersect and examine the output
Using the options we've specified (-wao) 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.
