IanSudbery
commented
10 years ago I also have this problem.
I have a BigWig file I would like to visualize, which shows the mapped read depth from an RNA-seq experiment. This file is one of the ENCODE rnaseq files. Importantly the bigwig only contains entries where there is expression. Bases with zero mapped depth are excluded.
I load a part of this up into R using rtracklayer's import command:
>gr = GRanges(seqnames=c("chr1"), IRanges(from=159500770, to=160000770))
>expres <- import("bigwig.bw", which = gr, asRangedData = F)
>expres
GRanges with 218180 ranges and 1 metadata column
Note that the window is 1Mb big, but my GRanges object only has 218180 entries
I now want to plot this, but if I do it straight off, then the resulting computation takes forever, and would eventaully return more data than can be visualized. So aggregate over windows:
>expression_plot = ggplot(expres, aes(y=score)) + stat_aggregate(aes(y=score),window=1000, geom="bar") + scale_y_log10() + ylab("Read\nDepth") + xlim(gr)
This plots the expected coverage plot at a suitable resolution.
The problem comes when I try to combine this with another track:
> gene_plot <- ggplot(genes) + geom_alignment(aes(group=gene_id, fill = strand, color=strand),range.geom="arrowrect", group.selfish = F, extend.size =50000) + xlim(gr)
> tracks(Genes=gene_plot,, Expression=expression_plot, heights=c(3,1)) + theme_alignment(border = F, grid = F) + xlim(gr)
Error in if (zero_range(range)) { : missing value where TRUE/FALSE needed
>traceback()
5: scales::expand_range(getLimits(grob)$ylim, mul = 0.05)
4: FUN(X[[2L]], ...)
3: lapply(dots[!fixs & ideo], function(grob) {
scales::expand_range(getLimits(grob)$ylim, mul = 0.05)
})
2: lapply(dots[!fixs & ideo], function(grob) {
scales::expand_range(getLimits(grob)$ylim, mul = 0.05)
})
1: tracks(Genes = gene_plot, Expression = exp_plot, heights = c(3,
1))I've managed to trace the problem to the fact that stat_aggregate returns NA values if there are no features in a particular window. This NA then propagates through to the data of the layer, and so getLimits returns NA, which is a problem for expand_ranges, although plotting the grob on its own seems fine.
I note the patch above, but my R development skills arn't good enough to apply it myself to our system.
tengfei
Hi Michael,
I found a bug in the bbgio function getLimits: NA values present in the y values give NA for both min and max values. This leads to an error when combining plots using tracks. Here's the line in particular:
res <- suppressWarnings(list(xlim = c(min(c(l.res$xmin, x, xmin)), max(c(l.res$xmax, x, xmax, xend))), ylim = c(min(c(l.res$ymin, y, ymin)), max(c(l.res$ymax, y, ymax, yend)))))
here's my fix:
res <- suppressWarnings(list(xlim = c(min(c(l.res$xmin, x, xmin), na.rm=TRUE), max(c(l.res$xmax, x, xmax, xend), na.rm=TRUE)), ylim = c(min(c(l.res$ymin, y, ymin), na.rm=TRUE), max(c(l.res$ymax, y, ymax, yend), na.rm=TRUE))))
There's also some unusual behavior when combining tracks when you set legend.position="none". I've attached an example plot. There's excess space on the right margin where the legend should be but isn't.
J