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


###################################################
### chunk number 2: setup2
###################################################
## for debugging:
options(error=recover)


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


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


###################################################
### chunk number 5: showX
###################################################
x


###################################################
### chunk number 6: plateFileTable
###################################################
cellHTS:::tableOutput(file.path(dataPath, "Platelist.txt"), "plate list", preName="twoChannels")


###################################################
### chunk number 7: configure the data
###################################################
x <- configure(x, "Plateconf.txt", "Screenlog.txt",
  "Description.txt", path=dataPath) 


###################################################
### chunk number 8: plateConfscreenLogTable
###################################################
cellHTS:::tableOutput(file.path(dataPath, "Plateconf.txt"), 
    "plate configuration", selRows=1:4, preName="twoChannels")
cellHTS:::tableOutput(file.path(dataPath, "Screenlog.txt"), 
    "screen log", selRows=1:2, preName="twoChannels")


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


###################################################
### chunk number 10: define controls
###################################################
## Define the controls for the different channels:
negControls <- vector("character", length=dim(x$xraw)[4])

# channel 1 - gene A
negControls[1] <- "(?i)^geneA$" # case-insensitive and match the empty string at the beginning and end of a line (to distinguish between "geneA" and "geneAB", for example. Although it is not a problem for the present well annotation)
 
# channel 2 - gene A and geneB
negControls[2] <- "(?i)^geneA$|^geneB$" 

posControls <- vector("character", length=dim(x$xraw)[4])
# channel 1 - no controls
# channel 2 - geneC and geneD
posControls[2] <- "(?i)^geneC$|^geneD$"


###################################################
### chunk number 11: writeReport1Show eval=FALSE
###################################################
## out <- writeReport(x, outdir="raw",posControls=posControls, negControls=negControls)


###################################################
### chunk number 12: writeReport1Do
###################################################
out <- writeReport(x, force=TRUE, outdir="raw", posControls=posControls, negControls=negControls)


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


###################################################
### chunk number 14: plateMedianChannels
###################################################
x <- normalizePlates(x, normalizationMethod="median")


###################################################
### chunk number 15: set cut-off for R1
###################################################
ctoff <- apply(x$xnorm[,,,1], 3, quantile, probs=0.05, na.rm=TRUE)


###################################################
### chunk number 16: FvsRcorrected
###################################################
R <- getMatrix(x$xnorm, channel=1, na.rm=FALSE)
F <- getMatrix(x$xnorm, channel=2, na.rm=FALSE)

# Use the controls of R2 channel:
posC <- which(regexpr(posControls[2], as.character(x$wellAnno), perl=TRUE)>0)
negC <- which(regexpr(negControls[2], as.character(x$wellAnno), perl=TRUE)>0)

ylim <- range(F, na.rm=TRUE)
xlim <- range(R, na.rm=TRUE)

par(mfrow=c(1,ncol(F)), mai=c(1.15,1.15, 0.3,0.3))
for (r in 1:ncol(F)) {
ind <- apply(cbind(R[,r],F[,r]), 1, function(z) any(is.na(z)))
plot(R[!ind,r],F[!ind,r], col= densCols(cbind(R[!ind,r], F[!ind,r])), pch=16,cex=0.7,
     xlab="R1 (log scale)", ylab="R2 (log scale)", log="xy", 
     ylim=ylim, xlim=xlim)
abline(v=ctoff[r], col="grey", lty=2, lwd=2)
points(R[posC,r],F[posC,r], col="red", pch=20, cex=0.9)
points(R[negC,r],F[negC,r], col="green", pch=20, cex=0.9)
ind <- which(x$xnorm[,,r,1] <= ctoff[r])
points(R[ind,r],F[ind,r], col="grey", pch=20, cex=0.9)
#legend("topleft", col=c("grey", "red", "green"), legend=c("masked", "positive controls","negative controls"),  bty="n", pch=20, cex=0.9)
 }


###################################################
### chunk number 17: masking
###################################################
y <- x
for(r in 1:dim(y$xnorm)[3]) {
ind <- y$xnorm[,,r,1]<=ctoff[r]
y$xraw[,,r,][ind] <- NA
}


###################################################
### chunk number 18: Ratio R2/R1
###################################################
y <- normalizeChannels(y, adjustPlates="median")


###################################################
### chunk number 19: alternativeNormalization
###################################################
x = summarizeChannels(x, 
         fun = function(r1, r2, thresh=quantile(r1, probs=0.05, na.rm=TRUE)) 
           ifelse(r1>thresh, log2(r2/r1), as.numeric(NA)),
      adjustPlates="median") 


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


###################################################
### chunk number 21: scores
###################################################
par(mfrow=c(1,2))
ylim <- quantile(x$score, c(0.001, 0.999), na.rm=TRUE)
boxplot(x$score ~ x$wellAnno, col="lightblue", outline=FALSE, ylim=ylim)
imageScreen(x, zrange=c(-2,4))


###################################################
### chunk number 22: RedefineControls
###################################################
## Define the controls for the normalized intensities (only one channel):
# For the single channel, the negative controls are geneA and geneB 
negControls <- "(?i)^geneA$|^geneB$" 
posControls <- "(?i)^geneC$|^geneD$"


###################################################
### chunk number 23: report2Show eval=FALSE
###################################################
## out <- writeReport(x, outdir="logRatio", 
## imageScreenArgs=list(zrange=c(-4,4)),
## plotPlateArgs=list(xrange=c(-3,3)), 
## posControls=posControls, negControls=negControls)


###################################################
### chunk number 24: report2Do
###################################################
out <- writeReport(x, force=TRUE, outdir="logRatio", 
imageScreenArgs=list(zrange=c(-4,4)),
plotPlateArgs=list(xrange=c(-3,3)), 
posControls=posControls, negControls=negControls)


###################################################
### chunk number 25: browse2 eval=FALSE
###################################################
## browseURL(out)


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


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