## ---- echo = FALSE, message = FALSE---------------------------------------------------------------
library(markdown)
options(markdown.HTML.options = c(options('markdown.HTML.options')[[1]], "toc"))

library(knitr)
knitr::opts_chunk$set(
    error = FALSE,
    tidy  = FALSE,
    message = FALSE,
    fig.align = "center",
    fig.width = 6,
    fig.height = 6)
options(markdown.HTML.stylesheet = "custom.css")

options(width = 100)

library(HilbertCurve)

## -------------------------------------------------------------------------------------------------
library(HilbertCurve)
library(circlize)
set.seed(12345)

## ---- fig.width = 12, fig.height = 3, echo = FALSE------------------------------------------------
grid.newpage()
pushViewport(viewport(layout = grid.layout(nr = 1, nc = 4)))
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 1))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, newpage = FALSE, title = "level = 2")
hc_segments(hc, x1 = 1, x2 = 2, gp = gpar(col = "red", lwd = 2))
upViewport()
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 2))
hc = HilbertCurve(1, 64, level = 3, reference = TRUE, newpage = FALSE, title = "level = 3")
hc_segments(hc, x1 = 1, x2 = 2, gp = gpar(col = "red", lwd = 2))
upViewport()
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 3))
hc = HilbertCurve(1, 256, level = 4, reference = TRUE, newpage = FALSE, title = "level = 4")
hc_segments(hc, x1 = 1, x2 = 2, gp = gpar(col = "red", lwd = 2))
upViewport()
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 4))
hc = HilbertCurve(1, 1024, level = 5, reference = TRUE, reference_gp = gpar(col = "grey"), 
    arrow = FALSE, newpage = FALSE, title = "level = 5")
hc_segments(hc, x1 = 1, x2 = 2, gp = gpar(col = "red", lwd = 2))
upViewport()
upViewport()

## ---- eval = FALSE--------------------------------------------------------------------------------
#  for(i in 1:1024) {
#      hc = HilbertCurve(1, 1024, level = 5, reference = TRUE, arrow = FALSE)
#      hc_points(hc, x1 = i, np = NULL, pch = 16, size = unit(2, "mm"))
#  }

## ---- fig.width = 9, fig.height = 8---------------------------------------------------------------
library(HilbertVis)
pos = HilbertVis::hilbertCurve(5)
mat = as.matrix(dist(pos))
library(ComplexHeatmap)

ht = Heatmap(mat, name = "dist", cluster_rows = FALSE, cluster_columns = FALSE, 
    show_row_names = FALSE, show_column_names = FALSE, 
    heatmap_legend_param = list(title = "euclidean_dist"))
draw(ht, padding = unit(c(5, 5, 5, 2), "mm"))
decorate_heatmap_body("dist", {
    grid.segments(c(0.25, 0.5, 0.75, 0, 0, 0), c(0, 0, 0, 0.25, 0.5, 0.75), 
          c(0.25, 0.5, 0.75, 1, 1, 1), c(1, 1, 1, 0.25, 0.5, 0.75), gp = gpar(lty = 2))
    grid.text(rev(c(256, 512, 768, 1024)), 0, c(0, 256, 512, 768)/1024, just = "bottom", 
        rot = 90, gp = gpar(fontsize = 10))
    grid.text(c(1, 256, 512, 768, 1024), c(1, 256, 512, 768, 1024)/1024, 1, just = "bottom",
        gp = gpar(fontsize = 10))
})

## ---- eval = FALSE--------------------------------------------------------------------------------
#  hc = HilbertCurve(...)    # initialize the curve
#  hc_points(hc, ...)        # add points
#  hc_rect(hc, ...)          # add rectangles
#  hc_polygon(hc, ...)       # add polygons
#  hc_segments(hc, ...)      # add lines
#  hc_text(hc, ...)          # add text

## ---- eval = FALSE--------------------------------------------------------------------------------
#  hc = HilbertCurve(1, 100, level = 4, reference = TRUE)

## ---- fig.width = 12, fig.height = 6, echo = FALSE------------------------------------------------
grid.newpage()
pushViewport(viewport(layout = grid.layout(nr = 2, nc = 4)))
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 1))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'bottomleft', first_seg = "horizontal", 
    title = "start_from = 'bottomleft'\nfirst_seg = 'horizontal'")
upViewport()
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 2))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'topleft', first_seg = "horizontal", 
    title = "start_from = 'topleft'\nfirst_seg = 'horizontal'")
upViewport()
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 3))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'bottomright', first_seg = "horizontal",
    title = "start_from = 'bottomright'\nfirst_seg = 'horizontal'")
upViewport()
pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 4))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'topright', first_seg = "horizontal",
    title = "start_from = 'topright'\nfirst_seg = 'horizontal'")
upViewport()

pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 1))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'bottomleft', first_seg = "vertical", 
    title = "start_from = 'bottomleft'\nfirst_seg = 'vertical'")
upViewport()
pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 2))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'topleft',  first_seg = "vertical",
    title = "start_from = 'topleft'\nfirst_seg = 'vertical'")
upViewport()
pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 3))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'bottomright',  first_seg = "vertical",
    title = "start_from = 'bottomright'\nfirst_seg = 'vertical'")
upViewport()
pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 4))
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, reference_gp = gpar(col = "#999999", lwd = 1), 
    newpage = FALSE, start_from = 'topright',  first_seg = "vertical",
    title = "start_from = 'topright'\nfirst_seg = 'vertical'")
upViewport()

upViewport()

## -------------------------------------------------------------------------------------------------
library(IRanges)
x = sort(sample(100, 20))
s = x[1:10*2 - 1]
e = x[1:10*2]
ir = IRanges(s, e)
ir

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
hc_points(hc, ir)

## ---- fig.width = 3, fig.height = 3---------------------------------------------------------------
hc = HilbertCurve(1, 16, level = 2, reference = TRUE, title = "np = 3")
hc_points(hc, x1 = 1, x2 = 2, np = 3)

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
hc_points(hc, ir, np = 3, gp = gpar(fill = rand_color(length(ir))),
    shape = sample(c("circle", "square", "triangle", "hexagon", "star"), length(ir), replace = TRUE))

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
hc_points(hc, ir, np = NULL, size = unit(runif(length(ir)), "cm"), pch = 16)

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
hc_segments(hc, ir, gp = gpar(lwd = 5))

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
hc_rect(hc, ir, gp = gpar(fill = "#FF000080"))

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
hc_polygon(hc, ir, gp = gpar(fill = "red", col = "black"))

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
labels = sample(letters, length(ir), replace = TRUE)
hc_text(hc, ir, labels = labels, gp = gpar(fontsize = width(ir)*2+5))

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4, reference = TRUE)
hc_polygon(hc, ir)
hc_text(hc, ir, labels = "a")
hc_text(hc, ir, labels = "b", centered_by = "polygon")

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 4)
hc_segments(hc, IRanges(1, 100))   # This is an other way to add background line
hc_rect(hc, ir, gp = gpar(fill = rand_color(length(ir), transparency = 0.8)))
hc_polygon(hc, ir[c(1,3,5)], gp = gpar(col = "red"))
hc_points(hc, ir, np = 3, gp = gpar(fill = rand_color(length(ir))),
    shape = sample(c("circle", "square", "triangle", "hexagon", "star"), length(ir), replace = TRUE))
hc_text(hc, ir, labels = labels, gp = gpar(fontsize = width(ir)*2+5, col = "blue", font = 2))

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(0.1, 0.8, level = 4, reference = TRUE)
hc_points(hc, x1 = c(0.15, 0.55), x2 = c(0.25, 0.65))

## ---- eval = FALSE--------------------------------------------------------------------------------
#  # code not run when generating the vignette
#  hc = HilbertCurve(1, 1000000000000, level = 4, reference = TRUE)
#  hc_points(hc, x1 = 400000000000, x2 = 600000000000)

## ---- eval = FALSE--------------------------------------------------------------------------------
#  # code not run when generating the vignette
#  hc = HilbertCurve(-100, 100, level = 4, reference = TRUE)
#  hc_points(hc, x1 = -50, x2 = 50)

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 9, mode = "pixel")
hc_layer(hc, ir)

## -------------------------------------------------------------------------------------------------
col_fun = colorRamp2(c(0, 1), c("white", "red"))
x = seq(10, 960, length = 100)
x1 = x# start of each interval
x2 = x + 2    # end of each interval
value = runif(100)  # associated values
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode")
hc_layer(hc, x1 = x1, x2 = x2, mean_mode = "absolute", col = col_fun(value))
hc_layer(hc, x1 = 750, x2 = 850, col = "#00000010")

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode")
hc_layer(hc, x1 = x1, x2 = x2, mean_mode = "absolute", col = col_fun(value), grid_line = 3)

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode")
hc_layer(hc, x1 = x1, x2 = x2, mean_mode = "absolute", col = col_fun(value), border = TRUE)

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode")
hc_layer(hc, x1 = x1, x2 = x2, mean_mode = "absolute", col = col_fun(value))
hc_polygon(hc, x1 = x1, x2 = x2)

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode")
hc_layer(hc, x1 = x1, x2 = x2, mean_mode = "absolute", col = col_fun(value))
hc_polygon(hc, x1 = c(1, 200, 500), x2 = c(200, 500, 1000))
hc_text(hc, x1 = c(1, 200, 500), x2 = c(200, 500, 1000), 
    labels = c("A", "B", "C"), gp = gpar(fontsize = 20), centered_by = "polygon")

## ---- eval = FALSE--------------------------------------------------------------------------------
#  # code not run, only for demonstration
#  hc_png(hc, file = "test.png")

## ---- eval = FALSE--------------------------------------------------------------------------------
#  r = r*alpha + r0*(1-alpha)
#  g = g*alpha + g0*(1-alpha)
#  b = b*alpha + b0*(1-alpha)

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode")
hc_layer(hc, x1 = x1, x2 = x2, mean_mode = "absolute")
hc_layer(hc, x1 = c(300, 750), x2 = c(400, 850), col = "#0000FF20",
    overlay = function(r0, g0, b0, r, g, b, alpha) { # it's only applied in [300, 400] and [750, 850]
        l = is_white(r0, g0, b0) # `is_white` simple tests whether it is the white color
        
        # change the red background to green
        if(any(!l)) {
            r0[!l] = 0 
            g0[!l] = 1 
            b0[!l] = 0 
        }

        # overlay #0000FF20 to the background
        default_overlay(r0, g0, b0, r, g, b, alpha)
    })

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode")
hc_layer(hc, x1 = x1, x2 = x2, mean_mode = "absolute", col = col_fun(value))
hc_layer(hc, x1 = c(300, 750), x2 = c(400, 850), col = "#00000010",
    overlay = function(r0, g0, b0, r, g, b, alpha) {
        # non-white pixels
        l = !is_white(r0, g0, b0)

        # original value
        v = col2value(r0[l], g0[l], b0[l], col_fun = col_fun)

        # new color schema
        col_fun_new = colorRamp2(c(0, 1), c("white", "blue"))
        col_new = col_fun_new(v, return_rgb = TRUE)
        r0[l] = col_new[, 1]
        g0[l] = col_new[, 2]
        b0[l] = col_new[, 3]

        # overlay #00000010 to the background
        default_overlay(r0, g0, b0, r, g, b, alpha)
    })

## -------------------------------------------------------------------------------------------------
value = runif(length(ir))
col_fun = colorRamp2(c(0, 1), c("white", "red"))
legend1 = Legend(at = seq(0, 1, by = 0.2), col_fun = col_fun, title = "continuous")
legend2 = Legend(at = c("A", "B"), legend_gp = gpar(fill = c("#00FF0080", "#0000FF80")), 
    title = "discrete")

legend = list(legend1, legend2)

hc = HilbertCurve(1, 100, reference = TRUE, title = "points", legend = legend)
hc_points(hc, ir, np = 3, gp = gpar(fill = col_fun(value)))
hc_rect(hc, ir, gp = gpar(fill = sample(c("#00FF0020", "#0000FF20"), length(ir), replace = TRUE)))

## -------------------------------------------------------------------------------------------------
col_fun = colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
breaks = seq(-2, 2, by = 2)
lgd1 = Legend(at = breaks, col_fun = col_fun, title = "style 1")
lgd2 = Legend(at = rev(breaks), legend_gp = gpar(fill = col_fun(rev(breaks))), title = "style 2")
lgd3 = Legend(col_fun = col_fun, title = "style 3")

value = rnorm(length(x))
hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode", 
    legend = list(lgd1, lgd2, lgd3))
hc_layer(hc, x1 = x, x2 = x+2, mean_mode = "absolute", col = col_fun(value))

## ---- fig.height = 3, echo = FALSE----------------------------------------------------------------
library(grid)
grid.lines(c(0.1, 0.9), c(0.3, 0.3), gp = gpar(col = "red", lwd = 4))
grid.text("segment on the curve", 0.5, 0.25, just = "top", gp = gpar(col = "red"))
grid.lines(c(0, 0.2), c(0.6, 0.6))
grid.lines(c(0.3, 0.5), c(0.6, 0.6))
grid.lines(c(0.7, 1), c(0.6, 0.6))
grid.lines(c(0.1, 0.2), c(0.6, 0.6), gp = gpar(lwd = 4)); grid.text("x_1", 0.15, 0.58, just = "top")
grid.lines(c(0.3, 0.5), c(0.6, 0.6), gp = gpar(lwd = 4)); grid.text("x_2", 0.4, 0.58, just = "top")
grid.lines(c(0.7, 0.9), c(0.6, 0.6), gp = gpar(lwd = 4)); grid.text("x_3", 0.8, 0.58, just = "top")
grid.lines(c(0.1, 0.1), c(0, 1), gp = gpar(lty = 2, col = "grey"))
grid.lines(c(0.9, 0.9), c(0, 1), gp = gpar(lty = 2, col = "grey"))
grid.text("input intervals", 0.5, 0.65, just = "bottom")

## -------------------------------------------------------------------------------------------------
col = rainbow(100)
hc = HilbertCurve(1, 100, level = 5)
hc_points(hc, x1 = 1:99, x2 = 2:100, np = 3, gp = gpar(col = col, fill = col))

## -------------------------------------------------------------------------------------------------
hc = HilbertCurve(1, 100, level = 5)
hc_rect(hc, x1 = 1:99, x2 = 2:100, gp = gpar(col = col, fill = col))

## ---- fig.width = 7, fig.height = 6---------------------------------------------------------------
library(ComplexHeatmap)
library(RColorBrewer)

load(system.file("extdata", "cidr_list.RData", package = "HilbertCurve"))

cidr_list = cidr_list[sapply(cidr_list, length) > 0]

country = rep(names(cidr_list), times = sapply(cidr_list, length))

ip = unlist(cidr_list, "r")

# convert ip address to numbers
mat = t(as.matrix(data.frame(lapply(strsplit(ip, "\\.|/"), as.numeric))))
start = mat[, 1]*256^3 + mat[, 2]*256^2 + mat[, 3]*256 + mat[, 4]
width = sapply(mat[, 5], function(x) strtoi(paste(rep(1, 32 - x), collapse = ""), base = 2))

# top 8 countries
col = structure(rep("grey", length(cidr_list)), names = names(cidr_list))
top8_rate = sort(tapply(width, country, sum), decreasing = TRUE)[1:8]/256^4
top8 = names(top8_rate)
col[top8] = brewer.pal(8, "Set1")
top8_rate = paste0(round(top8_rate*100), "%")

# this is the part of using HilbertCurve package
lgd = Legend(at = c(top8, "Others"), labels = c(paste(top8, top8_rate), "Others"),
    legend_gp = gpar(fill = c(col[top8], "grey")), title = "Top8 contries")
hc = HilbertCurve(0, 256^4, start_from = "topleft", first_seg = "horizontal",
    mode = "pixel", level = 9, legend = lgd)
hc_layer(hc, x1 = start, x2 = start + width, col = col[country])

## ---- eval = FALSE, echo = -c(20, 30)-------------------------------------------------------------
#  load(system.file("extdata", "chinese_dynasty.RData", package = "HilbertCurve"))
#  
#  detect_os = function() {
#      if (grepl('w|W', .Platform$OS.type)) {
#          os = "Windows"
#      } else {
#          if (grepl('darwin', version$os)) {
#              os = "MacOS"
#          } else {
#              os = "Linux"
#          }
#      }
#      return(os)
#  }
#  # default font family for Chinese under different OS
#  fontfamily = switch(detect_os(),
#      Windows = "SimSun",
#      MacOS = "Songti SC",
#      Linux = "Songti SC")
#  
#  png("chinese_dynasty.png", width = 500, height = 500)
#  hc = HilbertCurve(min(chinese_dynasty[[2]]), max(chinese_dynasty[[3]]),
#      title = "Chinese dynasty (1122 B.C. ~ 1911 A.D.)",
#      title_gp = gpar(fontsize = 16, fontfamily = fontfamily))
#  hc_segments(hc, x1 = chinese_dynasty[[2]], x2 = chinese_dynasty[[3]],
#      gp = gpar(col = rand_color(nrow(chinese_dynasty), transparency = 0.2),
#          lwd = runif(nrow(chinese_dynasty), min = 1, max = 10)))
#  hc_text(hc, x1 = chinese_dynasty[[2]], x2 = chinese_dynasty[[3]], labels = chinese_dynasty[[1]],
#      gp = gpar(fontsize = (chinese_dynasty[[3]] - chinese_dynasty[[2]])/500 * 10 + 8,
#          fontfamily = fontfamily))
#  
#  year = seq(-1000, 2000, by = 100)
#  hc_text(hc, x1 = year, labels = ifelse(year < 0, paste0(-year, "BC"), year),
#      gp = gpar(fontsize = 8))
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='chinese_dynasty.png' /></center></p>\n")

## ---- eval = FALSE--------------------------------------------------------------------------------
#  hc = GenomicHilbertCurve(chr, species, ...)
#  hc = GenomicHilbertCurve(background, ...)

## ---- eval = FALSE--------------------------------------------------------------------------------
#  hc_points(hc, gr, ...)
#  hc_rect(hc, gr, ...)
#  hc_polygon(hc, gr, ...)
#  hc_segments(hc, gr, ...)
#  hc_text(hc, gr, ...)
#  hc_layer(hc, gr, ...)

## -------------------------------------------------------------------------------------------------
library(GenomicRanges)
load(system.file("extdata", "refseq_chr1.RData", package = "HilbertCurve"))
hc = GenomicHilbertCurve(chr = "chr1", level = 5, reference = TRUE, 
    reference_gp = gpar(lty = 1, col = "grey"), arrow = FALSE)
hc_segments(hc, g, gp = gpar(lwd = 6, col = rand_color(length(g))))

## -------------------------------------------------------------------------------------------------
# for generating the legend
lgd = Legend(col_fun = colorRamp2(c(0, 1), c("yellow", "red")), 
    title = "Conservation",
    at = c(0, 0.2, 0.4, 0.6, 0.8, 1), 
    labels = c("0%", "20%", "40%", "60%", "80%", "100%"))
chr1_len = 249250621

## ---- fig.width = 7.5, fig.height = 7-------------------------------------------------------------
load(system.file("extdata", "mouse_net.RData", package = "HilbertCurve"))
seqlengths(mouse) = chr1_len # it is only used to extract the complement
nonmouse = gaps(mouse); nonmouse = nonmouse[strand(nonmouse) == "*"]
gr = c(mouse, nonmouse)
col = c(rep("red", length(mouse)), rep("yellow", length(nonmouse)))
hc = GenomicHilbertCurve(chr = "chr1", level = 6, 
    title = "Conservation between mouse and human on chr1",
    legend = lgd)
hc_points(hc, gr, np = 3, gp = gpar(col = NA, fill = col))

## ---- fig.width = 7.5, fig.height = 7-------------------------------------------------------------
load(system.file("extdata", "zebrafish_net.RData", package = "HilbertCurve"))
seqlengths(zebrafish) = chr1_len
nonzebrafish = gaps(zebrafish); nonzebrafish = nonzebrafish[strand(nonzebrafish) == "*"]
gr = c(zebrafish, nonzebrafish)
col = c(rep("red", length(zebrafish)), rep("yellow", length(nonzebrafish)))
hc = GenomicHilbertCurve(chr = "chr1", level = 6, 
    title = "Conservation between zebrafish and human on chr1",
    legend = lgd)
hc_points(hc, gr, np = 3, gp = gpar(col = NA, fill = col))

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, echo = 2:10---------------------------
#  png('gc_percent_chr1_points.png', width = 500, height = 500)
#  df = read.table(pipe("awk '$1==\"chr1\"' ~/HilbertCurveTest/hg19_gc_percent_window1000.bed"))
#  col_fun = colorRamp2(quantile(df[[5]], c(0.1, 0.5, 0.9)), c("green", "#FFFFCC", "red"))
#  lgd = Legend(col_fun = col_fun, title = "GC percent",
#      at = c(300, 400, 500, 600),
#      labels = c("30%", "40%", "50%", "60%"))
#  
#  hc = GenomicHilbertCurve(chr = "chr1", level = 6, legend = lgd)
#  hc_points(hc, df, np = 3, gp = gpar(fill = col_fun(df[[5]]), col = col_fun(df[[5]])))
#  hc_rect(hc, reduce(g), gp = gpar(fill = "#00000020", col = NA))
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='gc_percent_chr1_points.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = 2:4---------
#  png("gc_percent_chr1.png", width = 500, height = 500)
#  hc = GenomicHilbertCurve(chr = "chr1", level = 9, mode = "pixel", legend = lgd)
#  hc_layer(hc, df, col = col_fun(df[[5]]))
#  hc_layer(hc, reduce(g), col = "#00000020")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='gc_percent_chr1.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = 2:6---------
#  png("gc_percent_half_chr1.png", width = 500, height = 500)
#  background = GRanges(seqnames = "chr1", ranges = IRanges(1, ceiling(chr1_len/2)))
#  hc = GenomicHilbertCurve(background = background, level = 9, mode = "pixel", legend = lgd,
#      title = "First half of chromosome 1")
#  hc_layer(hc, df, col = col_fun(df[[5]]))
#  hc_layer(hc, reduce(g), col = "#00000020")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='gc_percent_half_chr1.png' /></center></p>\n")

## -------------------------------------------------------------------------------------------------
library(GetoptLong)
plot_histone_mark = function(mark) {
    df = read.table(pipe(qq("awk '$5>0 && $1==\"chr1\"' ~/HilbertCurveTest/UCSD.Lung.@{mark}.STL002.bed")), 
        sep = "\t")
    col_fun = colorRamp2(c(0, quantile(df[[5]], 0.99)), c("white", "red"))
    lgd1 = Legend(col_fun = col_fun, title = "Intensity")
    lgd2 = Legend(labels = c(mark, "gene"), legend_gp = gpar(fill = c("#FF0000", "#CCCCCC")),
        title = "Layer")

    hc = GenomicHilbertCurve(chr = "chr1", level = 9, mode = "pixel", 
        title = qq("Intensity of @{mark} mark on chr1"), legend = list(lgd1, lgd2))
    hc_layer(hc, df, col = col_fun(df[[5]]))
    hc_layer(hc, reduce(g), col = "#00000010")
}

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = 2-----------
#  png("H3K27ac_chr1.png", width = 500, height = 500)
#  plot_histone_mark("H3K27ac")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='H3K27ac_chr1.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = 2-----------
#  png("H3K36me3_chr1.png", width = 500, height = 500)
#  plot_histone_mark("H3K36me3")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='H3K36me3_chr1.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = 2-----------
#  png("H3K4me3_chr1.png", width = 500, height = 500)
#  plot_histone_mark("H3K4me3")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='H3K4me3_chr1.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = 2-----------
#  png("H3K9me3_chr1.png", width = 500, height = 500)
#  plot_histone_mark("H3K9me3")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='H3K9me3_chr1.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = c(1:9, 11:25)----
#  df = read.table(pipe("awk '$5>0 && $1==\"chr1\"' ~/HilbertCurveTest/UCSD.Lung.H3K36me3.STL002.bed"),
#      sep = "\t")
#  col_fun = colorRamp2(c(0, quantile(df[[5]], 0.99)), c("white", "red"))
#  col_fun_new = colorRamp2(c(0, quantile(df[[5]], 0.99)), c("white", "purple"))
#  lgd1 = Legend(col_fun = col_fun, title = "Intensity")
#  lgd2 = Legend(labels = c("H3K36me3", "overlap", "gene"),
#      legend_gp = gpar(fill = c("#FF0000", "purple", "#CCCCCC")), title = "Layer")
#  
#  png("H3K36me3_chr1_overlap.png", width = 500, height = 500)
#  hc = GenomicHilbertCurve(chr = "chr1", level = 9, mode = "pixel",
#      title = "Intensity of H3K36me3 mark on chr1", legend = list(lgd1, lgd2))
#  hc_layer(hc, df, col = col_fun(df[[5]]))
#  hc_layer(hc, reduce(g), col = "#00000010",
#      overlay = function(r0, g0, b0, r, g, b, alpha) {
#          l = !is_white(r0, g0, b0)
#          v = col2value(r0[l], g0[l], b0[l], col_fun = col_fun)
#  
#          col_new = col_fun_new(v, return_rgb = TRUE)
#          r0[l] = col_new[, 1]
#          g0[l] = col_new[, 2]
#          b0[l] = col_new[, 3]
#  
#          default_overlay(r0, g0, b0, r, g, b, alpha)
#      })
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='H3K36me3_chr1_overlap.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = c(1:4, 6:7)----
#  df = read.table(pipe("awk '$1==\"chr1\"' ~/HilbertCurveTest/UCSD.Lung.Bisulfite-Seq.STL002.bed"),
#      sep = "\t")
#  col_fun = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))
#  lgd = Legend(col_fun = col_fun, title = "Methylation")
#  png("methylation_chr1.png", width = 500, height = 500)
#  hc = GenomicHilbertCurve(chr = "chr1", level = 9, mode = "pixel", legend = lgd)
#  hc_layer(hc, df, col = col_fun(df[[5]]), mean_mode = "absolute")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='methylation_chr1.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = -c(1, 16)----
#  png("cnv_all_chromosomes.png", width = 500, height = 500)
#  hc = GenomicHilbertCurve(chr = paste0("chr", 1:22), level = 9, mode = "pixel")
#  df_gain = read.table("~/HilbertCurveTest/Stringent.Gain.hg19.2015-02-03.txt",
#      header = TRUE, stringsAsFactors = FALSE)
#  hc_layer(hc, df_gain, col = "red")
#  df_loss = read.table("~/HilbertCurveTest/Stringent.Loss.hg19.2015-02-03.txt",
#      header = TRUE, stringsAsFactors = FALSE)
#  hc_layer(hc, df_loss, col = "green", grid_line = 3, grid_line_col = "grey",
#      overlay = function(r0, g0, b0, r, g, b, alpha) {
#          l = !is_white(r0, g0, b0)
#          r[l] = 160/255
#          g[l] = 32/255
#          b[l] = 240/255
#          list(r, g, b)
#  })
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='cnv_all_chromosomes.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, echo = 2, fig.keep = "none"-----------
#  png('map_all_chromosomes.png', width = 500, height = 500)
#  hc_map(hc)
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='map_all_chromosomes.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = -c(1, 13)----
#  png('cnv_all_chromosomes_overlay.png', width = 500, height = 500)
#  hc = GenomicHilbertCurve(chr = paste0("chr", 1:22), level = 9, mode = "pixel")
#  hc_layer(hc, df_gain, col = "red")
#  hc_layer(hc, df_loss, col = "green",
#      overlay = function(r0, g0, b0, r, g, b, alpha) {
#          l = !is_white(r0, g0, b0)
#          r[l] = 160/255
#          g[l] = 32/255
#          b[l] = 240/255
#          list(r, g, b)
#  })
#  hc_map(hc, add = TRUE, fill = rand_color(22, transparency = 0.8))
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='cnv_all_chromosomes_overlay.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = -c(1, 13)----
#  png('cnv_all_chromosomes_border.png', width = 500, height = 500)
#  hc = GenomicHilbertCurve(chr = paste0("chr", 1:22), level = 9, mode = "pixel")
#  hc_layer(hc, df_gain, col = "red")
#  hc_layer(hc, df_loss, col = "green",
#      overlay = function(r0, g0, b0, r, g, b, alpha) {
#          l = !is_white(r0, g0, b0)
#          r[l] = 160/255
#          g[l] = 32/255
#          b[l] = 240/255
#          list(r, g, b)
#  })
#  hc_map(hc, add = TRUE, fill = NA, border = "grey")
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='cnv_all_chromosomes_border.png' /></center></p>\n")

## ---- eval = grepl("odcf", Sys.info()["nodename"]) & FALSE, fig.keep = "none", echo = 2:32--------
#  png('cnv_compare.png', width = 900, height = 300)
#  plot_curve = function(column, ...) {
#      cm = ColorMapping(levels = c("gain", "loss", "both"), color = c("red", "green", "purple"))
#      hc = GenomicHilbertCurve(chr = paste0("chr", 1:22), level = 9, mode = "pixel",
#          title = qq("CNV for @{column}"), ...)  ## `qq()` is from the GetoptLong package
#  
#      df = read.table("~/HilbertCurveTest/Stringent.Gain.hg19.2015-02-03.txt", header = TRUE,
#          stringsAsFactors = FALSE)
#      df = df[df[[column]] > 0, , drop = FALSE]
#      hc_layer(hc, df, col = "red")
#  
#      df = read.table("~/HilbertCurveTest/Stringent.Loss.hg19.2015-02-03.txt", header = TRUE,
#          stringsAsFactors = FALSE)
#      df = df[df[[column]] > 0, , drop = FALSE]
#      hc_layer(hc, df, col = "green",
#          overlay = function(r0, g0, b0, r, g, b, alpha) {
#              l = !is_white(r0, g0, b0)
#              r[l] = 160/255
#              g[l] = 32/255
#              b[l] = 240/255
#              list(r, g, b)
#      })
#      hc_map(hc, labels = 1:22, add = TRUE, fill = NA, border = "grey")
#  }
#  
#  pn = c("African", "Asian", "European")
#  pushViewport(viewport(layout = grid.layout(nrow = 1, ncol = 3)))
#  for(i in seq_along(pn)) {
#      pushViewport(viewport(layout.pos.row = 1, layout.pos.col = i))
#      plot_curve(pn[i], newpage = FALSE)
#      upViewport()
#  }
#  upViewport()
#  
#  invisible(dev.off())

## ---- echo = FALSE, results = "asis"--------------------------------------------------------------
cat("<p><center><img src='cnv_compare.png' /></center></p>\n")

## -------------------------------------------------------------------------------------------------
sessionInfo()