library(DiffBind)
diffbind <- system.file(package="DiffBind", "extra")
dir(diffbind)
csv <- read.csv(file.path(diffbind, "tamoxifen.csv"))

fls <- dir(file.path(diffbind, "peaks"), full=TRUE)
names(fls) <- sub(".bed.gz", "", basename(fls))
length(fls)

library(rtracklayer)
peak1 <- import(fls[1])
peak1
hist(width(peak1))

peaks <- lapply(fls, import)
all <- reduce(do.call(c, unname(peaks)))

occupied <- sapply(peaks, countOverlaps, query=all)
head(occupied)
colSums(occupied)
plot(table(rowSums(occupied != 0)))

library(TxDb.Hsapiens.UCSC.hg18.knownGene)
gn <- genes(TxDb.Hsapiens.UCSC.hg18.knownGene)
chr18 <- gn[seqnames(gn) == "chr18"]

sum(countOverlaps(all, chr18) != 0)
## same as sum(all %over% chr18)
sum(countOverlaps(all, chr18, type="within") != 0)
## same as sum(all %within% chr18)

upstream <- flank(chr18, 5000)
sum(all %over% upstream)

downstream <- flank(chr18, 5000, start=FALSE)
sum(all %over% downstream)

near <- resize(chr18, width(chr18) + 10000, fix="center")
sum(all %over% near)

dist <- distanceToNearest(all[all %over% upstream], starts)
plot(density(mcols(dist)$distance))

mcols(all)$is_near_gene <- all %over% nearby