### R code from vignette source 'vignettes/compEpiTools/inst/doc/compEpiTools.rnw'

###################################################
### code chunk number 1: options
###################################################
options(width=72)


###################################################
### code chunk number 2: libs
###################################################
library(compEpiTools)
require(BSgenome.Mmusculus.UCSC.mm9)
require(TxDb.Mmusculus.UCSC.mm9.knownGene)
require(org.Mm.eg.db)
require(rtracklayer)
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene


###################################################
### code chunk number 3: GRcounting
###################################################
bampath <- system.file("extdata", "ex1.bam", package="Rsamtools")
ir <- IRanges(start=c(1000, 100), end=c(2000, 1000))
gr <- GRanges(seqnames=Rle(c('seq1','seq2')), ranges=ir)
res <- GRbaseCoverage(Object=gr, bam=bampath)
GRcoverage(Object=gr, bam=bampath, Nnorm=FALSE, Snorm=FALSE)[1]
GRcoverageInbins(Object=gr, bam=bampath, Nnorm=FALSE, Snorm=FALSE, Nbins=5)[1,]
summit <- GRcoverageSummit(Object=gr, bam=bampath)
plot(res[[1]], type='l', xlab='bp', ylab='reads count')
abline(v=start(summit[1])-start(gr[1])+1, lty=2, lwd=2)


###################################################
### code chunk number 4: GRann
###################################################
TSSpos <- TSS(txdb)
gr <- TSSpos[1:5]
start(gr) <- start(gr) - 1000
end(gr) <- end(gr) - 600
mcols(gr) <- NULL
# retrieving CGI mm9 islands from UCSC annotation tables
session <- browserSession()
genome(session) <- 'mm9'
query <- ucscTableQuery(session, 'cpgIslandExt')
CGIgr <- as(track(query), 'GRanges')
res <- GRannotate(Object=GRmidpoint(gr), txdb=txdb, EG2GS=org.Mm.egSYMBOL, 
	upstream=2000, downstream=1000, userAnn=GRangesList(CGI=CGIgr))
show(res)


###################################################
### code chunk number 5: GRheat1
###################################################
gr <- TSSpos[1:50]
start(gr) <- start(gr) - 1000
end(gr) <- end(gr) - 600
extgr <- GRanges(seqnames(gr), ranges=IRanges(start(gr) - 1000, end(gr) + 1000))
data <- heatmapData(grl=list(ChIPseq=gr), refgr=extgr, type='gr', nbins=20, txdb=txdb)
pvalues <- c(runif(20,1e-20,1e-8), runif(15,1e-4,1e-2), runif(15,0.5,1))
pvalues <- cbind(pvalues, rep(0, 50), rep(0, 50))
rownames(data[[1]][[1]]) <- paste(1:50, signif(pvalues[,1],1), sep=' # ')
heatmapPlot(matList=data[[1]], clusterInds=1:3)


###################################################
### code chunk number 6: GRheat2
###################################################
heatmapPlot(matList=data[[1]], sigMat=pvalues, clusterInds=1:3)


###################################################
### code chunk number 7: info
###################################################
sessionInfo()