--- a
+++ b/b_DownstreamAnalysisScript/bulkATACana_3_annotatePeak.R
@@ -0,0 +1,54 @@
+
+# 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)
+
+# analysis
+library(ChIPseeker)
+library(GenomicFeatures)
+
+# * 2. Load data ----------------------------------------------------------
+
+txdb <- makeTxDbFromGFF("../../data/hg19genes.gtf")
+
+peakMeta <- readRDS("peakMeta.rds")
+peakMeta[, id := paste(Chr, Start, End, sep = "_")][]
+
+peakGroup <- readRDS("peakGroup.rds")
+peakIndex <- map(peakGroup, ~ {match(.x, peakMeta$id)})
+
+# * 3. Analyze ------------------------------------------------------------
+
+homerPeak <- peakMeta
+homerPeak[, Chr := str_c("chr", Chr)][]
+
+fwrite(homerPeak, "allPeak.homer", sep = "\t", col.names = F)
+
+peakGRlist <- map(peakIndex, ~ {
+  GRanges(
+    seqnames = peakMeta[.x, Chr],
+    ranges = IRanges(
+      start = peakMeta[.x, Start],
+      end = peakMeta[.x, End]))})
+
+peakAnnoList <- map(peakGRlist, annotatePeak, TxDb = txdb, tssRegion = c(-2000, 2000))
+
+map(peakGRlist, ~ {as.data.table(.x)[, c("width", "strand") := NULL][, id := 1:.N][, seqnames := str_c("chr", seqnames)]}) %>%
+  iwalk(~ fwrite(.x, str_c("peakGroup_", .y, ".bed"), sep = "\t", col.names = F))
+
+saveRDS(peakAnnoList, "peakAnnoList.rds")