Plot Reads Against Coordinates

Question

[Vikas Bansal]

Hi,

I have 2 files - bed file (chromosome, start and end) and second file is mapped reads which are in sam file. I want to plot reads against coordinates from my bed file. Example my bed file contains-

chr4   50    100
chr4   100   150
chr4   150   200
chr4   200   250

With this bed file, I can extract reads from sam file which are mapped between chr4: 50 and 250 positions. In the end, I want a plot like this (image taken from internet but I modified it little bit). Each read is shown according to its start position and length of the bar represents read length. I found this image and about python and matlab here. But I want to do this in R. The thing is, I have about 40 bed files and with the help of par(mfrow) in R, I will make 40 plots on same sheet. Any suggestions how would you do it or any R package will be great. I would really appreciate your help.

Thanks and Best regards,

Vikas

Answer 1

[Neilfws]

(1) There are very similar plots in the documentation (PDF) for the R/Bioconductor package ggbio. I've linked to v2.9 since documentation for the newer 2.10 is not yet written in full.

(2) I would not put 40 plots on the same "sheet".

[Vikas Bansal]

Thanks. You are right, it is not a good idea to put 40 plots on same sheet.

Answer 2

[Sean Davis]

This does not produce a "publication quality" set of figures, but you could use the SRAdb bioconductor package to control the IGV browser to automatically make snapshots of any regions of interest. There is an example in the SRAdb vignette or in the "Accomplishing Exome Sequencing Data Tasks with Bioconductor" tutorial on this page.