### R code from vignette source 'chipseq_thurs_pm.Rnw'

###################################################
### code chunk number 1: adjust_prompt_1
###################################################
options(prompt=" ", continue=" ")


###################################################
### code chunk number 2: load_libraries
###################################################
library("GenomicRanges")
library("Repitools")
library("BSgenome.Hsapiens.UCSC.hg18")
library("chipseq")
library("matrixStats")


###################################################
### code chunk number 3: load_chipseq_data
###################################################
load("d.Rdata")


###################################################
### code chunk number 4: load_mbdseq_data
###################################################
load("methGR.Rdata")


###################################################
### code chunk number 5: explore_chip_seq_data
###################################################
class(d)
class(d[[1]])
names(d)
d
table( seqnames(d[[1]])  )
table( width(d)[[1]] )
elementLengths(d)
lapply( d, function(u) table( width(u) ))


###################################################
### code chunk number 6: load_anno_affy_expression
###################################################
load("anno.Rdata")
load("expr.Rdata")


###################################################
### code chunk number 7: enrichment_plot
###################################################
cols <- c("black","red","blue","darkgreen","orange","grey")
e <- enrichmentPlot(d, seq.len=300, lwd=4, xlim=c(0,200), col=cols)


###################################################
### code chunk number 8: enrichment_plot
###################################################
lapply(e[1:2],head)


###################################################
### code chunk number 9: variation_on_enrichment_plot
###################################################
# set ranges to plot
xr <- range( unlist(lapply(e,".subset",1))+1 )
yr <- quantile(unlist(lapply(e,".subset",2)),p=c(.02,.98))

# Plot coverage frequency by ChIP experiment, manually
plot(xr*10,type="n",xlim=xr, ylim=yr,
     xlab="Coverage+1", ylab="Frequency", log="xy")
dummy <- sapply(1:length(e), function(u) {
  points( e[[u]]$coverage, e[[u]]$bases, col=cols[u], lwd=4, type="b")
})
legend("topright",names(d),col=cols,lwd=4)


###################################################
### code chunk number 10: cpg_density_plot
###################################################
cpgDensityPlot(methGR, seq.len=300, organism=Hsapiens, 
               lwd=4, w.function="linear", window=600, verbose=TRUE)


###################################################
### code chunk number 11: find_regions_simple
###################################################
findRegionsSimple <- function(u,frag.len=300,lower=15) {
  # u is a 'GRanges' of reads
  u <- resize(u,frag.len)
  cv <- coverage(u)
  scv <- slice(cv, lower=lower)
  gr <- as(scv,"GRanges")
  mcols(gr)$view <- mcols(gr)$view + 1e6*as.numeric(seqnames(gr))
  gr <- unlist(reduce(split(gr,mcols(gr)$view)))
  mcols(gr)$mean <- unlist(viewMeans(scv))
  mcols(gr)$max <- unlist(viewMaxs(scv))
  gr
}


###################################################
### code chunk number 12: estimate_frag_size
###################################################
load("d.Rdata")
fr <- sapply(d, estimate.mean.fraglen)
fr
frm <- colMedians(fr)
frm


###################################################
### code chunk number 13: region_call
###################################################
regs <- mapply( function(u,v) {
  cat(".")
  findRegionsSimple(u,frag.len=round(v))
}, d, frm)


###################################################
### code chunk number 14: extract_and_get_limits
###################################################
whs <- sapply(regs, function(u) {
  cbind(width=width(u),mean=mcols(u)$mean)
})
xlim <- quantile( unlist(lapply(whs,function(u) u[,1])), p=c(.01,.99) )
ylim <- quantile( unlist(lapply(whs,function(u) u[,2])), p=c(.01,.99) )


###################################################
### code chunk number 15: peak_widths_heights
###################################################
par(mfrow=c(3,2))
invisible(mapply( function(u,v) {
  plot(u,main=v,pch=19,cex=.4,log="xy", xlim=xlim, ylim=ylim); grid();
}, whs, names(whs)))


###################################################
### code chunk number 16: gather_feature_scores
###################################################
covs <- featureScores(d, anno, up=10000, down=10000, freq=100, s.width=500)


###################################################
### code chunk number 17: bin_plots_heatmap
###################################################
binPlots(covs[3],ordering=data.frame(affy=expr),
         ord.label="Expression",n.bins=50,plot.type="heatmap")


###################################################
### code chunk number 18: bin_plots_line
###################################################
binPlots(covs[3],ordering=data.frame(affy=expr),
         ord.label="Expression",n.bins=20,plot.type="line")


###################################################
### code chunk number 19: cluster_plots
###################################################
clusterPlots(covs[-c(4,6)], function(x) sqrt(x), expr=expr, plot.type="heatmap",
             t.name="ChIP-seq related to Affy expression", n.clusters=8)


###################################################
### code chunk number 20: annotation_blocks_counts
###################################################
ac <- annotationCounts(d, anno, up=500, down=500, seq.len=200)
head(ac)


###################################################
### code chunk number 21: sessInfo
###################################################
toLatex(sessionInfo())


###################################################
### code chunk number 22: adjust_prompt_1
###################################################
options(prompt="> ", continue="+ ")