###################################################
### chunk number 1: setup1
###################################################
library("cellHTS")


###################################################
### chunk number 2: setup2
###################################################
## for debugging:
options(error=recover)
## for software development, when we do not want to install 
## the package after each minor change:
##   for(f in dir("~/huber/projects/Rpacks/cellHTS/R", full.names=TRUE, pattern=".R$"))source(f)


###################################################
### chunk number 3: dataPath
###################################################
experimentName = "KcViab"
dataPath=system.file(experimentName, package="cellHTS") 


###################################################
### chunk number 4: dirDataPath
###################################################
dataPath
rev(dir(dataPath))[1:12]


###################################################
### chunk number 5: readPlateData
###################################################
x = readPlateData("Platelist.txt", name=experimentName, path=dataPath)


###################################################
### chunk number 6: showX
###################################################
x


###################################################
### chunk number 7: plateFileTable
###################################################
cellHTS:::tableOutput(file.path(dataPath, "Platelist.txt"), "plate list")
cellHTS:::tableOutput(file.path(dataPath, names(x$intensityFiles)[1]), "signal intensity",
        header=FALSE, dropColumns=1)


###################################################
### chunk number 8: writeReport1Show eval=FALSE
###################################################
## out = writeReport(x)


###################################################
### chunk number 9: writeReport1Do eval=FALSE
###################################################
## out = writeReport(x, force=TRUE)


###################################################
### chunk number 10: printout eval=FALSE
###################################################
## out


###################################################
### chunk number 11: browseReport1 eval=FALSE
###################################################
## browseURL(out)


###################################################
### chunk number 12: annotatePlateRes
###################################################
x = configure(x, "Plateconf.txt", "Screenlog.txt", 
       "Description.txt", path=dataPath)


###################################################
### chunk number 13: plateConfscreenLogTable
###################################################
cellHTS:::tableOutput(file.path(dataPath, "Plateconf.txt"), 
  "plate configuration", selRows=25:28)
cellHTS:::tableOutput(file.path(dataPath, "Screenlog.txt"), 
  "screen log", selRows=1:3)


###################################################
### chunk number 14: 
###################################################
table(x$plateConf$Content)


###################################################
### chunk number 15: normalizePlateMedian
###################################################
x = normalizePlates(x, normalizationMethod="median") 


###################################################
### chunk number 16: summarizeReplicates
###################################################
x = summarizeReplicates(x, zscore="-", summary="mean") 


###################################################
### chunk number 17: boxplotzscore
###################################################
ylim = quantile(x$score, c(0.001, 0.999), na.rm=TRUE)
boxplot(x$score ~ x$wellAnno, col="lightblue", outline=FALSE, ylim=ylim)


###################################################
### chunk number 18: normalizePlateMedianWithZscore
###################################################
xalt = normalizePlates(x, normalizationMethod="median", zscore="-") 
xalt = summarizeReplicates(xalt, summary="mean")


###################################################
### chunk number 19: geneIDs
###################################################
x = annotate(x, geneIDFile="GeneIDs_Dm_HFA_1.1.txt", path=dataPath)


###################################################
### chunk number 20: geneIDsTable
###################################################
cellHTS:::tableOutput(file.path(dataPath, "GeneIDs_Dm_HFA_1.1.txt"), 
     "gene ID", selRows = 3:6)


###################################################
### chunk number 21: printxagain
###################################################
x


###################################################
### chunk number 22: savex
###################################################
save(x, file=paste(experimentName, ".rda", sep=""), compress=TRUE)


###################################################
### chunk number 23: writeReport2 eval=FALSE
###################################################
## out = writeReport(x, force=TRUE, 
##   plotPlateArgs   = list(xrange=c(0.5, 1.5)),
##   imageScreenArgs = list(zrange=c( -2, 6.5), ar=1)) 


###################################################
### chunk number 24: browseReport2 eval=FALSE
###################################################
## browseURL(out)


###################################################
### chunk number 25: imageScreen eval=FALSE
###################################################
## imageScreen(x, ar=1, zrange=c(-3,4))


###################################################
### chunk number 26: exportData eval=FALSE
###################################################
## writeTab(x, file="Data.txt")


###################################################
### chunk number 27: exportOtherData eval=FALSE
###################################################
## # determine the ratio between each well and the plate median
## y = array(as.numeric(NA), dim=dim(x$xraw))
## nrWell = dim(x$xraw)[1]
## for(p in 1:(dim(x$xraw)[2])) {
##   samples = (x$wellAnno[(1:nrWell)+nrWell*(p-1)]=="sample")
##   y[, p, , ] = apply(x$xraw[, p, , , drop=FALSE], 3:4, 
##      function(w) w/median(w[samples], na.rm=TRUE)) }
## y=signif(y, 4)
## out = matrix(y, nrow=prod(dim(y)[1:2]), ncol=dim(y)[3:4])
## out = cbind(x$geneAnno, x$wellAnno, out)
## colnames(out) = c(names(x$geneAnno), "wellAnno", 
## sprintf("Well/Median_r%d_ch%d", rep(1:dim(y)[3], dim(y)[4]), 
## rep(1:dim(y)[4], each=dim(y)[3])))
## write.tabdel(out, file="WellMedianRatio.txt")


###################################################
### chunk number 28: transfplots
###################################################
library("vsn")
par(mfcol=c(3,2))
myPlots=function(z,...) {
  hist(z[,1], 100, col="lightblue", xlab="",...)
  meanSdPlot(z, ylim=c(0, quantile(abs(z[,2]-z[,1]), 0.95, na.rm=TRUE)), ...)
  qqnorm(z[,1], pch='.', ...)
  qqline(z[,1], col='blue')
}
dv = matrix(x$xnorm, nrow=prod(dim(x$xnorm)[1:2]), ncol=dim(x$xnorm)[3])
myPlots(dv, main="untransformed")
xlog = normalizePlates(x, normalizationMethod="median", transform=log2)
dvlog = matrix(xlog$xnorm, nrow=prod(dim(xlog$xnorm)[1:2]), ncol=dim(xlog$xnorm)[3])
myPlots(dvlog, main="log2")


###################################################
### chunk number 29: sessionInfo
###################################################
toLatex(sessionInfo())