### R code from vignette source 'vignettes/chipseq/inst/doc/Workflow.Rnw'

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

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



###################################################
### code chunk number 2: preprocess (eval = FALSE)
###################################################
## qa_list <- lapply(sampleFiles, qa)
## report(do.call(rbind, qa_list))
## ## spend some time evaluating the QA report, then procede
## filter <- compose(chipseqFilter(), alignQualityFilter(15))
## cstest <- seqapply(sampleFiles, function(file) {
##   as(readAligned(file, filter), "GRanges")
## })
## cstest <- cstest[seqnames(cstest) %in% c("chr10", "chr11", "chr12")]


###################################################
### code chunk number 3: Workflow.Rnw:85-87
###################################################
data(cstest)
cstest


###################################################
### code chunk number 4: convert-cstest (eval = FALSE)
###################################################
## ## code used to convert the GenomeDataList to a GRangesList
## cstest <- seqapply(cstest, function(gd) {
##   gr <- do.call(c, lapply(names(gd), function(chr) {
##     pos <- gd[[chr]]
##     starts <- c(pos[["-"]] - 23L, pos[["+"]])
##     GRanges(chr, IRanges(starts, width = 24), 
##             rep(c("-", "+"), elementLengths(pos)))
##   }))
## })


###################################################
### code chunk number 5: Workflow.Rnw:105-106
###################################################
cstest$ctcf


###################################################
### code chunk number 6: Workflow.Rnw:115-116
###################################################
library(BSgenome.Mmusculus.UCSC.mm9)


###################################################
### code chunk number 7: estimate.mean.fraglen
###################################################
fraglen <- estimate.mean.fraglen(cstest$ctcf)
fraglen[!is.na(fraglen)]


###################################################
### code chunk number 8: Workflow.Rnw:151-153
###################################################
ctcf.ext <- resize(cstest$ctcf, width = 200)
ctcf.ext


###################################################
### code chunk number 9: Workflow.Rnw:164-166
###################################################
cov.ctcf <- coverage(ctcf.ext)
cov.ctcf


###################################################
### code chunk number 10: Workflow.Rnw:176-178
###################################################
islands <- slice(cov.ctcf, lower = 1)
islands


###################################################
### code chunk number 11: Workflow.Rnw:183-191
###################################################
viewSums(islands)
viewMaxs(islands)

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

nread.tab[,1:10]
depth.tab[,1:10]


###################################################
### code chunk number 12: Workflow.Rnw:207-217
###################################################
islandReadSummary <- function(x)
{
    g <- resize(x, 200)
    s <- slice(coverage(g), lower = 1)
    tab <- table(viewSums(s) / 200)
    df <- DataFrame(tab)
    colnames(df) <- c("chromosome", "nread", "count")
    df$nread <- as.integer(df$nread)
    df
}


###################################################
### code chunk number 13: Workflow.Rnw:221-222
###################################################
head(islandReadSummary(cstest$ctcf))


###################################################
### code chunk number 14: Workflow.Rnw:227-230
###################################################
nread.islands <- seqapply(cstest, islandReadSummary)
nread.islands <- stack(nread.islands, "sample")
nread.islands


###################################################
### code chunk number 15: Workflow.Rnw:236-239
###################################################
xyplot(log(count) ~ nread | sample, as.data.frame(nread.islands), 
       subset = (chromosome == "chr10" & nread <= 40), 
       layout = c(1, 2), pch = 16, type = c("p", "g"))


###################################################
### code chunk number 16: Workflow.Rnw:242-243
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 17: Workflow.Rnw:259-266
###################################################
xyplot(log(count) ~ nread | sample, as.data.frame(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, ...)
       })


###################################################
### code chunk number 18: Workflow.Rnw:269-270
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 19: Workflow.Rnw:277-306
###################################################
islandDepthSummary <- function(x)
{
  g <- resize(x, 200)
  s <- slice(coverage(g), lower = 1)
  tab <- table(viewMaxs(s) / 200)
  df <- DataFrame(tab)
  colnames(df) <- c("chromosome", "depth", "count")
  df$depth <- as.integer(df$depth)
  df
}

depth.islands <- seqapply(cstest, islandDepthSummary)
depth.islands <- stack(depth.islands, "sample")

xyplot(log(count) ~ depth | sample, as.data.frame(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)



###################################################
### code chunk number 20: Workflow.Rnw:309-310
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 21: islandDepthPlot
###################################################
islandDepthPlot(cov.ctcf)


###################################################
### code chunk number 22: peakCutoff
###################################################
peakCutoff(cov.ctcf, fdr = 0.0001)


###################################################
### code chunk number 23: Workflow.Rnw:337-339
###################################################
peaks.ctcf <- slice(cov.ctcf, lower = 8)
peaks.ctcf


###################################################
### code chunk number 24: peakSummary
###################################################
peaks <- peakSummary(peaks.ctcf)


###################################################
### code chunk number 25: Workflow.Rnw:357-368
###################################################
peak.depths <- viewMaxs(peaks.ctcf)

cov.pos <- coverage(ctcf.ext[strand(ctcf.ext) == "+"])
cov.neg <- coverage(ctcf.ext[strand(ctcf.ext) == "-"])

peaks.pos <- Views(cov.pos, ranges(peaks.ctcf))
peaks.neg <- Views(cov.neg, ranges(peaks.ctcf))

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



###################################################
### code chunk number 26: Workflow.Rnw:375-376
###################################################
coverageplot(peaks.pos$chr10[wpeaks[1]], peaks.neg$chr10[wpeaks[1]])


###################################################
### code chunk number 27: Workflow.Rnw:378-379
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 28: Workflow.Rnw:382-383
###################################################
coverageplot(peaks.pos$chr10[wpeaks[2]], peaks.neg$chr10[wpeaks[2]])


###################################################
### code chunk number 29: Workflow.Rnw:385-386
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 30: Workflow.Rnw:391-392
###################################################
coverageplot(peaks.pos$chr10[wpeaks[3]], peaks.neg$chr10[wpeaks[3]])


###################################################
### code chunk number 31: Workflow.Rnw:394-395
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 32: Workflow.Rnw:398-399
###################################################
coverageplot(peaks.pos$chr10[wpeaks[4]], peaks.neg$chr10[wpeaks[4]])


###################################################
### code chunk number 33: Workflow.Rnw:401-402
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 34: Workflow.Rnw:414-429
###################################################

## find peaks for GFP control
cov.gfp <- coverage(resize(cstest$gfp, 200))
peaks.gfp <- slice(cov.gfp, lower = 8)

peakSummary <- diffPeakSummary(peaks.gfp, peaks.ctcf)

xyplot(asinh(sums2) ~ asinh(sums1) | space,
       data = as.data.frame(peakSummary), 
       panel = function(x, y, ...) {
           panel.xyplot(x, y, ...)
           panel.abline(median(y - x), 1)
       },
       type = c("p", "g"), alpha = 0.5, aspect = "iso")



###################################################
### code chunk number 35: Workflow.Rnw:431-432
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 36: Workflow.Rnw:437-446
###################################################
peakSummary <- 
    within(peakSummary,
       {
           diffs <- asinh(sums2) - asinh(sums1)
           resids <- (diffs - median(diffs)) / mad(diffs)
           up <- resids > 2
           down <- resids < -2
           change <- ifelse(up, "up", ifelse(down, "down", "flat"))
       })


###################################################
### code chunk number 37: Workflow.Rnw:462-465
###################################################
library(TxDb.Mmusculus.UCSC.mm9.knownGene)
gregions <- transcripts(TxDb.Mmusculus.UCSC.mm9.knownGene)
gregions


###################################################
### code chunk number 38: Workflow.Rnw:469-470
###################################################
promoters <- flank(gregions, 1000, both = TRUE)


###################################################
### code chunk number 39: Workflow.Rnw:474-476
###################################################
peakSummary$inPromoter <- peakSummary %in% promoters
xtabs(~ inPromoter + change, peakSummary)


###################################################
### code chunk number 40: Workflow.Rnw:479-481
###################################################
peakSummary$inUpstream <- peakSummary %in% flank(gregions, 20000)
peakSummary$inGene <- peakSummary %in% gregions


###################################################
### code chunk number 41: Workflow.Rnw:484-492
###################################################
sumtab <- 
    as.data.frame(rbind(total = xtabs(~ change, peakSummary),
                        promoter = xtabs(~ change, 
                          subset(peakSummary, inPromoter)),
                        upstream = xtabs(~ change, 
                          subset(peakSummary, inUpstream)),
                        gene = xtabs(~ change, subset(peakSummary, inGene))))
##cbind(sumtab, ratio = round(sumtab[, "down"] /  sumtab[, "up"], 3))


###################################################
### code chunk number 42: rtracklayer-upload (eval = FALSE)
###################################################
## library(rtracklayer)
## session <- browserSession()
## genome(session) <- "mm9"
## session$gfpCov <- cov.gfp
## session$gfpPeaks <- peaks.gfp
## session$ctcfCov <- cov.ctcf
## session$ctcfPeaks <- peaks.ctcf


###################################################
### code chunk number 43: rtracklayer-view (eval = FALSE)
###################################################
## peak.ord <- order(unlist(peak.depths), decreasing=TRUE)
## peak.sort <- as(peaks.ctcf, "GRanges")[peak.ord]
## view <- browserView(session, peak.sort[1], full = c("gfpCov", "ctcfCov"))


###################################################
### code chunk number 44: rtracklayer-view-5 (eval = FALSE)
###################################################
## views <- browserView(session, head(peak.sort, 5), full = c("gfpCov", "ctcfCov"))


###################################################
### code chunk number 45: Workflow.Rnw:534-535
###################################################
sessionInfo()