###################################################
### chunk number 1: setup
###################################################
library(ChIPseqTutorial)


###################################################
### chunk number 2: 
###################################################
data(ctcf, package="ChIPseqTutorial")
data(gfp, package="ChIPseqTutorial")


###################################################
### chunk number 3: 
###################################################
ctcf
gfp


###################################################
### chunk number 4: 
###################################################
table(strand(ctcf))
table(chromosome(gfp))


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


###################################################
### chunk number 6: coverage
###################################################
cov.ctcf <- coverage(ctcf, width = mouse.chromlens, extend = 126L)
cov.ctcf
cov.ctcf$chr10


###################################################
### chunk number 7: 
###################################################
islands <- slice(cov.ctcf$chr10, lower = 1)
islands


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

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

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



###################################################
### chunk number 9: 
###################################################
islandReadSummary <- function(cov)
{
    s <- slice(cov, lower = 1)
    tab <- table(viewSums(s) / 150)
    ans <- data.frame(nread = as.numeric(names(tab)), 
                      count = as.numeric(tab))
    ans
}


###################################################
### chunk number 10: 
###################################################
head(islandReadSummary(cov.ctcf$chr10))


###################################################
### chunk number 11: make.groups
###################################################
nread.islands <- lapply(cov.ctcf, islandReadSummary)
nread.islands <- do.call(make.groups, nread.islands)
head(nread.islands)


###################################################
### chunk number 12: xyplot
###################################################
xyplot(log(count) ~ nread | which, data = nread.islands, 
       subset = (nread <= 20), pch = 16, type = c("p", "g"))


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


###################################################
### chunk number 14: 
###################################################
xyplot(log(count) ~ nread | which, data = nread.islands, 
       subset = (nread <= 20),
       pch = 16, 
       panel = function(x, y, ...) {
           panel.grid(h = -1, v = -1)
           panel.lmline(x[1:2], y[1:2], col = "black")
           panel.xyplot(x, y, ...)
       })


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


###################################################
### chunk number 16: 
###################################################
islandDepthSummary <- function(cov)
{
    s <- slice(cov, lower = 1)
    tab <- table(viewMaxs(s))
    ans <- data.frame(depth = as.numeric(names(tab)), count = as.numeric(tab))
    ans
}

depth.islands <- lapply(cov.ctcf, islandDepthSummary)
depth.islands <- do.call(make.groups, depth.islands)

xyplot(log(count) ~ depth | which, depth.islands, 
       subset = (depth <= 20), pch = 16,
       panel = function(x, y, ...) {
           panel.grid(h = -1, v = -1)
           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, ...)
       })



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


###################################################
### chunk number 18: peaks
###################################################
peaks <- slice(cov.ctcf$chr10, lower = 8)
peaks


###################################################
### chunk number 19: 
###################################################
peak.depths <- viewMaxs(peaks)
peak.areas <- viewSums(peaks)
xyplot(peak.areas ~  peak.depths)


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


###################################################
### chunk number 21: 
###################################################
wpeaks <- tail(order(peak.depths), 4)
peaks[wpeaks]


###################################################
### chunk number 22: 
###################################################
coverageplot <- function (peaks, xlab = "Position", ylab = "Coverage", ...)
{
    pos1 <- seq(start(peaks[1]), end(peaks[1]))
    cov1 <- as.integer(peaks[[1]])
    pos1 <- c(head(pos1, 1), pos1, tail(pos1, 1))
    cov1 <- c(0, cov1, 0)
    xyplot(cov1 ~ pos1, ..., panel = panel.polygon, 
           col = "lightgrey", xlab = xlab, ylab = ylab)

}


###################################################
### chunk number 23: 
###################################################
coverageplot(peaks[wpeaks[1]])


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


###################################################
### chunk number 25: all-peaks
###################################################
allPeaks <- slice(cov.ctcf, lower = 8)


###################################################
### chunk number 26: ctcf-counts
###################################################
ctcf.counts <- aggregate(cov.ctcf, allPeaks, sum)


###################################################
### chunk number 27: gfp-counts
###################################################
cov.gfp <- coverage(gfp, width = mouse.chromlens, extend = 126L)
gfp.counts <- aggregate(cov.gfp, allPeaks, sum)


###################################################
### chunk number 28: countsdf
###################################################
counts <- RangedData(allPeaks,
                     CTCF=as.integer(ctcf.counts),
                     GFP=as.integer(gfp.counts))


###################################################
### chunk number 29: 
###################################################
xyplot(asinh(CTCF)~asinh(GFP)|space, as.data.frame(counts))


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


###################################################
### chunk number 31: ctcf-tops
###################################################
rd0 <- counts[counts[["GFP"]] == 0,]["chr10"]
head(rd0[order(rd0[["CTCF"]], decreasing=TRUE),])


###################################################
### chunk number 32: 
###################################################
toLatex(sessionInfo())