###################################################
### chunk number 1: sourceCodes
###################################################

source("calculateStatistics.R")
source("calculatePValues.R")



###################################################
### chunk number 2: packageInstallation eval=FALSE
###################################################
## 
##   source("http://www.bioconductor.org/biocLite.R")
##   biocLite("Biobase")
##   biocLite("macat")
##   biocLite("hgu95av2")
## 


###################################################
### chunk number 3: packageLoading
###################################################

library(Biobase)
library(macat)
library(hgu95av2)
loaddatapkg("stjudem")



###################################################
### chunk number 4: macatDemo eval=FALSE
###################################################
## demo(macatdemo)
## konqueror("Results6_T.html")


###################################################
### chunk number 5: preparation
###################################################

## 'inconsistency' of stjude (duplicated 'geneName's )
expr <- stjude$expr[stjude$geneName,]
selSix <- stjude$chromosome == "6"
## removes duplicated gene symbols
ids <- rownames(expr[selSix,])
chip <- stjude$chip
thisEnv = paste(chip, "SYMBOL", sep = "") 
symbol = mget(ids, env = eval(as.symbol(thisEnv)))
selSixFinal <- which(selSix)[!duplicated(symbol)]



###################################################
### chunk number 6: paramaters
###################################################

## samples x genes
exprSix <- t(expr[selSixFinal,])
## paketSize controls memory usage
## (in this example about 1 GB)
## (trade off between time and memory)
paketSize <- 10000
## number of draws with replacement
## should translate to multiples of paketSize-1
nullLength <- 5*paketSize-1
## size of the sliding window (#genes)
windowSize <- 20
## restrict focus on sample subsets
## of size k=1,...,43 to save computing time
maxK <- 43



###################################################
### chunk number 7: calculation eval=FALSE
###################################################
## ## computational expense: about 1d
## ## s _ { k x G' x regulation }
## p <- calculatePValues(expr=exprSix,
##                       windowSize=windowSize,
##                       nullLength=nullLength,
##                       paketSize=paketSize,
##                       maxCluster=maxK,
##                       verbose=TRUE
##                       )
## save( p, file="p_cached.RData", compress=TRUE)


###################################################
### chunk number 8: resultLoading
###################################################
load("p_cached.RData")


###################################################
### chunk number 9: analysis
###################################################

p43up <- p[43,,1]
topSel <- order(p43up)[1:10]
## complete windows
tx <- sort(unique(as.vector(sapply(topSel,function(x) seq(x,x+windowSize-1)))))
numberOfRegions <- sum((tx[-1]-tx[-length(tx)]) > 1) + 1
print(numberOfRegions)



###################################################
### chunk number 10: comparison
###################################################
genes <- names(p43up[tx])
## obtain gene annotation
thisEnv <- paste(chip, "MAP", sep = "")
cytoband <- mget(genes, env = eval(as.symbol(thisEnv)))
thisEnv <- paste(chip, "SYMBOL", sep = "") 
symbol <- mget(genes, env = eval(as.symbol(thisEnv)))

print(cbind(GDash=tx,symbol=unlist(symbol),cytoband=unlist(cytoband)))