### R code from vignette source 'vignettes/ggbio/inst/doc/chip-seq.Rnw'

###################################################
### code chunk number 1: load
###################################################
library(chipseq)
library(GenomicFeatures)
data(cstest)
unique(seqnames(cstest))
## subset
chrs <- c("chr10", "chr11", "chr12")
cstest <- keepSeqlevels(cstest, chrs)
## estimate fragment length
fraglen <- estimate.mean.fraglen(cstest$ctcf)
fraglen[!is.na(fraglen)]
## that's around 200
## cstest.gr <- stack(cstest)
## head(cstest.gr)
## cstest.ext <- resize(cstest.gr, width = 200)
## extending them
ctcf.ext <- resize(cstest$ctcf, width = 200)
cov.ctcf <- coverage(ctcf.ext)
gfp.ext <- resize(cstest$gfp, width = 200)
cov.gfp <- coverage(gfp.ext)
## estimate peak cutoff
peakCutoff(cov.ctcf, fdr = 0.0001)
## we can use 8


###################################################
### code chunk number 2: find-close
###################################################
c.p <- cstest$ctcf[seqnames(cstest$ctcf) == "chr10" & 
                   strand(cstest$ctcf) == "+",]

c.n <- cstest$ctcf[seqnames(cstest$ctcf) == "chr10" & 
                   strand(cstest$ctcf) == "-",]


cov.p <- coverage(c.p)
cov.n <- coverage(c.n)
v1 <- viewWhichMaxs(slice(cov.p, lower = 8))$chr10
v2 <- viewWhichMaxs(slice(cov.n, lower = 8))$chr10
res <- expand.grid(v1, v2)
wh <- as.numeric(res[order(abs(res[,1] - res[, 2]))[1], ])
gr.wh <- GRanges("chr10", IRanges(wh[1], wh[2]))
gr.wh <- resize(gr.wh, width(gr.wh) + 200, fix = "center")



###################################################
### code chunk number 3: reads-close
###################################################
library(ggbio)
ctcf.sub <- subsetByOverlaps(cstest$ctcf, gr.wh)
p1 <- autoplot(ctcf.sub, aes(fill = strand), facets = strand ~ .)
ctcf.ext.sub <- subsetByOverlaps(ctcf.ext, gr.wh)
p2 <- autoplot(ctcf.ext.sub, aes(fill = strand), facets = strand ~ .)
tracks("original" = p1, "extending" = p2, heights = c(3, 5))


###################################################
### code chunk number 4: coverage-close
###################################################
ctcf.sub <- subsetByOverlaps(cstest$ctcf, gr.wh)
p1 <- autoplot(ctcf.sub, aes(fill = strand), facets = strand ~ ., stat = "coverage", 
               geom = "area")
ctcf.ext.sub <- subsetByOverlaps(ctcf.ext, gr.wh)
p2 <- autoplot(ctcf.ext.sub, aes(fill = strand), facets = strand ~ ., 
               stat = "coverage", geom = "area")
tracks("original" = p1, "extending" = p2)


###################################################
### code chunk number 5: genome-bin-no-ylim
###################################################
library(ggbio)
p1 <- autoplot(cov.ctcf)
p2 <- autoplot(cov.gfp)
tracks(ctcf = p1, gfp = p2)


###################################################
### code chunk number 6: genome-bin-ylim
###################################################
library(ggbio)
p1 <- autoplot(cov.ctcf)
p2 <- autoplot(cov.gfp)
tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 2e6))


###################################################
### code chunk number 7: genome-bin-maxs
###################################################
p1 <- autoplot(cov.ctcf, type = "viewMaxs")
p2 <- autoplot(cov.gfp,type = "viewMaxs")
tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 2e6))


###################################################
### code chunk number 8: genome-bin-100
###################################################
p1 <- autoplot(cov.ctcf, type = "viewMaxs", nbin = 100)
p2 <- autoplot(cov.gfp,type = "viewMaxs", nbin = 100)
tracks("ctcf" = p1, "gfp" = p2)  + coord_cartesian(ylim = c(0, 1e6))


###################################################
### code chunk number 9: genome-bin-heat
###################################################
p1 <- autoplot(cov.ctcf, type = "viewMeans", nbin = 100, geom = "heatmap")
p2 <- autoplot(cov.gfp,type = "viewMeans", nbin = 100, geom = "heatmap")
tracks(ctcf = p1, gfp = p2)  + scale_fill_continuous(limits = c(0, 8e05)) + 
                                 scale_color_continuous(limits = c(0, 8e05))


###################################################
### code chunk number 10: genome-slice
###################################################
p1 <- autoplot(cov.ctcf, type = "viewMaxs", stat = "slice", lower = 8)
p2 <- autoplot(cov.gfp,type = "viewMaxs", stat = "slice", lower = 8)
tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 15000))


###################################################
### code chunk number 11: genome-slice-drop
###################################################
p1 <- autoplot(cov.ctcf, type = "viewMaxs", stat = "slice", lower = 8)
p2 <- autoplot(cov.gfp,type = "viewMaxs", stat = "slice", lower = 8, drop = FALSE)
tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 15000))


###################################################
### code chunk number 12: genome-slice-rect
###################################################
p1 <- autoplot(cov.ctcf, stat = "slice", lower = 5, geom = "rect")
p2 <- autoplot(cov.gfp, stat = "slice", lower = 5, geom = "rect")
tracks(ctcf = p1, gfp = p2)


###################################################
### code chunk number 13: genome-slice-heatmap
###################################################
p1 <- autoplot(cov.ctcf, type = "viewMaxs", stat = "slice", lower = 8, geom = "heatmap")
p2 <- autoplot(cov.gfp,type = "viewMaxs", stat = "slice", lower = 8, geom = "heatmap", 
               drop = FALSE)
tracks("ctcf" = p1, "gfp" = p2) + scale_fill_continuous(limits = c(1000, 6000)) + 
                                 scale_color_continuous(limits = c(1000, 6000))


###################################################
### code chunk number 14: region
###################################################
peaks.ctcf <- slice(cov.ctcf, lower = 8)
peaks.gfp <- slice(cov.gfp, lower = 8)
## this function is from vignette of chipseq
peakSummary <- diffPeakSummary(peaks.gfp, peaks.ctcf)
 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"))
})
ps.down <- peakSummary[peakSummary$change == "down" & peakSummary$space == "chr11", ]
pk.down <- ps.down[order(ps.down$diffs),]
pk.down
## 
library(TxDb.Mmusculus.UCSC.mm9.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene
tx <- transcripts(txdb)
gn <- transcriptsBy(txdb, by = "gene")
fu <- fiveUTRsByTranscript(txdb)

idx <- which(countOverlaps(as(pk.down, "GRanges"), flank(fu, width = 100)) == 1)
wh.p <- as(pk.down[idx[2], ], "GRanges")
wh.pw <- resize(wh.p, width = 30000, fix = "center")



###################################################
### code chunk number 15: ideogram
###################################################
library(biovizBase)
mm9 <- getIdeogram("mm9")
cyto.def <- getOption("biovizBase")$cytobandColor
cyto.new <- c(cyto.def, c(gpos33 = "grey80", gpos66 = "grey60"))
p.ideo <- plotIdeogram(mm9, "chr10", zoom = c(start(wh.pw),end(wh.pw)))  + 
  scale_fill_manual(values = cyto.new) 
print(p.ideo)


###################################################
### code chunk number 16: txdb
###################################################
p.gene <- autoplot(txdb, which = wh.pw)


###################################################
### code chunk number 17: ideo-features
###################################################
## coverage
cstest.s <- stack(cstest)
cstest.s <- resize(cstest.s, width = 200)
cstest.sub <- subsetByOverlaps(cstest.s, wh.pw)
p.cov <- autoplot(cstest.sub, stat = "coverage", facets = sample ~ ., 
                  geom = "area")
## ideogram
tracks("ideogram" = p.ideo, "coverage" = p.cov, "gene" = p.gene, xlim = as(wh.pw, "GRanges"), 
       xlim.change = c(FALSE, TRUE, TRUE), 
       heights = c(1, 5, 5))