Search
Question: Gviz: ylim of the AlignmentsTrack
0
gravatar for wsciekly.maciek
11 weeks 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?

ADD COMMENTlink modified 11 weeks ago by florian.hahne@novartis.com1.6k • written 11 weeks ago by wsciekly.maciek0
0
gravatar for florian.hahne@novartis.com
11 weeks ago by
Switzerland

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 11 weeks 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]
        print(max(x$"22"))

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 11 weeks 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.

Florian

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

Help
Access

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