###################################################
### chunk number 1: setup
###################################################
library("RbcBook1")


###################################################
### chunk number 2: startup
###################################################
library("prada")
library("facsDorit")
library("geneplotter")


###################################################
### chunk number 3: checkVersions
###################################################
stopifnot(package.version("prada") >= package_version("1.1.13"))
stopifnot(package.version("facsDorit") >= package_version("1.1.0"))


###################################################
### chunk number 4: readViability
###################################################
viab <- read.table(system.file("extdata", "BoutrosKiger.tab",
           package="facsDorit"), header=TRUE, as.is=TRUE, sep="\t")
viab[1:2, ]


###################################################
### chunk number 5: checkViab
###################################################
fourcolumns <- c("Kc1", "Kc2", "S2R1", "S2R2")
stopifnot(all(fourcolumns %in% colnames(viab)))
stopifnot("Plate" %in% colnames(viab))


###################################################
### chunk number 6: 
###################################################
viab$Plate384 <- with(viab, (Plate-min(Plate)) %/% 4)


###################################################
### chunk number 7: bpprun
###################################################
rg  <- 1000*c(250, 2000) 
boxplot(log(Kc1) ~ Plate384, data=viab, outline=FALSE, col="#A6CEE3", 
        ylim=log(rg))


###################################################
### chunk number 8: scprun
###################################################
fac    <- factor(viab$Plate384)
colors <- rainbow(nlevels(fac))[as.integer(fac)]
perm   <- sample(nrow(viab))
plot(viab[perm, c("Kc1", "Kc2")], col=colors[perm], pch=".", 
     log="xy", xlim=rg, ylim=rg)


###################################################
### chunk number 9: normViability
###################################################
library("MASS")
expts  <- c("Kc1", "Kc2", "S2R1", "S2R2")
allMedian <- log(median(as.matrix(viab[, expts])))
for(ex in expts) {
  lmRes <- rlm(log(viab[[ex]]) ~ fac)
  viab[[paste("norm", ex, sep="")]] <- residuals(lmRes) + allMedian
}


###################################################
### chunk number 10: select
###################################################
score <- rowMeans(viab[, paste("normKc", 1:2, sep="")])
sel <- order(score)[1:3] 
viab[sel,1:7]


###################################################
### chunk number 11: bpprnr
###################################################
boxplot(normKc1 ~ Plate384, data=viab, outline=FALSE, col="#1F78B4", 
        ylim=log(rg))


###################################################
### chunk number 12: scprnr
###################################################
plot(viab[perm, c("normKc1", "normKc2")], pch=".", col=colors[perm],
     xlim=log(rg), ylim=log(rg))


###################################################
### chunk number 13: altModel eval=FALSE
###################################################
## lmRes <- rlm(log(luc) ~ expt * fac, data=rViab)


###################################################
### chunk number 14: readFCS1
###################################################
sampleDir <- system.file("extdata", "map", package="facsDorit")
B05 <- readFCS(file.path(sampleDir,  "060304MAPK controls.B05"))
B05


###################################################
### chunk number 15: readFCS2
###################################################
description(B05)[c(130,137,139)]


###################################################
### chunk number 16: readFCS3
###################################################
exprs(B05)[1:2,]


###################################################
### chunk number 17: readCytoSet
###################################################
mapk <- readCytoSet(path=sampleDir, phenoData="plateIndex.txt")
pData(mapk)[1:2,]


###################################################
### chunk number 18: readCytoSet3
###################################################
mapk[[1]]
exprs(mapk[[1]])[1:2,]


###################################################
### chunk number 19: FSCvsSSC
###################################################
x <- exprs(B05)[, c("FSC-H","SSC-H")]
plot(x, pch=20, col=densCols(x))


###################################################
### chunk number 20: fitNorm2
###################################################
nfit <- fitNorm2(x, scalefac=2)


###################################################
### chunk number 21: mahalanobis eval=FALSE
###################################################
## mh <- mahalanobis(x, center=nfit$mu, cov=nfit$S)
## plot(log(nfit$p), mh)


###################################################
### chunk number 22: plotNorm2-2
###################################################
plotNorm2(nfit, ellipse=TRUE)


###################################################
### chunk number 23: plotNorm2-3
###################################################
nfit3 <- fitNorm2(x, scalefac=3)
plotNorm2(nfit3, ellipse=TRUE)


###################################################
### chunk number 24: B05
###################################################
B05.sel <- B05[nfit$sel,]


###################################################
### chunk number 25: FL1vsFL4do1
###################################################
myPlot <- function(x) {
  ex <- exprs(x)[,c("FL1-H","FL4-H")]
  plot(ex, pch=20, col=densCols(ex),
       xlim=range(exprs(B05)[,"FL1-H"]),
       ylim=range(exprs(B05)[,"FL4-H"]))
}


###################################################
### chunk number 26: FL1vsFL4do2
###################################################
myPlot(B05)


###################################################
### chunk number 27: FL1vsFL4do3
###################################################
myPlot(B05.sel)


###################################################
### chunk number 28: FL1vsFL4show eval=FALSE
###################################################
## myPlot <- function(x) {
##   ex <- exprs(x)[,c("FL1-H","FL4-H")]
##   plot(ex, pch=20, col=densCols(ex))
## }
## myPlot(B05)
## myPlot(B05.sel)


###################################################
### chunk number 29: platePlot1
###################################################
nrCells <- csApply(mapk, nrow)
plotPlate(nrCells, nrow=8, ncol=12, main="Cell number",
          col=brewer.pal(9, "YlOrBr"), width=6.3)


###################################################
### chunk number 30: Rggobi eval=FALSE
###################################################
## library("Rggobi")
## x <- exprs(B05)
## gg <- ggobi(x)
## gg$setGlyphs(5,1, 1:nrow(x))
## gg$setColors(rep(9, nrow(x)))
## gg$scatterplot("FL1-H", "FL4-H")
## gg$setBrushColor(5)
## gg$setMode("Brush")


###################################################
### chunk number 31: preprocApo
###################################################
preprocess <- function(x) {
  for(i in 1:length(x)) {
    dat <- exprs(x[[i]])
    fn  <- fitNorm2(dat[, c("FSC-H", "SSC-H")], scalefac=1.5)
    x[[i]] <- dat[fn$sel,]
  }
  return(x)
}
apo <- readCytoSet(path=system.file("extdata", "apoptosis",
         package="facsDorit"),  phenoData="plateIndex.txt")
apoP <- preprocess(apo)


###################################################
### chunk number 32: apo-schematic
###################################################
plot(0,0, type="n", xlab="protein expression", ylab="caspase-3 activation",
xaxt="n", yaxt="n")
abline(v=0, h=0)
text(x=c(-0.5, 0.5, -0.5, 0.5), y=c(0.5, 0.5, -0.5, -0.5),
labels=c("untransfected\napoptotic cells\n(UA)", 
         "transfected\napoptotic cells\n(TA)",
         "untransfected\nnon-apoptotic cells\n(UN)",
         "transfected\nnon-apoptotic cells\n(TN)"), adj=0.5, cex=1.2)


###################################################
### chunk number 33: apo-mock
###################################################
calcthr <- function(x) {h <- hubers(x) ; h$mu+2.5*h$s}
mock <- exprs(apoP[[1]])[,c("FL1-H", "FL4-H")]
plot(mock, pch=".",
     xlab="protein expression", xlim=c(0, 1023), 
     ylab="caspase-3 activation", ylim=c(0, 1023))
thrYFP   <- calcthr(mock[,1])
thrCASP3 <- calcthr(mock[,2])
abline(v=thrYFP, h=thrCASP3, col="red", lty=2)


###################################################
### chunk number 34: cide
###################################################
cide <- exprs(apoP[[6]])[,c("FL1-H", "FL4-H")]


###################################################
### chunk number 35: FisherTest
###################################################
ct   <- thresholds(cide, xthr=thrYFP, ythr=thrCASP3)
fisher.test(ct)


###################################################
### chunk number 36: apo-bcl2
###################################################
plot(cide, pch=".",
     xlab="protein expression", ylab="caspase-3 activation",
     xlim=c(0, 1023), ylim=c(0, 1023))
abline(v=thrYFP, h=thrCASP3, col="red", lty=2)


###################################################
### chunk number 37: calcOdds
###################################################
calcOdds <- function(x){
   ct  <- thresholds(x[,c("FL1-H", "FL4-H")], xthr=thrYFP, ythr=thrCASP3)
   f <- fisher.test(ct)
   res <- -log10(f$estimate)
   return(ifelse((f$p.value > 0.01 | is.infinite(res)), 0, res))
}
odds <- csApply(apoP, calcOdds)


###################################################
### chunk number 38: platePlotApoOdds
###################################################
cols <-  brewer.pal(9, "Reds")[c(rep(1,4), 2:9)]
plotPlate(odds, nrow=8, ncol=12, main="log odds ratios", desc=c("act", ""),
          col=cols, width=6.3, na.action="omit", ind=pData(apo)$wellnr)


###################################################
### chunk number 39: dopreprocessMap
###################################################
mapkP <- preprocess(combineFrames(mapk, factor(pData(mapk)$clone)))


###################################################
### chunk number 40: plotMap
###################################################
par(mfrow=c(1,3))
groups <- c("MEK", "DSPP", "ERK")
names(groups) <- c("a)", "b)", "c)")
library("locfit")
for(i in seq(along=groups)) {
  dat <- exprs(mapkP[[groups[i]]])[, c("FL1-H", "FL4-H")]
  lcft <- locfit.robust(x=dat[, 1], y=dat[, 2], deg=1, alpha=1, maxk=512)
  sel <- sample(1:nrow(exprs(mapkP[[groups[i]]])), 2000)
  plot(dat[sel,], pch=20, col=densCols(dat[sel,]), main=paste(names(groups)[i], groups[i]))
  lines(lcft, col="red", lwd=2)
}


###################################################
### chunk number 41: plotMapSimp eval=FALSE
###################################################
## groups <- c("MEK", "DSPP", "ERK")
## for(i in groups) {
##   dat <- exprs(mapkP[[i]])[, c("FL1-H", "FL4-H")]
##   lcft <- locfit.robust(x=dat[, 1], y=dat[, 2], deg=1, alpha=1, maxk=512)
##   plot(dat, pch=20, col=densCols(dat), main=i)
##   lines(lcft, col="red", lwd=2)
## }


###################################################
### chunk number 42: renameGroups
###################################################
names(groups) <- groups


###################################################
### chunk number 43: calcMapSimp
###################################################
sapply(groups, function (i) {
  dat   <- exprs(mapkP[[i]])[, c("FL1-H", "FL4-H")]
  dlcft <- locfit.robust(x=dat[,1], y=dat[,2], deg=1, alpha=1, deriv=1,
                            maxk=512)
  pp    <- preplot(dlcft, newdata=600, band="local")
  c(delta=pp$fit, zscore=pp$fit/pp$se.fit)
})