%\VignetteIndexEntry{Main vignette (complete version): End-to-end analysis of cell-based screens}
%\VignetteKeywords{Cell based assays}
%\VignettePackage{cellHTS2}


% ---- This is just a stub. -----
% The real .Rnw source for this vignette is in inst/scripts/cellhts2Complete.Rnw% It cannot be put here (inst/doc/) since it takes several minutes to run.
% So it needs to be run manually, and the resulting PDF is zipped and checked into the 
% Bioconductor SVN repository. See the 'Makefile' in this directory, which will 
% copy the zipped PDF from ../scripts/ to here and then unzip it.


\documentclass[11pt]{article}
\usepackage{Sweave}
\begin{document}


\begin{Schunk}
\begin{Sinput}
> library("cellHTS2")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> workPath <- getwd()
> if (!("cellhts.bib" %in% dir())) system(sprintf("cp %s/cellhts.bib .", 
+     system.file("doc", package = "cellHTS2")))
> options(error = recover)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> experimentName <- "KcViab"
> dataPath <- system.file(experimentName, package = "cellHTS2")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> dataPath
> rev(dir(dataPath))[1:12]
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> x <- readPlateList("Platelist.txt", name = experimentName, 
+     path = dataPath)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> x
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> cellHTS2:::tableOutput(file.path(dataPath, "Platelist.txt"), 
+     "plate list")
> cellHTS2:::tableOutput(file.path(dataPath, names(intensityFiles(x))[1]), 
+     "signal intensity", header = FALSE)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> state(x)
\end{Sinput}
\end{Schunk}




\begin{Schunk}
\begin{Sinput}
> out <- writeReport(raw = x, force = TRUE, outdir = tempdir())
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> out
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> browseURL(out)
\end{Sinput}
\end{Schunk}



\begin{Schunk}
\begin{Sinput}
> x <- configure(x, descripFile = "Description.txt", confFile = "Plateconf.txt", 
+     logFile = "Screenlog.txt", path = dataPath)
\end{Sinput}
\end{Schunk}



\begin{Schunk}
\begin{Sinput}
> cellHTS2:::tableOutputWithHeaderRows(file.path(dataPath, 
+     "Plateconf.txt"), "plate configuration", selRows = NULL)
> cellHTS2:::tableOutput(file.path(dataPath, "Screenlog.txt"), 
+     "screen log", selRows = 1:3)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> table(wellAnno(x))
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> png("cellhts2Complete-configurationplot.png", width = 324, 
+     height = 324)
> configurationAsScreenPlot(x)
> dev.off()
\end{Sinput}
\end{Schunk}




\begin{Schunk}
\begin{Sinput}
> xn <- normalizePlates(x, scale = "multiplicative", log = FALSE, 
+     method = "median", varianceAdjust = "none")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> compare2cellHTS(x, xn)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> xsc <- scoreReplicates(xn, sign = "-", method = "zscore")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> xsc <- summarizeReplicates(xsc, summary = "mean")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> scores <- Data(xsc)
> ylim <- quantile(scores, c(0.001, 0.999), na.rm = TRUE)
> boxplot(scores ~ wellAnno(x), col = "lightblue", outline = FALSE, 
+     ylim = ylim)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> y <- scores2calls(xsc, z0 = 1.5, lambda = 2)
> png("cellhts2Complete-callvalues.png")
> plot(Data(xsc), Data(y), col = "blue", pch = ".", xlab = "z-scores", 
+     ylab = "calls", main = expression(1/(1 + e^{
+         -lambda * (z - z[0])
+     })))
> dev.off()
\end{Sinput}
\end{Schunk}




\begin{Schunk}
\begin{Sinput}
> xsc <- annotate(xsc, geneIDFile = "GeneIDs_Dm_HFA_1.1.txt", 
+     path = dataPath)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> cellHTS2:::tableOutput(file.path(dataPath, "GeneIDs_Dm_HFA_1.1.txt"), 
+     "gene ID", selRows = 3:6)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> data("bdgpbiomart")
> fData(xsc) <- bdgpbiomart
> fvarMetadata(xsc)[names(bdgpbiomart), "labelDescription"] <- sapply(names(bdgpbiomart), 
+     function(i) sub("_", " ", i))
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> rnwPath <- system.file("doc/Rnw", package = "cellHTS2")
> setwd(rnwPath)
> system(sprintf("cp biomart.tex %s", workPath))
> setwd(workPath)
\end{Sinput}
\end{Schunk}




\begin{Schunk}
\begin{Sinput}
> xsc
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> save(xsc, file = paste(experimentName, ".rda", sep = ""))
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> setSettings(list(plateList = list(reproducibility = list(include = TRUE, 
+     map = TRUE), intensities = list(include = TRUE, map = TRUE)), 
+     screenSummary = list(scores = list(range = c(-4, 
+         8), map = TRUE))))
> out <- writeReport(raw = x, normalized = xn, scored = xsc, 
+     force = TRUE)
\end{Sinput}
\end{Schunk}










\begin{Schunk}
\begin{Sinput}
> library("Category")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> obsolete <- c("GO:0005489", "GO:0015997", "GO:0045034", 
+     "GO:0005660")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> scores <- as.vector(Data(xsc))
> names(scores) <- geneAnno(xsc)
> sel <- !is.na(scores) & (!is.na(fData(xsc)$go))
> goids <- strsplit(fData(xsc)$go[sel], ", ")
> goids <- lapply(goids, function(x) x[!(x %in% obsolete)])
> genes <- rep(geneAnno(xsc)[sel], listLen(goids))
> cache(categs <- cateGOry(genes, unlist(goids, use.names = FALSE)))
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> nrMem <- rowSums(categs)
> remGO <- which(nrMem < 3 | nrMem > 1000)
> categs <- categs[-remGO, , drop = FALSE]
> nrMem <- rowSums(t(categs))
> rem <- which(nrMem == 0)
> if (length(rem) != 0) categs <- categs[, -rem, drop = FALSE]
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> stats <- scores[sel & (names(scores) %in% colnames(categs))]
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> isDup <- duplicated(names(stats))
> table(isDup)
> dupNames <- names(stats)[isDup]
> sp <- stats[names(stats) %in% dupNames]
> sp <- split(sp, names(sp))
> table(sapply(sp, length))
> aux <- stats[!isDup]
> aux[names(sp)] <- sapply(sp, max)
> stats <- aux
> rm(aux)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> m <- match(colnames(categs), names(stats))
> stats <- stats[m]
> stopifnot(colnames(categs) == names(stats))
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> acMean <- applyByCategory(stats, categs)
> acTtest <- applyByCategory(stats, categs, FUN = function(v) t.test(v, 
+     stats)$p.value)
> acNum <- applyByCategory(stats, categs, FUN = length)
> isEnriched <- (acTtest <= 0.001) & (acMean > 0.5)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> png("cellhts2Complete-volcano.png", width = 180, height = 180)
> par(mai = c(0.9, 0.9, 0.1, 0.1))
> px <- cbind(acMean, -log10(acTtest))
> plot(px, main = "", xlab = expression(z[mean]), ylab = expression(-log[10] ~ 
+     p), pch = ".", col = "black")
> points(px[isEnriched, ], pch = 16, col = "red", cex = 0.7)
> dev.off()
> stopifnot(identical(names(acMean), names(acTtest)), identical(names(acMean), 
+     names(acNum)))
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> enrichedGOCateg <- names(which(isEnriched))
> require("GO.db")
> res <- data.frame(`$n$` = acNum[isEnriched], `$z_{\mbox{\scriptsize mean}}$` = signif(acMean[isEnriched], 
+     2), `$p$` = signif(acTtest[isEnriched], 2), GOID = I(enrichedGOCateg), 
+     Ontology = I(sapply(enrichedGOCateg, function(x) Ontology(get(x, 
+         GOTERM)))), description = I(sapply(enrichedGOCateg, 
+         function(x) Term(get(x, GOTERM)))), check.names = FALSE)
> mt <- match(res$Ontology, c("CC", "BP", "MF"))
> stopifnot(!any(is.na(mt)))
> res <- res[order(mt, res$"$p$"), ]
> cellHTS2:::dataframeOutput(res, header = TRUE, caption = sprintf("Top %d Gene Ontology categories with respect to $z$-score.", 
+     nrow(res)), label = "enrichedGoCateg", gotable = TRUE)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> data2003 <- read.table(file.path(dataPath, "Analysis2003.txt"), 
+     header = TRUE, as.is = TRUE, sep = "\t")
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> i <- data2003$Position + 384 * (data2003$Plate - 1)
> data2003$ourScore <- as.vector(Data(xsc))[i]
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> png("cellhts2Complete-scoresComparison.png", width = 324, 
+     height = 324)
> par(mfrow = c(7, 9), mai = c(0, 0, 0, 0))
> for (i in 1:max(data2003$Plate)) {
+     sel <- (data2003$Plate == i)
+     plot(data2003$ourScore[sel], data2003$Score[sel], 
+         pch = 19, cex = 0.6)
+ }
> dev.off()
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> out <- templateDescriptionFile("template-Description.txt", 
+     force = TRUE)
> out
> readLines(out)
\end{Sinput}
\end{Schunk}


\begin{Schunk}
\begin{Sinput}
> cellHTS2:::tableOutput(file.path(dataPath, "old-Plateconf.txt"), 
+     "cellHTS package-specific plate configuration", selRows = 1:28)
> cellHTS2:::tableOutput(file.path(dataPath, "old-Screenlog.txt"), 
+     "cellHTS package-specific screen log", selRows = 1:3)
\end{Sinput}
\end{Schunk}




\begin{Schunk}
\begin{Sinput}
> library("vsn")
> png("cellhts2Complete-transfplots.png", width = 324, 
+     height = 474)
> 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 <- Data(xn)[, , 1]
> myPlots(dv, main = "untransformed")
> xlog <- normalizePlates(x, scale = "multiplicative", 
+     log = TRUE, method = "median", varianceAdjust = "byExperiment")
> dvlog <- Data(xlog)[, , 1]
> myPlots(dvlog, main = "log2")
> dev.off()
\end{Sinput}
\end{Schunk}




\begin{Schunk}
\begin{Sinput}
> toLatex(sessionInfo())
\end{Sinput}
\end{Schunk}





\end{document}

