## ---- 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"), 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"), 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"), 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"), 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"), 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"), 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"), 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"), newpage = FALSE, start_from = 'topright', first_seg = "vertical", title = "start_from = 'topright'\nfirst_seg = 'vertical'") upViewport() upViewport() ## ------------------------------------------------------------------------------------------------- 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 = "#FF0000", 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) 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 = sample(1000, 100) value = runif(length(x)) hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode") hc_layer(hc, x1 = x, x2 = x+2, 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 = x, x2 = x+2, 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 = x, x2 = x+2, mean_mode = "absolute", col = col_fun(value), border = TRUE) ## ---- 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 = x, x2 = x+2, mean_mode = "absolute") hc_layer(hc, x1 = 750, x2 = 850, col = "#0000FF20", overlay = function(r0, g0, b0, r, g, b, alpha) { l = r0 > 0 # if the pixel is red # red -> green, just switch red and green channel tmp = r0[l] r0[l] = g0[l] g0[l] = tmp # add to the curve default_overlay(r0, g0, b0, r, g, b, alpha) }) ## ------------------------------------------------------------------------------------------------- hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode") hc_layer(hc, x1 = x, x2 = x+2, mean_mode = "absolute", col = col_fun(value)) hc_layer(hc, x1 = 750, x2 = 850, col = "#00000010", # remember r0, g0, ... only correspond to the pixels in the grey area overlay = function(r0, g0, b0, r, g, b, alpha) { # non-white areas l = !(r0 == 1 & g0 == 1 & b0 == 1) # original value v = col2value(r0[l], g0[l], b0[l], col_fun = col_fun) # new color theme 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] 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") cm = ColorMapping(col_fun = col_fun) lgd3 = color_mapping_legend(cm, title = "style 3", plot = FALSE) 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("window 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.lines(c(0.3, 0.5), c(0.6, 0.6), gp = gpar(lwd = 4)) grid.lines(c(0.7, 0.9), c(0.6, 0.6), gp = gpar(lwd = 4)) 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("genomic regions", 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 = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", echo = -c(21, 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("

\n") ## ---- 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, ...) ## ------------------------------------------------------------------------------------------------- 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 cm = ColorMapping(col_fun = colorRamp2(c(0, 1), c("yellow", "red"))) lgd = color_mapping_legend(cm, title = "Conservation", plot = FALSE, 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("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", echo = 2:11------ # 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")) # cm = ColorMapping(col_fun = col_fun) # lgd = color_mapping_legend(cm, title = "GC percent", plot = FALSE, # 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", fig.keep = "none", echo = c(2:4,6:7)---- # 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()) # # or you can only save the Hilbert curve part as a PNG file # # hc_png(hc, file = "gc_percent_chr1.png") ## ---- echo = FALSE, results = "asis"-------------------------------------------------------------- cat("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", fig.keep = "none", echo = c(2:6,8:9)---- # 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()) # # or you can only save the Hilbert curve part as a PNG file # # hc_png(hc, file = "gc_percent_half_chr1.png") ## ---- echo = FALSE, results = "asis"-------------------------------------------------------------- cat("

\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")) cm = ColorMapping(col_fun = col_fun) lgd = color_mapping_legend(cm, title = "Intensity", plot = FALSE) cm2 = ColorMapping(level = c(mark, "gene"), color = c("#FF0000", "#CCCCCC")) lgd2 = color_mapping_legend(cm2, title = "Layer", plot = FALSE) hc = GenomicHilbertCurve(chr = "chr1", level = 9, mode = "pixel", title = qq("Intensity of @{mark} mark on chr1"), legend = list(lgd, lgd2)) hc_layer(hc, df, col = col_fun(df[[5]])) hc_layer(hc, reduce(g), col = "#00000010") # or you can only save the Hilbert curve part as a PNG file # hc_png(hc, qq("@{mark}_chr1.png")) } ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", fig.keep = "none", echo = c(1:9, 11:25, 27:28)---- # 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")) # cm = ColorMapping(col_fun = col_fun) # lgd = color_mapping_legend(cm, title = "Intensity", plot = FALSE) # cm2 = ColorMapping(level = c("H3K36me3", "overlap", "gene"), color = c("#FF0000", "purple", "#CCCCCC")) # lgd2 = color_mapping_legend(cm2, title = "Layer", plot = FALSE) # # 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(lgd, 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 = !(r0 == 1 & g0 == 1 & b0 == 1) # 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()) # # or you can only save the Hilbert curve part as a PNG file # # hc_png(hc, "H3K36me3_chr1_overlap.png") ## ---- echo = FALSE, results = "asis"-------------------------------------------------------------- cat("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", fig.keep = "none", echo = c(1:5, 7:9, 11:12)---- # 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")) # cm = ColorMapping(col_fun = col_fun) # lgd = color_mapping_legend(cm, title = "Methylation", plot = FALSE) # 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") # hc_layer(hc, reduce(g), col = "#00000010") # invisible(dev.off()) # # or you can only save the Hilbert curve part as a PNG file # # hc_png(hc, file = 'methylation_chr1.png') ## ---- echo = FALSE, results = "asis"-------------------------------------------------------------- cat("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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 = !(r0 == 1 & g0 == 1 & b0 == 1) # r[l] = 160/255 # g[l] = 32/255 # b[l] = 240/255 # list(r, g, b) # }) # invisible(dev.off()) ## ---- echo = FALSE, results = "asis"-------------------------------------------------------------- cat("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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 = !(r0 == 1 & g0 == 1 & b0 == 1) # 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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 = !(r0 == 1 & g0 == 1 & b0 == 1) # 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("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", 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 = !(r0 == 1 & g0 == 1 & b0 == 1) # 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") # } # # 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("

\n") ## ------------------------------------------------------------------------------------------------- sessionInfo()