Question: Gviz: ylim of the AlignmentsTrack
gravatar for wsciekly.maciek
8 months ago by
wsciekly.maciek0 wrote:

I'm having a problem with plotting the AlignmentsTrack in Gviz. My input is a bam file, I want to plot the coverage of a specified region. (~200nt). The problem is that I would like my histogram to be stretched to the ymax, therefore I want to set something like ylim=c(0,max_coverage_per_position).

Could you let me know how to solve this most efficiently?

gviz gviz-package • 284 views
ADD COMMENTlink modified 8 months ago by florian.hahne@novartis.com1.6k • written 8 months ago by wsciekly.maciek0
Answer: Gviz: ylim of the AlignmentsTrack
gravatar for
8 months ago by

I am not quite sure I understand. Unless provided explicitly, the ylim values should already be computed dynamically. If that is not the desired range, you will need to compute coverage at the region to be plotted from you bam file manually and explicitly set ylim with the max coverage value. There should be loads of instruction around explaining how to do that. E.g., Any package to calculate NGS coverage and plot it?

Bioconductor's Rle vectors support most arithmetic operation, including min and max.


ADD COMMENTlink written 8 months ago by florian.hahne@novartis.com1.6k

Dear Florian,
I did just as you said and something is inconsistent here...
I select the maximum coverage as:

        aln <- readGAlignments(opt$bam)
        cvg = coverage(aln)
        x = cvg[gr]

And it prints 352; gr holds the interesting genomic region for me.

Later I construct the Alignments track with ylim = c(0,max(x$"22"))). Please take a look at the attached end plot. The ylim is set correctly, but the coverage never seem to get close to 300...

Also, please notice a weird, isolated grey dot in the middle of the plot.
Could you please help me out how to get rid of it and why the max coverage on the plot looks different from the calculated one?

ADD REPLYlink written 8 months ago by wsciekly.maciek0

Gviz applies a bit of smoothing in order to make the coverage plot look a bit more appealing:

res <- .pxResolution(coord = "x")
brks <- ceiling((maxBase - mminBase)/res)
x <- seq(mminBase, maxBase, len = brks)
y <- tapply(as.integer(cov[mminBase:maxBase]), cut(mminBase:maxBase, breaks = brks), mean)

There's a bit of logic flaw in that it uses the actual coverage range for ylim, and not the smoothed range. Also I can see that having the ability to turn the smoothing off would be beneficial. Not currently do-able, but we can expose that with the next release.


ADD REPLYlink written 8 months ago by florian.hahne@novartis.com1.6k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 287 users visited in the last hour