twoChannels.R
In cellHTS2: Analysis of cell-based screens - revised version of cellHTS

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

###################################################
### code chunk number 1: setup1
###################################################
library("cellHTS2")


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


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


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


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


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


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


###################################################
### code chunk number 8: plateConfscreenLogTable
###################################################
cellHTS2:::tableOutputWithHeaderRows(file.path(dataPath, "Plateconf.txt"), 
                                     "plate configuration", selRows=NULL, 
                                     preName="twoChannels")

cellHTS2:::tableOutput(file.path(dataPath, "Screenlog.txt"), 
                       "screen log", selRows=1:2, preName="twoChannels")


###################################################
### code chunk number 9: twoChannels.Rnw:198-199
###################################################
table(wellAnno(x))


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

# 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(Data(x))[3])
# channel 1 - no controls
# channel 2 - geneC and geneD
posControls[2] <- "(?i)^geneC$|^geneD$"


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


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


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


###################################################
### code chunk number 14: plateMedianChannels
###################################################
xn <- normalizePlates(x, scale="multiplicative", method="median", 
                      varianceAdjust="none")


###################################################
### code chunk number 15: set cut-off for R1
###################################################
ctoff <- quantile(Data(xn)[,,1], probs=0.05, na.rm=TRUE)


###################################################
### code chunk number 16: FvsRcorrected
###################################################
R <- Data(xn)[,,1]
F <- Data(xn)[,,2]

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

plot(R[,r],F[,r], col= densCols(cbind(R[,r], F[,r])), pch=16,cex=0.7,
     xlab="R1 (log scale)", ylab="R2 (log scale)", log="xy", 
     ylim=ylim, xlim=xlim)

abline(v=ctoff, 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(Data(xn)[,r,1] <= ctoff)
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)
 }


###################################################
### code chunk number 17: summarizeChannels
###################################################
xn1 <- summarizeChannels(xn, fun = function(r1, r2, 
                             thresh=quantile(r1, probs=0.05, na.rm=TRUE)) 
                         ifelse(r1>thresh, r2/r1, as.numeric(NA)))


###################################################
### code chunk number 18: show summarized object
###################################################
dim(Data(xn1))


###################################################
### code chunk number 19: redo plate correction
###################################################
xn1 <- normalizePlates(xn1, scale="multiplicative", log=TRUE, method="median", 
                       varianceAdjust="none") 


###################################################
### code chunk number 20: score and summarize replicates
###################################################
xsc <- scoreReplicates(xn1, sign="-", method="zscore") 
xsc <- summarizeReplicates(xsc, summary="mean") 


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


###################################################
### code 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$"


###################################################
### code chunk number 23: report2Show (eval = FALSE)
###################################################
## setSettings(list(platelist=list(intensities=list(include=TRUE)),
##                  screenSummary=list(scores=list(range=c(-4,4)))))
## out <- writeReport(raw=x, normalized=xn1, scored=xsc, 
##                    outdir="logRatio", 
##                    map=TRUE, 
##                    posControls=posControls, negControls=negControls)


###################################################
### code chunk number 24: report2Do
###################################################
setSettings(list(platelist=list(intensities=list(include=TRUE)),
                 screenSummary=list(scores=list(range=c(-4,4)))))
out <- writeReport(raw=x, normalized=xn1, scored=xsc, 
                   force=TRUE, outdir="logRatio", 
                   map=TRUE,
                   posControls=posControls, negControls=negControls)


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


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


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