I’m calling variants from exome sequencing data and I need to evaluate the efficiency of the capture and the coverage along the target regions.

This sounds like a great use case for bedtools, your swiss-army knife for genomic arithmetic and interval manipulation. I’m lucky enough to be able to walk down the hall and bug Aaron Quinlan, creator of bedtools, whenever I have a “how do I do X with bedtools?” question (which happens frequently).

As open-source bioinformatics software documentation goes, bedtools’ documentation is top-notch. In addition, Aaron recently pointed out a work-in-progress bedtools cookbook that he’s putting together, giving code examples for both typical and clever uses of bedtools.

Getting back to my exome data, one way to visualize this is to plot the cumulative distribution describing the fraction of targeted bases that were covered by >10 reads, >20 reads, >80 reads, etc. For example, covering 90% of the target region at 20X coverage may be one metric to assess your ability to reliably detect heterozygotes. Luckily for me, there’s a bedtools protocol for that.

The basic idea is that for each sample, you’re using bedtools coverage to read in both a bam file containing your read alignments and a bed file containing your target capture regions (for example, you can download NimbleGen’s V3 exome capture regions here). The -hist option outputs a histogram of coverage for each feature in the BED file as well as a summary histogram for all of the features in the BED file. Below I’m using GNU parallel to run this on all 6 of my samples, piping the output of bedtools to grep out only the lines starting with “all.”



Now that I have text files with coverage histograms for all the regions in the capture target, I can now plot this using R.



You can see that sample #2 had some problems. Relative to the rest of the samples you see it take a quick nose-dive on the left side of the plot relative to the others. Running this through Picard showed that 86% of the reads from this sample were duplicates.



Thanks, Aaron, for the tips.

Bedtools protocols: https://github.com/arq5x/bedtools-protocols