Author: Zuguang Gu ( z.gu@dkfz.de )
Date: 2016-02-27
This supplementary contains source code that generates the six plots from figure 1 in the main manuscript.
library(HilbertCurve)
library(ComplexHeatmap)
library(circlize)
Figure 1A:
set.seed(123)
x = sort(sample(100, 20))
s = x[1:10*2 - 1]
e = x[1:10*2]
ir = IRanges(s, e)
labels = sample(letters, length(ir), replace = TRUE)
hc = HilbertCurve(1, 100, level = 4, title = "A) Combine low-level graphics")
hc_segments(hc, IRanges(1, 100))
hc_rect(hc, ir, gp = gpar(fill = rand_color(length(ir), transparency = 0.8)))
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))
Figure 1B, how to get the data is explained in SupplS3.
chr1_len = 249250621
load(system.file("extdata", "zebrafish_net.RData", package = "HilbertCurve"))
# extract the complement regions that are not conserved
seqlengths(zebrafish) = chr1_len
nonzebrafish = gaps(zebrafish); nonzebrafish = nonzebrafish[strand(nonzebrafish) == "*"]
gr = c(zebrafish, nonzebrafish)
# define associated color for regions
col = c(rep("red", length(zebrafish)), rep("yellow", length(nonzebrafish)))
cm = ColorMapping(col_fun = colorRamp2(c(0, 1), c("yellow", "red")))
lgd = color_mapping_legend(cm, title = "Conservation", at = c(0, 0.2, 0.4, 0.6, 0.8, 1),
labels = c("0%", "20%", "40%", "60%", "80%", "100%"), plot = FALSE)
hc = GenomicHilbertCurve(chr = "chr1", level = 6,
title = "B) Sequence conservation between zebrafish and human",
legend = lgd)
hc_points(hc, gr, np = 3, gp = gpar(col = NA, fill = col))
Figure 1C. UCSD.Lung.H3K36me3.STL002.bed is downloaded from Roadmap.
load(system.file("extdata", "refseq_chr1.RData", package = "HilbertCurve"))
# this line is only used to extract regions on chromosome 1
df = read.table(pipe("awk '$5>0 && $1==\"chr1\"' ~/HilbertCurveTest/UCSD.Lung.H3K36me3.STL002.bed"), sep = "\t")
# default color theme
col_fun = colorRamp2(c(0, quantile(df[[5]], 0.99)), c("white", "red"))
# color theme for the regions intersected with gene bodies
col_fun_new = colorRamp2(c(0, quantile(df[[5]], 0.99)), c("white", "purple"))
# generate legends
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)
hc = GenomicHilbertCurve(chr = "chr1", level = 9, mode = "pixel",
title = "C) Intensity of H3K36me3 mark, chromosome 1", 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)
# change to the white-purple color theme
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)
})
# you can save the curve as 512x512 pixel PNG file
# hc_png(hc, file = "H3K36me3_chr1.png")
Figure 1D. UCSD.Lung.Bisulfite-Seq.STL002.bed is downloaded from Roadmap.
df = read.table(pipe("awk '$1==\"chr1\"' ~/HilbertCurveTest/UCSD.Lung.Bisulfite-Seq.STL002.bed"), sep = "\t")
# color theme
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)
hc = GenomicHilbertCurve(chr = "chr1", level = 9, mode = "pixel",
title = "D) Methylation, chromosome 1", legend = lgd)
# note 'absolute' mode is used here
hc_layer(hc, df, col = col_fun(df[[5]]), mean_mode = "absolute")
# you can save the curve as 512x512 pixel PNG file
# hc_png(hc, file = "methylation_chr1.png")
Figure 1E and 1F. Copy number variation data is downloaded from Database of Genomic Variants, CNV map.
cm = ColorMapping(levels = c("gain", "loss", "both"), color = c("red", "green", "purple"))
lgd = color_mapping_legend(cm, title = "CNV", plot = FALSE)
hc = GenomicHilbertCurve(chr = paste0("chr", 1:22), level = 9, mode = "pixel",
title = "E) Copy number variation, 22 chromosomes", legend = lgd)
# layer for gain events
df = read.table("~/HilbertCurveTest/Stringent.Gain.hg19.2015-02-03.txt",
header = TRUE, stringsAsFactors = FALSE)
hc_layer(hc, df, col = "red")
# layer for loss events
df = read.table("~/HilbertCurveTest/Stringent.Loss.hg19.2015-02-03.txt",
header = TRUE, stringsAsFactors = FALSE)
hc_layer(hc, df, col = "green",
overlay = function(r0, g0, b0, r, g, b, alpha) {
l = !(r0 == 1 & g0 == 1 & b0 == 1)
# change overlapped regions to purple
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")
# you can save the curve as 512x512 pixel PNG file
# hc_png(hc, file = "copy_number_variation.png")
# the map of chromosome regions to accompany figure 1E
hc_map(hc, title = "F) Map for 22 chromosomes")
sessionInfo()
## R version 3.2.2 (2015-08-14)
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## locale:
## [1] C
##
## attached base packages:
## [1] stats4 parallel methods grid stats graphics grDevices
## [8] utils datasets base
##
## other attached packages:
## [1] circlize_0.3.4 ComplexHeatmap_1.6.0 HilbertCurve_1.1.3
## [4] GenomicRanges_1.22.3 GenomeInfoDb_1.6.3 IRanges_2.4.6
## [7] S4Vectors_0.9.19 BiocGenerics_0.16.1
##
## loaded via a namespace (and not attached):
## [1] whisker_0.3-2 knitr_1.12.3 XVector_0.10.0
## [4] magrittr_1.5 zlibbioc_1.16.0 colorspace_1.2-6
## [7] lattice_0.20-33 rjson_0.2.15 stringr_1.0.0
## [10] tools_3.2.2 png_0.1-7 RColorBrewer_1.1-2
## [13] formatR_1.2.1 HilbertVis_1.28.0 GlobalOptions_0.0.8
## [16] dendextend_1.1.2 shape_1.4.2 evaluate_0.8
## [19] stringi_1.0-1 GetoptLong_0.1.1