4
$\begingroup$

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]: enter image description here

[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   -
Improve this question
$\endgroup$
2
  • $\begingroup$ I am not familiar with this field, so I would like to know the terms I should avoid/are not good enough. What kind of keywords have you searched? $\endgroup$
    – llrs
    Jan 28, 2019 at 13:28
  • $\begingroup$ @llrs, I'm not exactly sure what you mean by "terms I should avoid/are not good enough." - but I believe coverage plots are one class of visualisation which are more similar to binned histograms in general. Except a google search would yield in methods related to plotting coverage of transcripts or reads acquired from sequencing experiments [keywords: coverage plots, bedfiles, metaplots]. The key difference here is that I'm interested in scaled coverage over whole gene body (TSS-TES) instead of positions relative to TSS - similar to deeptools plotProfile. $\endgroup$
    – Siddharth
    Jan 28, 2019 at 13:53

2 Answers 2

Reset to default
3
$\begingroup$

If one assumes that repeat elements of a given type (e.g., LINEs) don't overlap each other, then the following will work:

  1. Split your BED file by repeat element, such that you have a LINE.bed, SINE.bed, etc.
  2. Convert those to bedGraph (e.g., awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"1.0"}' LINE.bed > LINE.bedGraph).
  3. Use UCSC tools to convert those bedGraph files to BigWig.
  4. 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
  5. 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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I still have to try it, but this might just work. I will mark it as answer unless something goes wrong. I would use scale-regions instead of reference-point since I'm interested in the whole region length. $\endgroup$
    – Siddharth
    Jan 31, 2019 at 12:13
0
$\begingroup$

There is a reasonably nice looking tutorial on metaplots in R here: https://rpubs.com/achitsaz/94710

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for answering this, usually we'd like to have the answer self-contained. Perhaps you could improve the question by including the relevant packages/code from that website into your answer. $\endgroup$
    – llrs
    Jan 29, 2019 at 13:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.