Bedtools: Analyzing your aligned experiment
- Amelia W Hall
- Anna Battenhouse
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.
idev -m 120 -N 1 -A OTH21164 -r CoreNGSday5 module load biocontainers module load bedtools mkdir -p $SCRATCH/core_ngs/bedtools_coverage cp $CORENGS/catchup/bedtools_merge/merged*bed $SCRATCH/core_ngs/bedtools_coverage/ cp $CORENGS/yeast_rnaseq/yeast_mrna.sort.filt.bam* $SCRATCH/core_ngs/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.
cds; cd core_ngs/bedtools_coverage bedtools coverage -s -d -a merged.good.sc_genes.bed -b yeast_mrna.sort.filt.bam > yeast_mrna.gene_coverage.txt wc -l yeast_mrna.gene_coverage.txt # 8,829,317 lines!
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:
cat yeast_mrna.gene_coverage.txt | grep -P 'YAL067C' | head -50
Will look like this:
chrI 7234 9016 YAL067C 1 - 1 0 chrI 7234 9016 YAL067C 1 - 2 0 chrI 7234 9016 YAL067C 1 - 3 0 chrI 7234 9016 YAL067C 1 - 4 0 chrI 7234 9016 YAL067C 1 - 5 0 chrI 7234 9016 YAL067C 1 - 6 0 chrI 7234 9016 YAL067C 1 - 7 0 chrI 7234 9016 YAL067C 1 - 8 0 chrI 7234 9016 YAL067C 1 - 9 0 chrI 7234 9016 YAL067C 1 - 10 0 chrI 7234 9016 YAL067C 1 - 11 0 chrI 7234 9016 YAL067C 1 - 12 0 chrI 7234 9016 YAL067C 1 - 13 0 chrI 7234 9016 YAL067C 1 - 14 0 chrI 7234 9016 YAL067C 1 - 15 0 chrI 7234 9016 YAL067C 1 - 16 0 chrI 7234 9016 YAL067C 1 - 17 1 chrI 7234 9016 YAL067C 1 - 18 1 chrI 7234 9016 YAL067C 1 - 19 1 chrI 7234 9016 YAL067C 1 - 20 1 chrI 7234 9016 YAL067C 1 - 21 1 chrI 7234 9016 YAL067C 1 - 22 1 chrI 7234 9016 YAL067C 1 - 23 1 chrI 7234 9016 YAL067C 1 - 24 1 chrI 7234 9016 YAL067C 1 - 25 1 chrI 7234 9016 YAL067C 1 - 26 1 chrI 7234 9016 YAL067C 1 - 27 1 chrI 7234 9016 YAL067C 1 - 28 1 chrI 7234 9016 YAL067C 1 - 29 1 chrI 7234 9016 YAL067C 1 - 30 1 chrI 7234 9016 YAL067C 1 - 31 1 chrI 7234 9016 YAL067C 1 - 32 1 chrI 7234 9016 YAL067C 1 - 33 1 chrI 7234 9016 YAL067C 1 - 34 1 chrI 7234 9016 YAL067C 1 - 35 1 chrI 7234 9016 YAL067C 1 - 36 1 chrI 7234 9016 YAL067C 1 - 37 1 chrI 7234 9016 YAL067C 1 - 38 2 chrI 7234 9016 YAL067C 1 - 39 2 chrI 7234 9016 YAL067C 1 - 40 2 chrI 7234 9016 YAL067C 1 - 41 3 chrI 7234 9016 YAL067C 1 - 42 3 chrI 7234 9016 YAL067C 1 - 43 3 chrI 7234 9016 YAL067C 1 - 44 3 chrI 7234 9016 YAL067C 1 - 45 4 chrI 7234 9016 YAL067C 1 - 46 4 chrI 7234 9016 YAL067C 1 - 47 4 chrI 7234 9016 YAL067C 1 - 48 4 chrI 7234 9016 YAL067C 1 - 49 4 chrI 7234 9016 YAL067C 1 - 50 4
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.
cat yeast_mrna.gene_coverage.txt | awk '
BEGIN{FS=OFS="\t"}
{if ($8>0) {print $1,$2-1+$7,$2+$7,$8}}' | \
sort -k1,1 -k2,2n -k3,3n > yeast_mrna.gene_coverage.almost.bedGraph
wc -l yeast_mrna.gene_coverage.almost.bedGraph # 5,710,186 -- better, but still big
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.
# per-base counts chrI 7271 7272 2 chrI 7272 7273 2 chrI 7273 7274 2 chrI 7274 7275 3 chrI 7275 7276 3 chrI 7276 7277 3 chrI 7277 7278 3 # corresponding region counts chrI 7271 7274 6 chrI 7274 7278 12
Here's some awk to do this:
cat yeast_mrna.gene_coverage.almost.bedGraph | awk '
BEGIN{FS=OFS="\t"; chr=""; start=-1; end=-1; count=0}
{if (chr != $1) { # new contig; finish previous
if (count > 0) { print chr,start,end,count }
chr=$1; start=$2; end=$3; count=$4
} else if (($2==end || $2==end+1) && ($4==count)) { # same or adjacent position with same count
end=$3;
} else { # new region on same contig; finish prev
if (count > 0) { print chr,start,end,count}
start=$2; end=$3; count=$4
}
}
END{ # finish last
if (count > 0) { print chr,start,end,count }
}' > yeast_mrna.gene_coverage.bedGraph
wc -l yeast_mrna.gene_coverage.bedGraph # 1,048,510 -- much better!
Make sure the total counts match!
cat yeast_mrna.gene_coverage.txt | awk '
BEGIN{tot=0}{tot=tot+$8}END{print tot}' # should be 86703686
cat yeast_mrna.gene_coverage.almost.bed | awk '
BEGIN{tot=0}{tot=tot+$4}END{print tot}' # should also be 86703686
cat yeast_mrna.gene_coverage.bedGraph | awk '
BEGIN{tot=0}{tot=tot+$4*($3-$2)}END{print tot}' # should also be 86703686
Now our yeast_mrna.gene_coverage.bedGraph file is a proper bedGraph, whose first lines look like this:
chrI 7250 7271 1 chrI 7271 7274 2 chrI 7274 7278 3 chrI 7278 7310 4 chrI 7310 7317 3 chrI 7317 7349 2 chrI 7349 7353 1 chrI 7500 7556 1 chrI 8851 8891 1 chrI 11919 11951 1
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:
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! cd $SCRATCH/core_ngs 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
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.
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.
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.
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.test/bedtools/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
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 [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!
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:
#after the creation of the first file
cat yeast_pairedend_sort.mapped.q1.merge.bed | awk '{print $1 "\t" $2 "\t" $3 "\t" $5 "\t" $6 "\t" $4}' > yeast_pairedend_sort.mapped.q1.merge.reorder.bed
#piped in-line
bedtools merge -s -c 4,5 -o count_distinct,sum -i yeast_pairedend_sort.mapped.q1.bed | awk '{print $1 "\t" $2 "\t" $3 "\t" $5 "\t" $6 "\t" $4}' > yeast_pairedend_sort.mapped.q1.merge.bed
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:
cd $SCRATCH/core_ngs/ mkdir intersect cd intersect cp /corral-repl/utexas/BioITeam/core_ngs_tools/alignment/bam/human_mirnaseq_hg19.bam . cp /corral-repl/utexas/BioITeam/core_ngs_tools/alignment/bam/human_rnaseq_bwa.bam . ls -lah
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:
bedtools intersect -wo -a human_rnaseq_bwa_sort.mapped.q1.merge.bed -b human_mirnaseq_hg19_sort.mapped.q1.merge.bed | 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 closest [OPTIONS] -a <FILE> \
-b <FILE1, FILE2, ..., FILEN>
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:
cd $SCRATCH/core_ngs mkdir closest cd closest cp /scratch/01786/awh394/core_ngs.test/closest/gencode.v19.proteincoding.genes.sort.merge.final . cp ../intersect/human_rnaseq_bwa_sort.mapped.q1.merge.bed . #or: cp /scratch/01786/awh394/core_ngs/closest/human_rnaseq_bwa_sort.mapped.q1.merge.bed .
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] -a <BED/GFF/VCF> -b <BED/GFF/VCF>
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:
cd $SCRATCH/core_ngs mkdir subtract cd subtract cp /scratch/01786/awh394/core_ngs.test/closest/gencode.v19.proteincoding.genes.sort.merge.final . cp /scratch/01786/awh394/core_ngs.test/closest/gencode.v19.genes.sort.merge.final .
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:
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.
cd $SCRATCH/core_ngs/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:
cd $SCRATCH/core_ngs/intersect/
cat hg19_rnaseq_mirnaseq_intersect.bed | awk 'BEGIN{FS="\t";OFS="\t";}{if ($6 == "+" && $12 == "+"){print}}' | more
Welcome to the University Wiki Service! Please use your IID (yourEID@eid.utexas.edu) when prompted for your email address during login or click here to enter your EID. If you are experiencing any issues loading content on pages, please try these steps to clear your browser cache.