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

###################################################
### code chunk number 1: prepare (eval = FALSE)
###################################################
## options(length=55, digits=3)
## options(SweaveHooks=list(
##    along=function() par(mar=c(2.5,4.2,4,1.5), font.lab=2),
##    boxplot=function() par(mar=c(5,5,1,1), font.lab=2)))
## set.seed(1)


###################################################
### code chunk number 2: loadpackage (eval = FALSE)
###################################################
## library("Ringo")
## library("biomaRt")
## library("topGO")
## library("xtable")
## library("ccTutorial")
## library("lattice")


###################################################
### code chunk number 3: locateData (eval = FALSE)
###################################################
## pairDir <- system.file("PairData",package="ccTutorial") 
## list.files(pairDir, pattern="pair$")


###################################################
### code chunk number 4: remark1 (eval = FALSE)
###################################################
## # the following chunk 'readNimblegen' requires at least 1GB of RAM
## # and takes about 10 minutes. If time and memory are issues, you can
## # skip this step, see chunk 'remark2' below.


###################################################
### code chunk number 5: readNimblegen (eval = FALSE)
###################################################
## RGs <- lapply(sprintf("files_array%d.txt",1:4),
##   readNimblegen, "spottypes.txt", path=pairDir)


###################################################
### code chunk number 6: loadProbeAnno (eval = FALSE)
###################################################
## data("probeAnno")
## allChrs <- chromosomeNames(probeAnno)


###################################################
### code chunk number 7: showmm9genes (eval = FALSE)
###################################################
## data("mm9genes")
## mm9genes[sample(nrow(mm9genes), 4), 
##    c("name", "chr", "strand", "start", "end", "symbol")]


###################################################
### code chunk number 8: loadMm9.gene2GO (eval = FALSE)
###################################################
## data("mm9.gene2GO")


###################################################
### code chunk number 9: loadGenesGOAnnotation (eval = FALSE)
###################################################
## data("mm9.g2p")


###################################################
### code chunk number 10: arrayGenes (eval = FALSE)
###################################################
## arrayGenes <- names(mm9.g2p)[listLen(mm9.g2p)>=5]
## arrayGenesWithGO <- intersect(arrayGenes, names(mm9.gene2GO))


###################################################
### code chunk number 11: remark2 (eval = FALSE)
###################################################
## # the following chunk 'preprocess' requires at least 1GB of RAM
## # and takes about 5 minutes. If time and memory are issues, 
## # instead of running that chunk, you can load the result 'X', the
## # ExpressionSet holding the fold changes after preprocessing, by
## data("X")


###################################################
### code chunk number 12: preprocess (eval = FALSE)
###################################################
## MAs <- lapply(RGs, function(thisRG)
##   preprocess(thisRG[thisRG$genes$Status=="Probe",], 
##              method="nimblegen", returnMAList=TRUE))
## MA <- do.call(rbind, MAs)
## X  <- asExprSet(MA)
## sampleNames(X) <- paste(X$Cy5, X$Tissue, sep=".")


###################################################
### code chunk number 13: chipAlongChromActc1 (eval = FALSE)
###################################################
## plot(X, probeAnno, chrom="2", xlim=c(113.8725e6,113.8835e6), 
##      ylim=c(-3,5), gff=mm9genes, paletteName='Set2')


###################################################
### code chunk number 14: smoothing (eval = FALSE)
###################################################
## smoothX <- computeRunningMedians(X, probeAnno=probeAnno, 
##   modColumn="Tissue", allChr=allChrs, winHalfSize=450, min.probes=5)
## sampleNames(smoothX) <- paste(sampleNames(X),"smoothed",sep=".")
## combX <- cbind2(X, smoothX)


###################################################
### code chunk number 15: smoothAlongChromActc1 (eval = FALSE)
###################################################
## plot(combX, probeAnno, chrom="2", xlim=c(113.8725e6,113.8835e6),
##      ylim=c(-3,5), gff=mm9genes, 
##      colPal=c(brewer.pal(8,"Set2")[1:2],brewer.pal(8,"Dark2")[1:2]))


###################################################
### code chunk number 16: computeY0 (eval = FALSE)
###################################################
## y0 <- apply(exprs(smoothX), 2, upperBoundNull, prob=0.99)


###################################################
### code chunk number 17: histogramsSmoothed (eval = FALSE)
###################################################
## myPanelHistogram <- function(x, ...){
##   panel.histogram(x, col=brewer.pal(8,"Dark2")[panel.number()], ...)
##   panel.abline(v=y0[panel.number()], col="red")
## }
## 
## h = histogram( ~ y | z, 
##       data = data.frame(
##         y = as.vector(exprs(smoothX)), 
##         z = rep(X$Tissue, each = nrow(smoothX))), 
##       layout = c(1,2), nint = 50, 
##       xlab = "smoothed reporter level [log2]",
##       panel = myPanelHistogram)
## 
## print(h)


###################################################
### code chunk number 18: computeY0Echo (eval = FALSE)
###################################################
## y0 <- apply(exprs(smoothX), 2, upperBoundNull, prob=0.99)


###################################################
### code chunk number 19: cherFinding (eval = FALSE)
###################################################
## chersX <- findChersOnSmoothed(smoothX, 
##    probeAnno = probeAnno, 
##    thresholds = y0, 
##    allChr = allChrs, 
##    distCutOff = 450, 
##    minProbesInRow = 5, 
##    cellType = X$Tissue)


###################################################
### code chunk number 20: relateChers (eval = FALSE)
###################################################
## chersX <- relateChers(chersX, mm9genes, upstream=5000)


###################################################
### code chunk number 21: loadCherFinding (eval = FALSE)
###################################################
## # since especially the call to relateChers takes some time, we load the
## ## pre-saved image here:
## data("chersX")


###################################################
### code chunk number 22: showChers (eval = FALSE)
###################################################
## chersXD <- as.data.frame(chersX)
## head(chersXD[
##   order(chersXD$maxLevel, decreasing=TRUE), 
##   c("chr", "start", "end", "cellType", "features", "maxLevel", "score")])


###################################################
### code chunk number 23: plotCher1 (eval = FALSE)
###################################################
## plot(chersX[[which.max(chersXD$maxLevel)]], smoothX, probeAnno=probeAnno, 
##      gff=mm9genes, paletteName="Dark2", ylim=c(-1,6))


###################################################
### code chunk number 24: showCellType (eval = FALSE)
###################################################
## table(chersXD$cellType)


###################################################
### code chunk number 25: getGenesEnrichedPerTissue (eval = FALSE)
###################################################
## brainGenes <- getFeats(chersX[sapply(chersX, cellType)=="brain"])
## heartGenes <- getFeats(chersX[sapply(chersX, cellType)=="heart"])
## brainOnlyGenes <- setdiff(brainGenes, heartGenes)
## heartOnlyGenes <- setdiff(heartGenes, brainGenes)


###################################################
### code chunk number 26: useTopGO (eval = FALSE)
###################################################
## brainRes <- sigGOTable(brainOnlyGenes, gene2GO=mm9.gene2GO,
##                        universeGenes=arrayGenesWithGO)
## print(brainRes)


###################################################
### code chunk number 27: useTopGOHeart (eval = FALSE)
###################################################
## heartRes <- sigGOTable(heartOnlyGenes,  gene2GO=mm9.gene2GO,
##                        universeGenes=arrayGenesWithGO)
## print(heartRes)


###################################################
### code chunk number 28: loadExpressionData (eval = FALSE)
###################################################
## data("barreraExpressionX")


###################################################
### code chunk number 29: loadArrayGenesToProbeSets (eval = FALSE)
###################################################
## data("arrayGenesToProbeSets")


###################################################
### code chunk number 30: compareChIPAndExpression (eval = FALSE)
###################################################
## bX <- exprs(barreraExpressionX)
## allH3K4me3Genes  <- union(brainGenes, heartGenes)
## allH3K4ProbeSets <- unlist(arrayGenesToProbeSets[allH3K4me3Genes])
## noH3K4ProbeSets  <- setdiff(rownames(bX), allH3K4ProbeSets)
## brainH3K4ExclProbeSets <- unlist(arrayGenesToProbeSets[brainOnlyGenes])
## heartH3K4ExclProbeSets <- unlist(arrayGenesToProbeSets[heartOnlyGenes])
## 
## brainIdx <- barreraExpressionX$Tissue=="Brain"
## 
## brainExpression <- list(
##   H3K4me3BrainNoHeartNo  = bX[noH3K4ProbeSets, brainIdx],
##   H3K4me3BrainYes        = bX[allH3K4ProbeSets, brainIdx],
##   H3K4me3BrainYesHeartNo = bX[brainH3K4ExclProbeSets, brainIdx],
##   H3K4me3BrainNoHeartYes = bX[heartH3K4ExclProbeSets, brainIdx]
## )


###################################################
### code chunk number 31: H3K4me3VsExpression (eval = FALSE)
###################################################
## boxplot(brainExpression, col=c("#666666","#999966","#669966","#996666"), 
##         names=NA, varwidth=TRUE, log="y", 
##         ylab='gene expression level in brain cells')
## mtext(side=1, at=1:length(brainExpression), padj=1, font=2, 
##       text=rep("H3K4me3",4), line=1)
## mtext(side=1, at=c(0.2, 1:length(brainExpression)), padj=1, font=2, 
##       text=c("brain/heart","-/-","+/+","+/-","-/+"), line=2)


###################################################
### code chunk number 32: testExpressionGreater (eval = FALSE)
###################################################
## with(brainExpression, 
##      wilcox.test(H3K4me3BrainYesHeartNo, H3K4me3BrainNoHeartNo, 
##                  alternative="greater"))


###################################################
### code chunk number 33: sessionInfo (eval = FALSE)
###################################################
## toLatex(sessionInfo())


###################################################
### code chunk number 34: remarkTables (eval = FALSE)
###################################################
## # the purpose of the following chunks is merely to provide pretty
## # latex-formated output of the tables generated in the tutorial


###################################################
### code chunk number 35: printMm9Genes (eval = FALSE)
###################################################
## print(xtable(mm9genes[sample(nrow(mm9genes), 4), 
##    c("name", "chr", "strand", "start", "end", "symbol")],
##    label="tab-mm9genes",
##    caption="\\sl An excerpt of the object 'mm9genes'."),
##    type="latex", table.placement="h!t", size="scriptsize",
##    include.rownames=FALSE)


###################################################
### code chunk number 36: printChersXD (eval = FALSE)
###################################################
## print(xtable(head(chersXD[order(chersXD$maxLevel, decreasing=TRUE), 
##   c("chr", "start", "end", "cellType", "features", "maxLevel", "score")]),
##   label="tab-chersXD",
##   caption="\\sl The first six lines of object 'chersXD'."),
##   type="latex", table.placement="h!t", size="scriptsize",
##   include.rownames=FALSE)


###################################################
### code chunk number 37: printBrainRes (eval = FALSE)
###################################################
## ## for having prettier tables in the PDF, we use 'xtable' here:
## print(xtable(brainRes, label="tab-brainResGO", 
##    caption="\\sl GO terms that are significantly over-represented among genes showing H3K4me3 enrichment specifically in brain cells"),
##    type="latex", table.placement="h!t", size="scriptsize",
##    include.rownames=FALSE)


###################################################
### code chunk number 38: printHeartRes (eval = FALSE)
###################################################
## print(xtable(heartRes, label="tab-heartResGO",
##    caption="\\sl GO terms that are significantly over-represented among genes showing H3K4me3 enrichment specifically in heart cells"),
##    type="latex", table.placement="h!b", size="scriptsize",
##    include.rownames=FALSE)