###################################################
### chunk number 1: 
###################################################
library(ShortRead)


###################################################
### chunk number 2: setup
###################################################

library(ShortRead)
library(lattice)



###################################################
### chunk number 3: 
###################################################
load("../data/ctcf.rda")
load("../data/gfp.rda")


###################################################
### chunk number 4: 
###################################################
ctcf
gfp


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


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


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


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


###################################################
### chunk number 9: 
###################################################
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 10: 
###################################################
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 11: 
###################################################
head(islandReadSummary(cov.ctcf$chr10))


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


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


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


###################################################
### chunk number 15: 
###################################################
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 16: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 17: 
###################################################
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 18: 
###################################################
plot(trellis.last.object())


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


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



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


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


###################################################
### chunk number 23: 
###################################################

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 24: 
###################################################
coverageplot(peaks[wpeaks[1]])


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


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