Hi all

I’ve been searching for a way to plot bigwig coverage from my ATAC-seq data and illustrate differences in open chromatin at particular loci. I want to avoid having to screenshot tracks in IGV browser because its tedious and laborious.

I have recently come across the wiggleplotr package, which I really like but some more tweaks to this tool would make it even more widely usable tool, especially for programming noobs like myself.

The documentation is clear and I managed to generate all the necessary metatables to get it to work in my genome of choice. However, the package is very much geared towards RNA-seq data and in particular plotting coverage over alternative splice forms of a single gene.

I would like a broader solution for plotting bigwig coverage at a genetic locus as I am not interested in alternative splice forms per se. I would like a way to centre the plot at a certain gene but also include any flanking genes in the plot for accuracy - wiggleplotr only seems to deal with one gene at a time.

Has anyone managed to use wiggleplotr for other NGS datasets or can someone recommend a viable alternative? 

Thanks

Michael

__________________
Dr Michael Borg
Gregor Mendel Institute
Dr. Bohr Gasse 3
1030 Vienna
Austria
 

Dear Michael,

Thank you for the question. Yes, although it is not explicitly described in the vignette, wiggleplotr can be used to display ATAC-seq read coverage in a genomic region. For an example, see Figure 3D or Figure 4C in this paper :https://www.biorxiv.org/content/early/2017/05/18/102392. I will try to add an example of this to the next version of the vignette.

As you can see, Figure 3D has two different genes. Unfortunately there is no way to do this automatically at the moment, you have to manually provide plotCoverage with the transcripts that you want to plot (which could be from different genes). I would suggest picking one representative transcript per gene, because otherwise there will be too many rows (each gene/transcript always assigned to a new row, even if it does not overlap any other transcript).

The simplest way to do this would be with the plotCoverageFromEnsembldb wrapper function where you can list your specify your gene names and transcript ids. See the vignette example for plotTranscriptsFromEnsembldb. You can also change either the flanking_length or region_coords parameters to change the width of the genomic region around the transcript annotations. See the documentation page of plotCoverage for more details, it has many more options that are not covered in the vignette.

Hope this helps. Also, if you post more concrete examples then I might be able to help you with the code.

Best wishes,
Kaur