EnrichedHeatmap example code

EnrichedHeatmap_example.R

#/ A minimal example on how to use EnrichedHeatmap together with rtracklayer.
#/ Required data are at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111902
#/ Go for the files :
#/ "GSM3045250_N_ATAC_Kaech_1_peaks.broadPeak.gz" and (genomic regions in BED-like format)
#/ "GSM3045250_N_ATAC_Kaech_1.bw"                     (the bigwig files with read counts)

require(EnrichedHeatmap)
require(rtracklayer)
require(circlize)
require(data.table)

#/ Load a BED file into R with data.table::fread(),
#/ then convert to GRanges format with GenomicRanges::GRanges():
targets <- makeGRangesFromDataFrame(
  df = fread("GSM3045250_N_ATAC_Kaech_1_peaks.broadPeak.gz", header = FALSE, data.table = FALSE),
  seqnames.field = "V1", start.field = "V2", end.field = "V3")

#/ For this tutorial we take the first 10000 regions rather than the full file:
targets <- head(targets, 10000)

#/ We take the center of each region/peak and want to extend by 5kb each direction:
ExtendSize <- 5000
targets.extended  <- resize(targets, fix = "center", width = ExtendSize*2)

#/ We load the relevant parts of the bigwig file into R. 
#/ This is more efficient than reading the entire file. 
BigWig <- rtracklayer::import("GSM3045250_N_ATAC_Kaech_1.bw", 
                              format = "BigWig", 
                              selection = BigWigSelection(targets.extended))

#/ Create the normalizedMatrix that EnrichedHeatmap accepts as input.
#/ We use the targets center (width=1) because from what I understand normalizeMatrix
#/ does not allow to turn off its "extend" option. Therefore we trick it by simply
#/ providing the peak centers and then let the function extend it by our predefined window size.
normMatrix <- normalizeToMatrix(signal = BigWig, 
                                target = resize(targets, fix = "center", width = 1), 
                                background = 0, 
                                keep = c(0, 0.99),      #/ minimal value to the 99th percentile
                                target_ratio = 0,
                                mean_mode = "w0",       #/ see ?EnrichedHeatmap on other options
                                value_column = "score", #/ = the name of the 4th column of the bigwig
                                extend = ExtendSize)

#/ Make a color gradient that covers the range of normMatrix from 0 to the 99th percentile.
#/ The percentile avoids outliers to skew the heatmap:
col_fun = circlize::colorRamp2(quantile(normMatrix, c(0, .99)), c("darkblue", "darkgoldenrod1"))

#/ heatmap function:
EH <- EnrichedHeatmap( mat = normMatrix, 
                       pos_line = FALSE, #/ no dashed lines around the start
                       border = FALSE,   #/ no box around heatmap
                       col = col_fun,    #/ color gradients from above
                       column_title = "Nice Heatmap", #/ column title 
                       column_title_gp = gpar(fontsize = 15, fontfamily = "sans"),
                       use_raster = TRUE, raster_quality = 10, raster_device = "png",
                       #/ turn off background colors
                       rect_gp = gpar(col = "transparent"), 
                       #/ legend options:
                       heatmap_legend_param = list(
                         legend_direction = "horizontal",
                         title = "normalized counts"),
                       #/ options for the profile plot on top of the heatmap:
                       top_annotation = HeatmapAnnotation(
                         enriched = anno_enriched(
                           gp = gpar(col = "black", lty = 1, lwd=2),
                           col="black")
                       )
) #/ end of EnrichedHeatmap function

#/ Save as pdf to disk:
pdf("EnrichedHeatmap.pdf")

draw(EH,                                 #/ plot the heatmap from above 
     heatmap_legend_side = "bottom",     #/ we want the legend below the heatmap
     annotation_legend_side = "bottom",  #/ legend on the bottom side
     padding = unit(c(4, 4, 4, 4), "mm") #/ some padding to avoid labels beyond plot borders
)

dev.off()