--- a
+++ b/c_VisualizationScript/Visulz_bulkATAC_heatmapPeak.R
@@ -0,0 +1,75 @@
+
+# MESSAGE -----------------------------------------------------------------
+#
+# author: Yulin Lyu
+# email: lvyulin@pku.edu.cn
+#
+# require: R whatever
+#
+# ---
+
+# * 1. Load packages ------------------------------------------------------
+
+setwd("exampleData/ATAC")
+
+# grammar
+library(tidyverse)
+library(magrittr)
+library(glue)
+library(data.table)
+library(GenomicRanges)
+
+# graphics
+library(ComplexHeatmap)
+library(circlize)
+library(RColorBrewer)
+library(ggsci)
+library(scales)
+
+# * 2. Load data ----------------------------------------------------------
+
+diffList <- readRDS("diffList.rds")
+peakMeta <- readRDS("peakMeta.rds")
+peakMtx <- readRDS("peakMtx.rds")
+
+diffPeak <- map(diffList, ~ {as.data.table(.x)[abs(Fold) > 1, paste(seqnames, start, end, sep = "_")]}) %>%
+  unlist() %>% unname() %>% unique()
+
+peakMeta[, id := paste(Chr, Start, End, sep = "_")][]
+
+diffMtx <- as.matrix(peakMtx[peakMeta$id %in% diffPeak]) %>% set_rownames(peakMeta[id %in% diffPeak, id])
+
+# * 3. Plot ---------------------------------------------------------------
+
+dir.create("graphics")
+
+heatData <- apply(diffMtx, 1, scale) %>% t() %>% set_colnames(colnames(diffMtx))
+
+colorFun <- colorRamp2(seq(2, -2, len = 9), brewer.pal(9, "RdBu"))
+
+p <- Heatmap(
+  heatData[, c(3:5, 1:2, 6:8)], col = colorFun, border = F,
+  cluster_rows = T, cluster_columns = F,
+  show_row_names = F, show_column_names = T,
+  show_row_dend = F, show_column_dend = F,
+  column_names_rot = 45,
+  row_km = 5,
+  row_km_repeats = 100,
+  left_annotation = rowAnnotation(
+    type = anno_block(
+      width = unit(.1, "in"), gp = gpar(fill = pal_nejm()(5)))),
+  heatmap_legend_param = list(
+    title = "Scaled expression",
+    title_position = "lefttop-rot",
+    legend_height = unit(2, "in")),
+  width = unit(3, "in"),
+  height = unit(6, "in"))
+
+png("graphics/heatmapPeak.png", width = 4, height = 7, units = "in", res = 300)
+p
+dev.off()
+
+peakOrder <- row_order(p)
+peakGroup <- map(peakOrder, ~ {rownames(diffMtx)[.x]})
+
+saveRDS(peakGroup, "peakGroup.rds")