Question: Gviz: ylim of the AlignmentsTrack
0
13 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 • 391 views
modified 13 months ago by florian.hahne@novartis.com1.6k • written 13 months ago by wsciekly.maciek0
Answer: Gviz: ylim of the AlignmentsTrack
0
13 months 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.

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?

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