###################################################
### chunk number 1: setup
###################################################

library(chipseq)
library(GenomicFeatures.Mmusculus.UCSC.mm9)
library(lattice)



###################################################
### chunk number 2: 
###################################################
data(cstest)
cstest


###################################################
### chunk number 3: 
###################################################
cstest$ctcf


###################################################
### chunk number 4: 
###################################################
str(cstest$ctcf$chr10)


###################################################
### chunk number 5: 
###################################################
library(BSgenome.Mmusculus.UCSC.mm9)
mouse.chromlens <- seqlengths(Mmusculus)
head(mouse.chromlens)


###################################################
### chunk number 6: 
###################################################
ext <- extendReads(cstest$ctcf$chr10, seqLen = 200)
head(ext)


###################################################
### chunk number 7: 
###################################################
cov <- coverage(ext, width = mouse.chromlens["chr10"])
cov


###################################################
### chunk number 8: 
###################################################
islands <- slice(cov, lower = 1)
islands


###################################################
### chunk number 9: 
###################################################
viewSums(head(islands))
viewMaxs(head(islands))

nread.tab <- table(viewSums(islands) / 200)
depth.tab <- table(viewMaxs(islands))

head(nread.tab, 10)
head(depth.tab, 10)



###################################################
### chunk number 10: 
###################################################
islandReadSummary <- function(x)
{
    g <- extendReads(x, seqLen = 200)
    s <- slice(coverage(g), lower = 1)
    tab <- table(viewSums(s) / 200)
    ans <- data.frame(nread = as.numeric(names(tab)), count = as.numeric(tab))
    ans
}


###################################################
### chunk number 11: 
###################################################
head(islandReadSummary(cstest$ctcf$chr10))


###################################################
### chunk number 12: 
###################################################
nread.islands <- gdapply(cstest, islandReadSummary)
nread.islands <- as(nread.islands, "data.frame")
head(nread.islands)


###################################################
### chunk number 13: 
###################################################
xyplot(log(count) ~ nread | sample, nread.islands, 
       subset = (chromosome == "chr10" & nread <= 40), 
       layout = c(1, 2), pch = 16, type = c("p", "g"))



###################################################
### chunk number 14: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 15: 
###################################################
xyplot(log(count) ~ nread | sample, nread.islands, 
       subset = (chromosome == "chr10" & nread <= 40), 
       layout = c(1, 2), pch = 16, type = c("p", "g"),
       panel = function(x, y, ...) {
           panel.lmline(x[1:2], y[1:2], col = "black")
           panel.xyplot(x, y, ...)
       })



###################################################
### chunk number 16: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 17: 
###################################################
islandDepthSummary <- function(x)
{
    g <- extendReads(x, seqLen = 200)
    s <- slice(coverage(g), lower = 1)
    tab <- table(viewMaxs(s))
    ans <- data.frame(depth = as.numeric(names(tab)), count = as.numeric(tab))
    ans
}

depth.islands <- gdapply(cstest, islandDepthSummary)
depth.islands <- as(depth.islands, "data.frame")

xyplot(log(count) ~ depth | sample, depth.islands, 
       subset = (chromosome == "chr10" & depth <= 20), 
       layout = c(1, 2), pch = 16, type = c("p", "g"),
       panel = function(x, y, ...) {
           lambda <- 2 * exp(y[2]) / exp(y[1])
           null.est <- function(xx) {
               xx * log(lambda) - lambda - lgamma(xx + 1)
           }
           log.N.hat <- null.est(1) - y[1]
           panel.lines(1:10, -log.N.hat + null.est(1:10), col = "black")
           panel.xyplot(x, y, ...)
       })

## depth.islands <- summarizeReads(cstest, summary.fun = islandDepthSummary)



###################################################
### chunk number 18: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 19: 
###################################################
peaks <- slice(cov, lower = 8)
peaks


###################################################
### chunk number 20: 
###################################################
peak.depths <- viewMaxs(peaks)

cov.pos <- coverage(extendReads(cstest$ctcf$chr10, strand = "+", seqLen = 200), 
                    width = mouse.chromlens["chr10"])
cov.neg <- coverage(extendReads(cstest$ctcf$chr10, strand = "-", seqLen = 200), 
                    width = mouse.chromlens["chr10"])

peaks.pos <- copyIRanges(peaks, cov.pos)
peaks.neg <- copyIRanges(peaks, cov.neg)

wpeaks <- tail(order(peak.depths), 4)
wpeaks



###################################################
### chunk number 21: 
###################################################
coverageplot(peaks.pos[wpeaks[1]], peaks.neg[wpeaks[1]])


###################################################
### chunk number 22: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 23: 
###################################################
coverageplot(peaks.pos[wpeaks[2]], peaks.neg[wpeaks[2]])


###################################################
### chunk number 24: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 25: 
###################################################
coverageplot(peaks.pos[wpeaks[3]], peaks.neg[wpeaks[3]])


###################################################
### chunk number 26: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 27: 
###################################################
coverageplot(peaks.pos[wpeaks[4]], peaks.neg[wpeaks[4]])


###################################################
### chunk number 28: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 29: 
###################################################

extRanges <- gdapply(cstest, extendReads, seqLen = 200)

peakSummary <-
    diffPeakSummary(extRanges$gfp, extRanges$ctcf,
                    chrom.lens = mouse.chromlens, lower = 10)

xyplot(asinh(sums2) ~ asinh(sums1) | chromosome,
       data = peakSummary, 
       ## subset = (chromosome %in% c("chr1", "chr2")),
       panel = function(x, y, ...) {
           panel.xyplot(x, y, ...)
           panel.abline(median(y - x), 1)
       },
       type = c("p", "g"), alpha = 0.5, aspect = "iso")



###################################################
### chunk number 30: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 31: 
###################################################
peakSummary <- 
    within(peakSummary,
       {
           diffs <- asinh(sums2) - asinh(sums1)
           resids <- (diffs - median(diffs)) / mad(diffs)
           up <- resids > 2
           down <- resids < -2
       })
head(peakSummary)


###################################################
### chunk number 32: 
###################################################
gregions <- transcripts(genes = geneMouse(), proximal = 2000)
## gregions$gene <- as.character(gregions$gene)
gregions


###################################################
### chunk number 33: 
###################################################
ans <- contextDistribution(peakSummary, gregions)
head(ans)


###################################################
### chunk number 34: 
###################################################
sumtab <- 
    as.data.frame(rbind(total = xtabs(total ~ type, ans),
                        promoter = xtabs(promoter ~ type, ans),
                        threeprime = xtabs(threeprime ~ type, ans),
                        upstream = xtabs(upstream ~ type, ans),
                        downstream = xtabs(downstream ~ type, ans),
                        gene = xtabs(gene ~ type, ans)))

cbind(sumtab, ratio = round(sumtab[, "down"] /  sumtab[, "up"], 3))



###################################################
### chunk number 35: 
###################################################
sessionInfo()