Home | Blog | Software | Publications | GitHub


Add labels to a genomic Hilbert curve under pixel mode

Making Hilbert curve for genomic data under “pixel” mode provides a high resolution way to visualize patterns both in a global and local scale. Under “pixel” mode, the curve is stored as an RGB matrix, and it is added to the graphic device as a raster image.

library(HilbertCurve)
library(circlize)
library(GenomicRanges)

set.seed(123)
bed = generateRandomBed(1000)
gr = GRanges(seqnames = bed[[1]], ranges = IRanges(bed[[2]], bed[[3]]), score = bed[[4]])
col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
hc = GenomicHilbertCurve(mode = "pixel", level = 10, title = "random bed")
hc_layer(hc, gr, col = col_fun(gr$score))

plot of chunk unnamed-chunk-2

Since there are multiple chromosomes, adding border for each chromosome helps to identify different chromosomes on the plot. Under “pixel” mode, hc_map() actually calculates the border of each chromosome and modify corresponding pixel in the RGB matrix to #808080. Then the raster image will be updated in the graphic device.

hc_map(hc, add = TRUE, fill = NA, border = "#808080")

plot of chunk unnamed-chunk-3

But still, it is not straightforward to tell which chromosome locates where. Also, it is impossible to add text directly to the plot because the Hilbert curve itself is stored as an RGB matrix. However, there is a workaround that we can add another Hilbert curve which only contains the labels for chromosomes on top of the first curve.

When making the plot, the curve itself belongs to a viewport with a name hilbert_curve_$i that the name can be obtained by paste0("hilbert_curve_", HilbertCurve:::.ENV$I_PLOT). Then we can go to that viewport by seekViewport() and add a second curve with the same setting as the first one expect the mode is set to “normal” and with a lower level (because we only want to locate each chromosome and a low level is sufficient for locating it, also lower level gives faster speed). Remember to set newpage = FALSE so that the second curve will not create a new graphic page.

seekViewport(paste0("hilbert_curve_", HilbertCurve:::.ENV$I_PLOT))

hc = GenomicHilbertCurve(mode = "normal", level = 6, newpage = FALSE)
hc_map(hc, add = TRUE, fill = NA, border = NA)

plot of chunk unnamed-chunk-4