Plotting coverage of annotation over collection of region I'm trying to plot "meta" coverage of annotation: i.e. features (eg. gene class) over certain regions. It is similar to read coverage plots over gene body, except my input is two bed files (both in BED6 format) - (A) one containing the regions for which the meta plot is needed, and another being the list of genes or any feature whose coverage over the regions in (A) is needed. 

Is there any package or tool which can create such plots (my domain is limited to python but I can try to work with R) ?

Something akin to this but for whole gene body[1]:


[1]- Kelley, D, Rinn, J (2012). Transposable elements reveal a stem cell-specific class of long noncoding RNAs. Genome Biol., 13, 11:R107.

Edit: The BED file containing list of regions looks like this:

chr3    39218734    40053659    region1 0.92426187419769    +
chr4    140163762   140453127   region2 0.896103896103896   -
chr7    40549151    41205036    region3 0.986072423398329   +
chr8    81291743    81963246    region4 0.94184168012924    -
chr9    12284032    12539789    region5 0.95539033457249    -


And the bed file containing features to be plotted looks like this:

chr3    39218100    40053200    LINE    1   +
chr4    140163962   140453027   LINE 1  -
chr7    40549002    41204999    SINE 1  +
chr8    81291143    81963846    LTR 1   -
chr9    12284332    12539720    LTR 1   -
If one assumes that repeat elements of a given type (e.g., LINEs) don't overlap each other, then the following will work:


Split your BED file by repeat element, such that you have a LINE.bed, SINE.bed, etc.
Convert those to bedGraph (e.g., awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"1.0"}' LINE.bed &gt; LINE.bedGraph).
Use UCSC tools to convert those bedGraph files to BigWig.
Install deepTools and run computeMatrix reference-point -b 2500 -a 20 -S LINE.bigWig SINE.bigWig LTR.bigWig -R Regions_Of_Interest.BED -o foo.mat.gz
Make a profile plot with plotProfile (plotProfile -m foo.mat.gz -o foo.png --perGroup)


If you do have overlapping regions you'll need to first make a disjoint set of intervals (you can use bedops for this), find the coverage of them (e.g., bedtools coverage) and then continue on with that.