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

###################################################
### code chunk number 1: installPackages (eval = FALSE)
###################################################
## source("http://www.bioconductor.org/biocLite.R")
## biocLite(c("Ringo", "biomaRt", "topGO", "ccTutorial"))


###################################################
### code chunk number 2: prepare (eval = FALSE)
###################################################
## options(length=60, "stringsAsFactors"=FALSE)
## set.seed(123)
## 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=4),
##    dens=function() par(mar=c(4.1, 4.1, 0.1, 0.1), font.lab=2)))


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


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


###################################################
### code chunk number 5: exampleFilesTxt (eval = FALSE)
###################################################
## read.delim(file.path(pairDir,"files_array1.txt"), header=TRUE)


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


###################################################
### code chunk number 7: showRG (eval = FALSE)
###################################################
## head(RGs[[1]]$R)
## head(RGs[[1]]$G)
## tail(RGs[[1]]$genes)


###################################################
### code chunk number 8: showProbeStatus (eval = FALSE)
###################################################
## table(RGs[[1]]$genes$Status)


###################################################
### code chunk number 9: imageRG (eval = FALSE)
###################################################
## RG1breaks <- c(0,quantile(RGs[[1]]$G, probs=seq(0,1,by=0.1)),2^16)
## png("ccTutorialArrayImages.png", units="in", res=200,
##      height=10.74*1.5, width=7.68*1.5)
## par(mar=c(0.01,0.01,2.2,0.01))
## layout(matrix(c(1,2,5,6,3,4,7,8,9,10,13,14,11,12,15,16),
##        ncol=4,byrow=TRUE))
## for (this.set in 1:4){
##   thisRG <- RGs[[this.set]]
##   for (this.channel in c("green","red")){
##     my.colors <- colorRampPalette(c("black",paste(this.channel,c(4,1),
##                                     sep="")))(length(RG1breaks)-1)
##     for (arrayno in 1:2){
##       image(thisRG, arrayno, channel=this.channel, 
##             mybreaks=RG1breaks, mycols=my.colors)
##       mtext(side=3, line=0.2, font=2, text=colnames(thisRG[[toupper(
##             substr(this.channel,1,1))]])[arrayno])
## }}}
## dev.off()


###################################################
### code chunk number 10: corPlotRG2G (eval = FALSE)
###################################################
## corPlot(log2(RGs[[2]]$G))


###################################################
### code chunk number 11: corPlotRG2R (eval = FALSE)
###################################################
## corPlot(log2(RGs[[2]]$R))


###################################################
### code chunk number 12: posToProbeAnno (eval = FALSE)
###################################################
## probeAnno <- posToProbeAnno(file.path(system.file("exonerateData",
##   package="ccTutorial"), "allChromExonerateOut.txt"))
## allChrs <- chromosomeNames(probeAnno)


###################################################
### code chunk number 13: posToProbeAnnoExtent (eval = FALSE)
###################################################
## genome(probeAnno) <- "M. musculus (mm9)"
## arrayName(probeAnno) <- "2005-06-17_Ren_MM5Tiling"


###################################################
### code chunk number 14: showProbeAnno (eval = FALSE)
###################################################
## show(probeAnno)
## ls(probeAnno)


###################################################
### code chunk number 15: showProbeAnno2 (eval = FALSE)
###################################################
## table(probeAnno["9.unique"])


###################################################
### code chunk number 16: reporterSpacing (eval = FALSE)
###################################################
## startDiffByChr <- lapply(as.list(allChrs), function(chr){
##   chrsta <- probeAnno[paste(chr,"start",sep=".")]
##   chruni <- probeAnno[paste(chr,"unique",sep=".")]
##   ## get start positions of unique reporter match positions
##   return(diff(sort(chrsta[chruni=="0"])))})
## startDiff <- unlist(startDiffByChr, use.names=FALSE)
## table(cut(startDiff, breaks=c(0,50,99,100,200,1000,max(startDiff))))


###################################################
### code chunk number 17: makeGffWithBiomaRt (eval = FALSE)
###################################################
## ensembl  <- useMart("ensembl", dataset="mmusculus_gene_ensembl")
## gene.ids <- unique(unlist(lapply(as.list(c(1:19,"X","Y")), 
##    function(this.chr) 
##     getBM(attributes="ensembl_gene_id", filters="chromosome_name", 
##           values=this.chr, mart=ensembl)[,1]), use.names=FALSE))
##     sel.attributes=c("ensembl_gene_id", "mgi_symbol", "chromosome_name",
##           "strand", "start_position","end_position", "description")
## mm9genes <- getBM(attributes=sel.attributes, filters="ensembl_gene_id",
##                   value=gene.ids, mart=ensembl)


###################################################
### code chunk number 18: replaceMm9GenesNames (eval = FALSE)
###################################################
## mm9genes$name    <- mm9genes$"ensembl_gene_id"
## mm9genes$gene    <- mm9genes$"ensembl_gene_id"
## mm9genes$chr     <- mm9genes$"chromosome_name"
## mm9genes$symbol  <- mm9genes$"mgi_symbol"
## mm9genes$start   <- mm9genes$"start_position"
## mm9genes$end     <- mm9genes$"end_position"
## mm9genes$feature <- rep("gene",nrow(mm9genes))


###################################################
### code chunk number 19: removeDuplicatedSymbols (eval = FALSE)
###################################################
## if (any(duplicated(mm9genes$name))){
##   dupl <- unique(mm9genes$name[duplicated(mm9genes$name)])
##   G <- lapply(as.list(dupl), function(this.gene){
##     this.gff <- subset(mm9genes,name == this.gene)
##     if (nrow(unique(this.gff[,c("name","chr","start","end",
##         "description")]))>1) return(this.gff[1,,drop=FALSE])
##     non.zero.gff <- subset(this.gff, nchar(symbol)>0)
##     this.other.sym <- NULL
##     if (nrow(non.zero.gff)> 0){
##       shortest <- which.min(nchar(non.zero.gff$symbol))
##       this.new.sym <- non.zero.gff$symbol[shortest]
##       if (nrow(non.zero.gff)>1)
##         this.other.sym <- paste("Synonyms", 
##            paste(non.zero.gff$symbol[-shortest],collapse=","),sep=":")
##     } else { this.new.sym <- "" }
##     this.gff$symbol[1] <- this.new.sym
##     if (!is.null(this.other.sym))
##       this.gff$description[1] <- paste(this.gff$description[1],
##                                        this.other.sym,sep=";")
##     return(this.gff[1,,drop=FALSE])
##   })
##   mm9genes <- rbind(mm9genes[-which(mm9genes$name %in% dupl),],
##                     do.call("rbind",G))
## }


###################################################
### code chunk number 20: reorderMm9 (eval = FALSE)
###################################################
## mm9genes <- mm9genes[order(mm9genes$chr, mm9genes$start),
##   c("name","chr","strand","start","end","symbol","description","feature")]
## rownames(mm9genes) <- NULL


###################################################
### code chunk number 21: loadMM9Genes (eval = FALSE)
###################################################
## data(mm9genes)


###################################################
### code chunk number 22: setSeed1 (eval = FALSE)
###################################################
## set.seed(1)


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


###################################################
### code chunk number 24: getGenesGOAnnotation (eval = FALSE)
###################################################
## ensembl  <- useMart("ensembl", dataset="mmusculus_gene_ensembl")
## ontoGOs <- lapply(as.list(c("biological_process","cellular_component", 
##                   "molecular_function")), function(onto){
##   ontoBM <- getBM(mart=ensembl, attributes=c("ensembl_gene_id", 
##                    paste("go",onto,"id", sep="_"), 
##                    paste("go",onto,"linkage_type", sep="_")), 
##                   filters="ensembl_gene_id", value=mm9genes$name)
##   names(ontoBM) <- c("ensembl_gene_id","go","evidence_code")
##   ontoBM <- subset(ontoBM,!( evidence_code %in% c("","IEA","NAS","ND")))
## })
## mm9GO <- do.call("rbind", ontoGOs)


###################################################
### code chunk number 25: getGOperGene (eval = FALSE)
###################################################
## mm9.gene2GO <- with(mm9GO, split(go, ensembl_gene_id))


###################################################
### code chunk number 26: loadGenesGOAnnotation (eval = FALSE)
###################################################
## data(mm9.gene2GO)
## data(mm9.g2p)


###################################################
### code chunk number 27: mappingGenesToProbes (eval = FALSE)
###################################################
## mm9.g2p <- features2Probes(gff=mm9genes, probeAnno=probeAnno)


###################################################
### code chunk number 28: showMm9g2p (eval = FALSE)
###################################################
## table(cut(listLen(mm9.g2p),breaks=c(-1,0,10,50,100,500,1200)))


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


###################################################
### code chunk number 30: 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 31: showX (eval = FALSE)
###################################################
## show(X)


###################################################
### code chunk number 32: chipAlongChromCrmp1 (eval = FALSE)
###################################################
## plot(X, probeAnno, chrom="5", xlim=c(37.63e6,37.64e6), ylim=c(-3,5),
##      gff=mm9genes, paletteName="Set2")


###################################################
### code chunk number 33: 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 34: smoothAlongChromCrmp1 (eval = FALSE)
###################################################
## plot(combX, probeAnno, chrom="5", xlim=c(37.63e6,37.64e6),
##      gff=mm9genes, ylim=c(-3,5),
##      colPal=c(brewer.pal(8,"Set2")[1:2],brewer.pal(8,"Dark2")[1:2]))


###################################################
### code chunk number 35: getChersXD (eval = FALSE)
###################################################
## data("chersX")
## chersXD <- as.data.frame(chersX)


###################################################
### code chunk number 36: tryCmarrt (eval = FALSE)
###################################################
## cmarrtDat <- do.call("rbind", lapply(as.list(allChrs), function(chr){
##   areUni <- probeAnno[paste(chr,"unique",sep=".")]==0
##   chrIdx <- match(probeAnno[paste(chr,"index",sep=".")][areUni],
##                   featureNames(X))
##   chrDat <- data.frame("chr"=rep(chr, sum(areUni)),
##      "start"=probeAnno[paste(chr,"start",sep=".")][areUni],
##      "stop"=probeAnno[paste(chr,"end",sep=".")][areUni],
##      "logR"=exprs(X)[chrIdx,1],
##      stringsAsFactors=FALSE)
## }))
## 
## cmarrtRes <- cmarrt.ma(cmarrtDat, M=0.5, frag.length=900,
##                        window.opt = "fixed.gen.dist")
## 
## cmarrtReg <- cmarrt.peak(cmarrtRes, alpha=0.05, method="BY", minrun=4) 
## cmarrtRegDf <- lapply(cmarrtReg, as.data.frame)$cmarrt.bound
## names(cmarrtRegDf)[1:3] <- c("chr","start","end")


###################################################
### code chunk number 37: showCmarrt (eval = FALSE)
###################################################
## head(cmarrtRegDf)


###################################################
### code chunk number 38: getRingoAndCmarrt (eval = FALSE)
###################################################
## ringoChersChr9 <- subset(chersXD, chr=="9" & cellType=="brain")
## cmarrtChersChr9 <- subset(cmarrtRegDf, chr=="9")
## dim(ringoChersChr9)
## dim(cmarrtChersChr9)


###################################################
### code chunk number 39: overlapRingoVsCmarrt (eval = FALSE)
###################################################
## chersChr9Overlap <- as.matrix(
##    regionOverlap(ringoChersChr9, cmarrtChersChr9))
## minRegChr9Len <- outer(with(ringoChersChr9, end-start+1), 
##                        with(cmarrtChersChr9, end-start+1), pmin)
## fracChr9Overlap <- chersChr9Overlap /minRegChr9Len


###################################################
### code chunk number 40: showOverlapRingoVsCmarrt (eval = FALSE)
###################################################
## summary(apply(fracChr9Overlap, 1, max))


###################################################
### code chunk number 41: loadCherFinding (eval = FALSE)
###################################################
## data(chersX)
## chersXD <- as.data.frame(chersX)


###################################################
### code chunk number 42: 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 43: useTopGO (eval = FALSE)
###################################################
## brainRes <- sigGOTable(brainOnlyGenes, gene2GO=mm9.gene2GO,
##                        universeGenes=arrayGenesWithGO)
## heartRes <- sigGOTable(heartOnlyGenes,  gene2GO=mm9.gene2GO,
##                        universeGenes=arrayGenesWithGO)


###################################################
### code chunk number 44: computeRegionsOverlap (eval = FALSE)
###################################################
## brainRegions <- subset(chersXD, cellType=="brain")
## heartRegions <- subset(chersXD, cellType=="heart")
## chersOBL <- as.matrix(regionOverlap(brainRegions, heartRegions))
## minRegLen <- outer(with(brainRegions, end-start+1),
##                    with(heartRegions, end-start+1), pmin)
## fracOverlap <- chersOBL/minRegLen


###################################################
### code chunk number 45: compTissueSpecificRegions (eval = FALSE)
###################################################
## brainSpecReg <-  brainRegions[rowMax(fracOverlap)<0.7,]
## heartSpecReg <-  heartRegions[rowMax(t(fracOverlap))<0.7,]
## mean(is.element(unlist(strsplit(brainSpecReg$features, 
##    split="[[:space:]]"), use.names=FALSE), brainOnlyGenes))
## selGenes <- intersect(unlist(strsplit(brainSpecReg$features, 
##    split="[[:space:]]"), use.names=FALSE), heartGenes)


###################################################
### code chunk number 46: targetPos (eval = FALSE)
###################################################
## targetPos <- seq(-5000, 10000, by=250)


###################################################
### code chunk number 47: newQuantsOverPostions (eval = FALSE)
###################################################
## selQop <- quantilesOverPositions(smoothX,
##    selGenes=selGenes, quantiles=c(0.5, 0.9),
##    g2p=mm9.g2p, positions=targetPos)


###################################################
### code chunk number 48: sepRegGenesSmoothedQuantiles (eval = FALSE)
###################################################
## plot(selQop, c("green","orange"))


###################################################
### code chunk number 49: topGOsepRegGenes (eval = FALSE)
###################################################
## sepRegRes <- sigGOTable(selGenes=selGenes, gene2GO=mm9.gene2GO,
##                         universeGenes=arrayGenesWithGO)
## print(sepRegRes)


###################################################
### code chunk number 50: printSepRegRes (eval = FALSE)
###################################################
## # this chunk only provides a prettier output of the table sepRegRes
## #  in latex format
## print(xtable(sepRegRes, label="tab-sepRegResGO", 
##    caption="\\sl GO terms that are significantly over-represented among genes that show different H3K4me3 regions in heart and brain cells"),
##    type="latex", table.placement="htb", size="scriptsize",
##    include.rownames=FALSE)


###################################################
### code chunk number 51: processBarerraExpressionData (eval = FALSE)
###################################################
## library("affy")
## library("mouse4302cdf")
## AB <- ReadAffy(celfile.path=system.file("expression", 
##                                         package="ccTutorial"))
## barreraExpressionX <- mas5(AB)
## barreraExpressionX$Tissue <- sapply(
##    strsplit(sampleNames(barreraExpressionX),split="\\."),"[",3)


###################################################
### code chunk number 52: mapEnsToAffy (eval = FALSE)
###################################################
## ensembl <- useMart("ensembl", dataset="mmusculus_gene_ensembl")
## bmRes <- getBM(attributes=c("ensembl_gene_id","affy_mouse430_2"), 
##                filters="ensembl_gene_id", value=arrayGenes, 
##                mart=ensembl)
## bmRes <- subset(bmRes, nchar(affy_mouse430_2)>0)
## arrayGenesToProbeSets <- split(bmRes[["affy_mouse430_2"]], 
##                                bmRes[["ensembl_gene_id"]])


###################################################
### code chunk number 53: lookArrayGenesToProbeSets (eval = FALSE)
###################################################
## table(listLen(arrayGenesToProbeSets))


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