Question

Gviz: ylim of the AlignmentsTrack

1

Entering edit mode

wsciekly.maciek ▴ 10

@wscieklymaciek-17486

Last seen 5.4 years ago

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 • 2.2k views

ADD COMMENT • link updated 2.3 years ago by Pariksheet Nanda ▴ 120 • written 5.6 years ago by wsciekly.maciek ▴ 10

score 1 · Answer 1 · 2021-12-24

It's been a few years and looks like there still isn't any option to disable smoothing for AlignmentsTrack. At least I don't see any way to do it from the Gviz documentation and the the R source code in R/Gviz-methods.R in the drawGD method for AlignmentsTrack still forces the smoothing logic pointed out by Florian's comment when the read size is more than 2 pixels.

As workaround, I use DataTrack as explained in the section "Building DataTrack objects from files" section of the Gviz vignette copied here for convenience:

bamFile <- system.file("extdata/test.bam", package = "Gviz")
dTrack4 <- DataTrack(range = bamFile, genome = "hg19", type = "l",
                     name = "Coverage", window = -1,
                     chromosome = "chr1")
plotTracks(dTrack4, from = 189990000, to = 190000000)

... you don't get the benefit of the rug plot, fill color, etc. as you would with AlignmentsTrack, nevertheless using DataTrack seems practical for simple use cases of plotting coverage.

I have sparse reads, and use the type = "h" to get a histogram similar to the UCSC browser.

score 0 · Answer 2 · 2018-09-24

0

Entering edit mode

florian.hahne@novartis.com ★ 1.6k

@florianhahnenovartiscom-3784

Last seen 5.6 years ago

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 COMMENT • link 5.6 years ago florian.hahne@novartis.com ★ 1.6k

0

Entering edit mode

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 REPLY • link 5.6 years ago wsciekly.maciek ▴ 10

2

Entering edit mode

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 REPLY • link 5.6 years ago florian.hahne@novartis.com ★ 1.6k