Hello,

I have a complex question for which I am trying to find a practical answer. I have some exome samples sequenced at various read depths (paired-end) using an exon target (capture) kit. I have pileup files for the samples and for each sample, I wish to see "what percentage of positions in the capture regions have what coverage".

So, for simplicity purposes, let's say the capture kit has 3 regions in BED format (chr, start, stop) :

chr1 1200 1300
chr5 2600 2700
chr7  3100 3200

I want to see the coverage of each base in each of these regions and put it into a category (let's say I have 10 categories - 1-10, 11-20...and so on till 91-100)

At the end I want to get a table for each region showing how many bases are in each category like :

chr    start    stop    1-10    11-20......... 91-100
chr1   1200   1300    25        57..............198
chr5   2600   2700   132       65............  123
chr7   3100   3200   189      156.........     203

Once I am done getting this table, I want to generate a "Percent Positions (Y axis) vs. Coverage (X-Axis)" plot where I plot each category of coverage against the percent of positions in the capture region that show this coverage.

If I am right, the formula to calculate percentage positions for each category would be :

(sum of number of bases in each region for that category) / (length of each region) * 100

So for 1-10 category, this would simply mean (25 + 132 + 189)/ (100 + 100 + 100) * 100

What are some of the fastest ways to do this analysis ? Does BEDTools have some tool which does this quickly ?