## ---- 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")) cm1 = ColorMapping(col_fun = col_fun) legend1 = color_mapping_legend(cm1, plot = FALSE, title = "continuous") cm2 = ColorMapping(levels = c("A", "B"), color = c("#00FF0080", "#0000FF80")) legend2 = color_mapping_legend(cm2, plot = FALSE, 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(0, 1), c("white", "red")) cm = ColorMapping(col_fun = col_fun) lgd = color_mapping_legend(cm, title = "value", plot = FALSE) value = runif(length(x)) hc = HilbertCurve(1, 1000, level = 9, mode = "pixel", title = "pixel mode", legend = lgd) 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 cm = ColorMapping(levels = c(top8, "Others"), color = c(col[top8], "grey")) lgd = color_mapping_legend(cm, title = "Top8 countries", labels = c(paste(top8, top8_rate), "Others"), plot = FALSE) 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]) ## ------------------------------------------------------------------------ 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") 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)) ## ---- 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"---- 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) }) hc_png(hc, file = 'cnv_all_chromosomes.png') ## ---- 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"---- 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)) hc_png(hc, file = 'cnv_all_chromosomes_overlay.png') ## ---- echo = FALSE, results = "asis"------------------------------------- cat("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", fig.keep = "none"---- 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") hc_png(hc, file = 'cnv_all_chromosomes_border.png') ## ---- echo = FALSE, results = "asis"------------------------------------- cat("

\n") ## ---- eval = grepl("tbi", Sys.info()["nodename"]) & Sys.info()["user"] == "guz", fig.keep = "none", echo = 2:35---- png('cnv_compare.png', width = 800, height = 800) 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 = 2, ncol = 2))) for(i in seq_along(pn)) { pushViewport(viewport(layout.pos.row = (i-1)%/%2 + 1, layout.pos.col = (i-1)%%2 + 1)) plot_curve(pn[i], newpage = FALSE) upViewport() } pushViewport(viewport(layout.pos.row = 2, layout.pos.col = 2)) hc_map(hc, title = "map for 22 chromosomes", newpage = FALSE) upViewport() upViewport() invisible(dev.off()) ## ---- echo = FALSE, results = "asis"------------------------------------- cat("

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