Home | Blog | Software | Publications | GitHub
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[], ranges = IRanges(bed[], bed[]), score = bed[]) 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))
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
Then the raster image will be updated in the graphic device.
hc_map(hc, add = TRUE, fill = NA, border = "#808080")
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
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)