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