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

library(chipseq)
library(lattice)



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

load("../data/alignedLocs.rda")
alignedLocs


###################################################
### chunk number 3: 
###################################################

alignedLocs$sample1


###################################################
### chunk number 4: 
###################################################

str(head(as(alignedLocs$sample1, "list"), 4))


###################################################
### chunk number 5:  eval=FALSE
###################################################

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


###################################################
### chunk number 6: 
###################################################

head(mouse.chromlens)


###################################################
### chunk number 7: 
###################################################

ext <- extendReads(alignedLocs$sample1$chr5, seqLen = 200)
head(ext)


###################################################
### chunk number 8: 
###################################################

cov <- coverage(ext, start = 1, end = mouse.chromlens["chr5"])
cov


###################################################
### chunk number 9: 
###################################################

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


###################################################
### chunk number 10: 
###################################################

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 11: 
###################################################

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


###################################################
### chunk number 12: 
###################################################

head(islandReadSummary(alignedLocs$sample1$chr5))


###################################################
### chunk number 13: 
###################################################

nread.islands <- gdApply(alignedLocs, islandReadSummary)
nread.islands <- as(nread.islands, "data.frame")
head(nread.islands)


###################################################
### chunk number 15: 
###################################################

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



###################################################
### chunk number 17: 
###################################################

xyplot(log(count) ~ nread | sample, nread.islands, 
       subset = (chromosome == "chr5" & nread <= 40), 
       layout = c(1, 3), 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 19: 
###################################################

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

depth.islands <- gdApply(alignedLocs, islandDepthSummary)
depth.islands <- as(depth.islands, "data.frame")

xyplot(log(count) ~ depth | sample, depth.islands, 
       subset = (chromosome == "chr5" & depth <= 20), 
       layout = c(1, 3), 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, ...)
       })


###################################################
### chunk number 21: 
###################################################

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


###################################################
### chunk number 22: 
###################################################

peak.depths <- viewMaxs(peaks)

cov.pos <- coverage(extendReads(alignedLocs$sample1$chr5, strand = "+", seqLen = 200), 
                    start = 1, end = mouse.chromlens["chr5"])
cov.neg <- coverage(extendReads(alignedLocs$sample1$chr5, strand = "-", seqLen = 200), 
                    start = 1, end = mouse.chromlens["chr5"])

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

which(peak.depths >= 40)



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

plotPeak(peaks.pos[6], peaks.neg[6])


###################################################
### chunk number 25: 
###################################################

plotPeak(peaks.pos[96], peaks.neg[96])


###################################################
### chunk number 27: 
###################################################
plotPeak(peaks.pos[126], peaks.neg[126])


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

plotPeak(peaks.pos[173], peaks.neg[173])


###################################################
### chunk number 31: 
###################################################

extRanges <- gdApply(alignedLocs, extendReads, seqLen = 200)

peakSummary <-
    diffPeakSummary(extRanges$sample1, extRanges$sample2,
                    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 33: 
###################################################

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


###################################################
### chunk number 34: 
###################################################

data(geneMouse)
gregions <- genomic_regions(genes = geneMouse, proximal = 2000)
gregions$gene <- as.character(gregions$gene)
str(gregions)


###################################################
### chunk number 35: 
###################################################

ans <- contextDistribution(peakSummary, gregions)
head(ans)


###################################################
### chunk number 36: 
###################################################

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 37: 
###################################################

sessionInfo()