library(chipseq)
library(GenomicFeatures)
library(lattice)

data(cstest)
cstest

cstest$ctcf

str(cstest$ctcf$chr10)

library(BSgenome.Mmusculus.UCSC.mm9)
mouse.chromlens <- seqlengths(Mmusculus)
head(mouse.chromlens)

bc <- basesCovered(cstest$ctcf$chr10, shift = 1:250, seqLen = 24)
xyplot(covered ~ mu, bc, type = "l")

ext <- extendReads(cstest$ctcf$chr10, seqLen = 150)
head(ext)

cov <- coverage(ext, width = mouse.chromlens["chr10"])
cov

islands <- slice(cov, lower = 1)
islands

viewSums(head(islands))
viewMaxs(head(islands))

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

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

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

head(islandReadSummary(cstest$ctcf$chr10))

nread.islands <- gdapply(cstest, islandReadSummary)
nread.islands <- as(nread.islands, "data.frame")
head(nread.islands)

xyplot(log(count) ~ nread | chromosome, nread.islands, 
       subset = (sample == "ctcf" & nread <= 20), 
       layout = c(3, 1), pch = 16, type = c("p", "g"))

plot(trellis.last.object())

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

plot(trellis.last.object())

islandDepthSummary <- function(x)
{
    g <- extendReads(x, seqLen = 150)
    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 | chromosome, depth.islands, 
       subset = (sample == "ctcf" & depth <= 20), 
       layout = c(3, 1), 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, ...)
       })

peaks <- slice(cov, lower = 8)
peaks

peak.depths <- viewMaxs(peaks)

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

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

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

coverageplot(peaks.pos[wpeaks[1]], peaks.neg[wpeaks[1]])
coverageplot(peaks.pos[wpeaks[2]], peaks.neg[wpeaks[2]])
coverageplot(peaks.pos[wpeaks[3]], peaks.neg[wpeaks[3]])
coverageplot(peaks.pos[wpeaks[4]], peaks.neg[wpeaks[4]])

subdata <- 
    subset(do.call(make.groups, cstest$ctcf$chr10),
           data > start(peaks)[wpeaks[1]] - 200 & data < end(peaks)[wpeaks[1]] + 200)
densityplot(~data, data = subdata, groups = which, auto.key = TRUE)
stripplot(which ~ data, data = subdata, jitter = TRUE)


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

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

xyplot(log1p(sums2) ~ log1p(sums1) | chromosome, data = peakSummary, 
       type = c("p", "g"), alpha = 0.5, aspect = "iso", pch = ".", cex = 3)

sessionInfo()