###################################################
### chunk number 1: load-genome
###################################################
library(BSgenome.Mmusculus.UCSC.mm9)


###################################################
### chunk number 2: define-patterns
###################################################
NN <- mkAllStrings(c("A","C","G","T"), 2)
motifs <- DNAStringSet(paste("CA",NN,"TG",sep=""))
pD <- PDict(motifs)
motifs <- as.character(motifs)


###################################################
### chunk number 3: matching-function
###################################################
findMotifs <- function(chr) {
  mindex <- matchPDict(pD, chr)
  seq <- rep(motifs, countIndex(mindex))
  gd <- RangedData(unlist(mindex), seq)
  gd[order(start(gd)),]
}


###################################################
### chunk number 4: bsapply
###################################################
params <- new("BSParams", X = Mmusculus, 
              FUN = findMotifs, 
              exclude = "[_MXY]")
motifLocs <- do.call("c", bsapply(params))


###################################################
### chunk number 5: find-peaks
###################################################
load("../data/alignedLocs.rda")
library(chipseq)
sample1 <- alignedLocs$sample1
extended <- extendReads(as.list(sample1))
callPeaks <- function(chr) {
  cov <- coverage(chr, start = 1, 
                  end = max(end(chr)))
  slice(cov, 8)
}
peaks <- lapply(extended, callPeaks)


###################################################
### chunk number 6: find-promoters
###################################################
data(geneMouse)
regions <- genomic_regions(geneMouse)
promRanges <- IRanges(regions$promoter.start, 
                      regions$promoter.end)
promoters <- split(promRanges, regions$chrom)


###################################################
### chunk number 7: define-filters
###################################################
overlapFilter <- function(x) { 
  function(rd) 
    ranges(rd)[[1]] %in% x[[names(rd)]]
}
promoterFilter <- overlapFilter(promoters)
peakFilter <- overlapFilter(peaks)
filters <- list(promoter = promoterFilter, 
                peak = peakFilter)
rules <- FilterRules(filters, active = FALSE)


###################################################
### chunk number 8: define-counter
###################################################
count_motifs <- function(rd) {
  nn <- substring(rd[["seq"]], 3, 4)
  df <- as.data.frame(table(factor(nn, NN)))
  colnames(df) <- c("seq", "count")
  df
}


###################################################
### chunk number 9: define-reducer
###################################################
reduce_counts <- function(counts) {
  counts <- do.call("rbind", counts)
  counts <- aggregate(counts[,2,drop=FALSE], 
                      list(seq = counts$seq), sum)
  counts$freq <- counts$count / sum(counts$count)
  counts
}


###################################################
### chunk number 10: invoke-rdapply
###################################################
rda <- RDApplyParams(motifLocs, count_motifs, 
                     filterRules = rules,
                     reducerFun = reduce_counts)


###################################################
### chunk number 11: count-all
###################################################
allCounts <- rdapply(rda)


###################################################
### chunk number 12: count-peaks
###################################################
active(filterRules(rda))["peak"] <- TRUE
peakCounts <- rdapply(rda)


###################################################
### chunk number 13: count-promoter
###################################################
active(filterRules(rda))["promoter"] <- TRUE
promoterCounts <- rdapply(rda)


###################################################
### chunk number 14: plot-all
###################################################
barplot(allCounts$freq, names.arg = allCounts$seq, main = "All motifs")


###################################################
### chunk number 15: plot-peaks
###################################################
barplot(peakCounts$freq, names.arg = peakCounts$seq, main = "All motifs")


###################################################
### chunk number 16: plot-promoter
###################################################
barplot(promoterCounts$freq, names.arg = promoterCounts$seq,
main = "Promoter peak motifs")


###################################################
### chunk number 17: plot-compare
###################################################
combinedFreqs <- rbind(all = allCounts$freq, peak = peakCounts$freq,
promoter = promoterCounts$freq)
barplot(combinedFreqs, beside = TRUE, names.arg = promoterCounts$seq,
main = "Motif freqs", legend.text = TRUE)


###################################################
### chunk number 18: session-info
###################################################
sessionInfo()