%\VignetteIndexEntry{End-to-end analysis of cell-based screens}
%\VignetteKeywords{Cell based assays}
%\VignettePackage{cellHTS}

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{color}
\definecolor{darkblue}{rgb}{0.0,0.0,0.75}
\usepackage[%
baseurl={http://www.bioconductor.org},%
pdftitle={End-to-end analysis of cell-based screens},%
pdfauthor={Wolfgang Huber},%
pdfsubject={cellHTS},%
pdfkeywords={Bioconductor},%
pagebackref,bookmarks,colorlinks,linkcolor=darkblue,citecolor=darkblue,%
pagecolor=darkblue,raiselinks,plainpages,pdftex]{hyperref}

\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\Rpackage}[1]{\textit{#1}}
\newcommand{\Rfunction}[1]{\textit{#1}}
\newcommand{\Rclass}[1]{\textit{#1}}

\newcommand{\myincfig}[3]{%
  \begin{figure}[tp]
    \begin{center}
      \includegraphics[width=#2]{#1}
      \caption{\label{#1}#3}
    \end{center}
  \end{figure}
}

\begin{document}

%------------------------------------------------------------
\title{End-to-end analysis of cell-based screens: from raw intensity readings
  to the annotated hit list}
%------------------------------------------------------------
\author{Michael Boutros, L\'igia Br\'as and Wolfgang Huber}
\maketitle
\tableofcontents

\section{Introduction}
This is a technical report that demonstrates the use of the
\Rpackage{cellHTS} package. It accompanies the paper \textit{Analysis
of cell-based RNAi screens} by Michael Boutros, Ligia Bras and
Wolfgang Huber.  This report explains all the steps necessary to run a
complete analysis of a cell-based high-throughput screen (HTS), from
raw intensity readings to an annotated hit list.

This text has been produced as a reproducible
document~\cite{Gentleman2004RepRes}. It contains the actual computer
instructions for the method it describes, and these in turn produce
all results, including the figures and tables that are shown here. The
computer instructions are given in the language R, thus, in order to
reproduce the computations shown here, you will need an installation
of R (version 2.2 or greater) together with the package
\Rpackage{cellHTS} and some other add-on packages.

To reproduce the computations shown here, you do not need to type them
or copy-paste them from the PDF file; rather, you can take the file
\textit{cellhts.Rnw} in the \textit{doc/Rnw} directory of the package,
open it in a text editor, run it using the R command
\Rfunction{Sweave}, and modify it to your needs.

First, we load the package.
%
<<setup1, results=hide>>=
library("cellHTS")
@ 
%
<<setup2, echo=FALSE, results=hide>>=
## 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)
@ 
%
%------------------------------------------------------------
\section{Reading the intensity data}
\label{sec:read}
%------------------------------------------------------------
We consider a cell-based screen that was conducted in microtiter plate
format, where a library of double-stranded RNAs was used to target the
corresponding genes in cultured \textit{Drosophila} $Kc_{167}$ 
cells~\cite{Boutros2004}. Each of the wells in the plates contains either a
gene-specific probe, a control, or it can be empty. 
The experiments were done in duplicate, and the viability of the cells after
treatment was recorded by a plate reader measuring luciferase activity, which
is indicative of ATP levels. Although this set of example data corresponds to a
single-channel screening assay, the \Rpackage{cellHTS} package can also 
deal with cases where there are readings from more color channels, corresponding 
to different reporters.
Usually, the measurements from each replicate and each color channel 
come in individual result files. The set of available result files and 
the information about them (which plate,
which replicate, which channel) is contained in a spreadsheet, which
we call the \emph{plate list file}. The first few lines of an example
plate list file are shown in Table~\ref{tab:platelist}.
\input{cellhts-platelist}

The first step of the analysis is to read the plate list file, to read
all the intensity files, and to assemble the data into a single R
object that is suitable for subsequent analyses.  The main component
of that object is one big table with the intensity readings of all
plates, channels, and replicates. We demonstrate the R instructions
for this step. First we define the path where the input files can be
found.
%
<<dataPath>>=
experimentName = "KcViab"
dataPath=system.file(experimentName, package="cellHTS") 
@ 
%
In this example, the input files are in the
\Robject{\Sexpr{experimentName}} directory of the \Rpackage{cellHTS}
package. To read your own data, modify \Robject{dataPath} to point to
the directory where they reside. We show the names of 12 files from
our example directory:
%
<<dirDataPath>>=
dataPath
rev(dir(dataPath))[1:12]
@ 
%
and read the data into the object \Robject{x}
%
<<readPlateData, results=hide>>=
x = readPlateData("Platelist.txt", name=experimentName, path=dataPath)
@ 
%
<<showX>>=
x
@ 
%% Create the example table:
%% it would have been nice to use the "xtable" package but it insists on 
%%   adding row numbers, which we don't like.
<<plateFileTable, results=hide, echo=FALSE>>=
cellHTS:::tableOutput(file.path(dataPath, "Platelist.txt"), "plate list")
cellHTS:::tableOutput(file.path(dataPath, names(x$intensityFiles)[1]), "signal intensity",
        header=FALSE, dropColumns=1)
@ 

%------------------------------------------------------------
\section{The \Rclass{cellHTS} class and reports}
%------------------------------------------------------------
The basic data structure of the package is the class
\Rclass{cellHTS}. In the previous section, we have created the
object \Robject{x}, which is an instance of this class. All subsequent
analyses, such as normalization, gene selection and annotation, will
add their results into this object. Thus, the complete analysis
project is contained in this object, and a complete dataset can be
shared with others and stored for subsequent computational analyses in
the form of such an object. In addition, the package offers export
functions for generating human-readable reports, which consist of
linked HTML pages with tables and plots. The final scored hit list
is written as a tab-delimited format suitable for reading
by spreadsheet programs.

To create a report, use the function \Rfunction{writeReport}. It will
create a directory of the name given by \Robject{x\$name} in the
working directory. Alternatively, the argument \Robject{outdir} can be
specified to direct the output to another directory. 
%
<<writeReport1Show, eval=FALSE>>=
writeReport(x)
@ 
<<writeReport1Do, echo=FALSE, results=hide>>=
out = writeReport(x, force=TRUE)
@ 
%
It can take a while to run this function, since it writes a large number
of graphics files. After this function has finished, the index page of
the report will be in the file \Sexpr{file.path(out,"index.html")}, 
and you can view it by directing a web browser to it.
%
<<browseReport1, eval=FALSE>>=
browseURL(file.path(out, "index.html"))
@ 

%------------------------------------------------------------
\section{Annotating the plate results}
%------------------------------------------------------------
\input{cellhts-plateconfiguration} 
\input{cellhts-screenlog} 


The next step of the analysis is to annotate the measured data with
information on controls and to flag invalid measurements. The software
expects the information on the controls in a so-called \emph{plate
configuration file} (see Section~\ref{sec:plateconf}).  This is a
tab-delimited file with one row per well.
%
<<confFile>>=
confFile = file.path(dataPath, "Plateconf.txt")
@ 
%
Selected lines of this file are shown in Table~\ref{tab:plateconfiguration}. 
Individual measurements can be flagged as invalid in the so-called
\emph{screen log file} (see Section~\ref{sec:screenlog}). 
%
<<logFile>>=
logFile  = file.path(dataPath, "Screenlog.txt")
@ 
%
The first 5 lines of this file are shown in Table~\ref{tab:screenlog}. 
The \emph{screen description} file contains a general description of
the screen, its goal, the conditions under which it was performed,
references, and any other information that is pertinent to the
biological interpretation of the experiments. 
<<descripFile>>=
descripFile  = file.path(dataPath, "Description.txt")
@ 
%
We now apply this information to the data object \Robject{x}.
<<annotatePlateRes>>=
x = configure(x, confFile, logFile, descripFile)
@ 
%
Note that the function \Rfunction{configure}\footnote{More precisely,
\Rfunction{configure} is a method for the S3 class \Rclass{cellHTS}.}
takes \Robject{x}, the result from Section~\ref{sec:read}, as an
argument, and we then overwrite \Robject{x} with the result of
this function.
%%
%% Create the example table for plateConf and screenLog
<<plateConfscreenLogTable, results=hide, echo=FALSE>>=
cellHTS:::tableOutput(confFile, "plate configuration", selRows=25:28)
cellHTS:::tableOutput(logFile, "screen log", selRows=1:3)
@ 

%%--------------------------------------------------
\subsection{Format of the plate configuration file}
\label{sec:plateconf}
%%--------------------------------------------------
The software expects this to be a rectangular table in a tabulator
delimited text file, with mandatory columns \emph{Batch}, \emph{Position},
\emph{Well}, \emph{Content}. The \emph{Position} column runs from 1 to the
number of wells in the plate (in the example,
\Sexpr{max(x$plateConf$Position)}), and \emph{Well} is the name of the
corresponding well in letter-number format (in this case, \Robject{A01} to \Robject{P24}). 
The \emph{Content} column
can contain one of the following: \textit{sample} (for wells that
contain genes of interest), \textit{pos} (for positive controls),
\textit{neg} (for negative controls), \textit{empty} (for empty
wells), and \textit{other} (for anything that does not fit into the
four other categories). Note that these annotations are used by the
software in the normalization, quality control, and gene selection
calculations. Data from wells that are annotated as \textit{empty} are
ignored, i.\,e.\ they are set to \Robject{NA}.  Here we look at the frequency of
each well annotation in the example data:
%
<<>>=
table(x$plateConf$Content)
@ 
%
\subsubsection{Multiple plate configurations}
Although it is good practice to use the same plate configuration for
the whole experiment, sometimes this does not work out, and there are
different parts of the experiment with different plate
configurations. It is possible to specify multiple plate
configurations simply by appending them to each other in the plate
configuration file, and marking them with different numbers in the
column \emph{Batch}. 

Note that replicated experiments per plate have to use the same plate
configuration.

%%--------------------------------------------------
\subsection{Format of the screen log file}
\label{sec:screenlog}
%%--------------------------------------------------
The screen log file is a tabulator delimited file with mandatory
columns \emph{Filename}, \emph{Well}, \emph{Flag}. In addition, it can 
contain arbitrary optional columns. Each row corresponds to one flagged
measurement, identified by the filename and the well identifier. The
type of flag is specified in the column \emph{Flag}. Most commonly,
this will have the value ``NA'', indicating that the measurement
should be discarded and regarded as missing.

%------------------------------------------------------------
\section{Normalization and summarization of replicates}
\label{sec:norm}
%------------------------------------------------------------ 
The function \Rfunction{normalizePlateMedian} adjusts for plate
effects by dividing each value in each plate by the median of values
in the plate:
\begin{eqnarray} 
x'_{ki} &=& \frac{x_{ki}}{M_i}\quad\quad\forall k,i \label{eq:normalizePlateMedian}\\
M_{i}&=&\mathop{\operatorname{median}}_{m\in\,\mbox{\scriptsize samples}} x_{mi}
\end{eqnarray}
where $x_{ki}$ is the raw intensity for the $k$-th well in the $i$-th
replicate file, and $x'_{ki}$ is the corresponding normalized intensity. 
The median is calculated across the wells annotated as \textit{sample} in 
the $i$-th result file. This is achieved by calling
%
<<normalizePlateMedian>>=
x = normalizePlateMedian(x) 
@ 
after which the normalized intensities are stored in the slot
\Robject{x\$xnorm}. This is an array of the same size as
\Robject{x\$xraw}. 

We can now summarize the replicates, calculating a single score for
each gene. One option is to take the minimum of the $z$-scores from the two
replicates:
\begin{eqnarray}
z_{ki} &=& \pm \frac{x'_{ki}-\hat{\mu}}{\hat{\sigma}} \label{eq:defz} \\
z_{k}  &=& \min_{i} z_{ki} \label{eq:scoresSummary}.
\end{eqnarray}
Here $\hat{\mu}$ and $\hat{\sigma}$ are estimators of location and
scale of the distribution of $x'_{ki}$.  We use robust estimators,
namely, median and median absolute deviation (mad).  As the values
$x'_{ki}$ were obtained using plate median
normalization~\eqref{eq:normalizePlateMedian}, it holds that
$\hat{\mu}=1$.  The symbol $\pm$ indicates that we allow for either
plus or minus sign in equation~\eqref{eq:defz}; the minus sign can be
useful in the application to an inhibitor assay, where an effect
results in a decrease of the signal and we may want to see this
represented by a large $z$-score.  Then, in
Equation~\eqref{eq:scoresSummary}, the summary is taken over all
replicates for probe $k$.  By using the minimun as the summary
function, the analysis is
particularly conservative: all replicate values have to be high in
order for $z_{k}$ to be high. Depending on the intended stringency of
the analysis, other plausible choices of summary function are the mean
and the maximum.
%
<<summarizeReplicates>>=
x = summarizeReplicates(x, zscore="-", summary="min") 
@ 
%
The resulting single $z$-score value per probe will be stored in the
slot \Robject{x\$score}.
Boxplots of the $z$-scores for the different types of probes are shown
in Figure~\ref{cellhts-boxplotzscore}.
%
<<boxplotzscore, fig=TRUE, include=FALSE, width=4.5, height=5.5>>=
ylim = quantile(x$score, c(0.001, 0.999), na.rm=TRUE)
boxplot(x$score ~ x$wellAnno, col="lightblue", outline=FALSE, ylim=ylim)
@ 
%
\myincfig{cellhts-boxplotzscore}{0.5\textwidth}{Boxplots of $z$-scores 
for the different types of probes.}
%

%------------------------------------------------------------
\subsection{Alternative processing strategies}
%------------------------------------------------------------
The HTML quality report will consider the values in the slot
\Robject{x\$xnorm} for the calculation of its quality metrics.  In the
example above, \Robject{x\$xnorm} contains the data after plate median
normalization, but before calculation of the $z$-scores and the
multiplication by $-1$. The package \Rpackage{cellHTS} allows some
flexibility with respect to these steps. We can already calculate the
$z$-scores and multiply by $-1$ in the function
\Rfunction{normalizePlateMedian}, and then do the summarization
between replicates, by calling the function \Rfunction{summarizeReplicates}
without the argument \Robject{zscore}.
%
<<normalizePlateMedianWithZscore>>=
xalt = normalizePlateMedian(x, zscore="-") 
xalt = summarizeReplicates(xalt, summary="min")
@ 
%
It is easy to define alternative normalization methods, for example,
to adjust for additional experimental biases besides the plate effect.
You might want to start by taking the source code of
\Rfunction{normalizePlateMedian} as a template.

%------------------------------------------------------------
\section{Annotation}
%------------------------------------------------------------
\input{cellhts-geneID} 

Up to now, the assayed genes have been identified solely by the
identifiers of the plate and the well that contains the probe for
them. The \emph{annotation file} contains additional annotation, such
as the probe sequence, references to the probe sequence in public
databases, the gene name, gene ontology annotation, and so forth.
Mandatory columns of the annotation file are \textit{Plate},
\textit{Well}, and \textit{GeneID}, and it has one row for each
well. The content of the \textit{GeneID} column will be species- or
project-specific. The first 5 lines of the example file are shown in
Table~\ref{tab:geneID}, where we have associated each probe with
CG-identifiers for the genes of \textit{Drosophila melanogaster}.
%
<<geneIDs>>=
geneIDFile = file.path(dataPath, "GeneIDs.txt")
x = annotate(x, geneIDFile)
@ 
%% Create the example table:
<<geneIDsTable, results=hide, echo=FALSE>>=
cellHTS:::tableOutput(geneIDFile, "gene ID", selRows = 3:6)
@ 
%
%---------------------------------------------------------------
\subsection{Adding additional annotation from public databases}
%---------------------------------------------------------------
The package \Rpackage{biomaRt} can be used to obtain additional
annotation from public databases~\cite{biomaRt2005}.  You will need
version 1.5.3 or later of the \Rpackage{biomaRt} package. This is
currently (March 2006) in the development branch of Bioconductor and
will be in the release 1.8.

%--------------------------- 
\subsubsection{Installation}
\label{sec:install}
%--------------------------- 
The installation of the \Rpackage{biomaRt} package can be a little bit
tricky, since it relies on the two packages \Rpackage{RCurl} and
\Rpackage{XML}, which in turn rely on the presence of the system
libraries \textit{libcurl} and \textit{libxml2} on your computer. If
you are installing the precompiled R packages (for example, this is
what most people do on Windows), then you need to make sure that the
system libraries on your computer are version-compatible with those on the
computer where the R packages were compiled, and that they are
found. If you are installing the R packages from source, then you need
to make sure that the library header files are available and that the
headers as well as the actual library is found by the compiler and
linker. Please refer to the \textit{Writing R Extensions} manual and
to the FAQ lists on www.r-project.org.

Naturally, before facing an extensive installation struggle, you will
want to explore the \Rpackage{cellHTS} package and see whether it is
worthwhile. We have made available the result of the annotation with
\Rpackage{biomaRt}, described below, as a precomputed R object:
%
<<bdgpbiomart>>=
data("bdgpbiomart")
x$geneAnno = bdgpbiomart
@ 

%------------------------------------------------------------------
\subsubsection{Using biomaRt to annotate the target genes online}
%------------------------------------------------------------------
The commands in this section are optional, you can move on to
Section~\ref{sec:report} if you do not have the \Rpackage{biomaRt}
package or do not want to use it. In the remainder of this section, we
will demonstrate how to obtain the dataframe \Robject{bdgpbiomart} by
querying the online webservice \textit{BioMart} and through it the
Ensembl genome annotation database~\cite{Ensembl2006}.
%
<<biomaRt, eval=FALSE, echo=FALSE, results=hide>>=
Sweave("biomart.Rnw")
@ 
\input{biomart.tex}

% This is how 'bdgpbiomart' it is created from 'x$geneAnno' 
% (before the assignment in the chunk <<bdgpbiomart>> above) from the 
% table 'bmAll' (which is produced by the call to Sweave("biomart.Rnw")
% below):
%
% bdgpbiomart = cbind(x$geneAnno, bmAll)
% save(bdgpbiomart, file="../../../data/bdgpbiomart.rda", compress=TRUE)

%------------------------------------------------------------
\section{Report}
\label{sec:report}
%------------------------------------------------------------
We have now completed the analysis tasks: the dataset has been read, 
configured, normalized, scored, and annotated:
%
<<printxagain>>=
x
@
%
We can now save the data set to a file. 
%
<<savex>>=
save(x, file=paste(experimentName, ".rda", sep=""), compress=TRUE)
@ 
% 
The dataset can be loaded again for subsequent analysis, or passed
on to others. To produce a comprehensive report, we can call the
function \Rfunction{writeReport} again,
%
<<writeReport2, results=hide>>=
writeReport(x, force=TRUE, 
  plotPlateArgs   = list(xrange=c(0.5, 1.5)),
  imageScreenArgs = list(zrange=c( -2, 6.5), ar=1)) 
@ 
%
and use a web browser to view the resulting report
<<browseReport2, eval=FALSE>>=
browseURL(file.path(x$name, "index.html"))
@ 
% 
Now, the report contains a quality report for each plate, and also for
the whole screening assays. The experiment-wide report presents the
$Z'$-factor determined for each experiment (replicate) using the
positive and negative controls~\cite{Oldenburg1999}, the boxplots with
raw and normalized intensities for the different plates, and the
screen-wide plot with the $z$-scores in every well position of each
plate. The latter image plot can also be produced separately using the
function \Rfunction{imageScreen} given in the \Rpackage{cellHTS}
package. This might be useful if we want to select the best display
for our data, namely, the aspect ratio for the plot and/or the range
of $z$-score values to be mapped into the color scale. These can be
passed to the function's arguments \Robject{ar} and \Robject{zrange},
respectively. For example,

<<imageScreen, eval=FALSE, results=hide>>=
imageScreen(x, ar=1, zrange=c(-3,4))
@ 

It should be noted that the per-plate and per-experiment quality
reports are constructed based on the content of \Robject{x\$xnorm}, if
it is present in the \Robject{x} object.  Otherwise, it uses the
content given in the slot \Robject{x\$xraw}.
In the case of dual-channel experiments,  the \Robject{x\$xnorm} slot
could also contain the ratio between the intensities in two different channels, etc. 
The main point that we want to highlight is that \Robject{x\$xnorm} should
contain the data that we want to visualize in the HTML quality reports.
On the other hand, \Robject{x\$score} should always contain the final list of 
scored probes (one value per probe).

At this point we are finished with the basic analysis of the
screen. As one example for how one could continue to further mine the
screen results for biologically relevant patterns, we demonstrate an
application of category analysis.

%------------------------------------------------------------
\section{Category analysis}
%------------------------------------------------------------
We would like to see whether there are Gene Ontology
categories~\cite{GO} overrepresented among the probes with a high
score. For this we use the category analysis from Robert Gentleman's
\Rpackage{Category} package~\cite{GentlemanCategories}. Similar
analyses could be done for other categorizations, for example
chromosome location, pathway membership, or categorical phenotypes
from other studies.

You will need version 1.3.2 or later of the \Rpackage{Category}
package. This is currently (March 2006) in the development branch of
Bioconductor and will be in the release 1.8.
%
<<category, results=hide>>=
library("Category")
stopifnot(package.version("Category") >= package_version("1.3.2"))
@ 
%
Now we can create the category matrix. This a matrix with one column for
each probe and one row for each category. The matrix element
\Robject{[i,j]} is \Robject{1} if probe \Robject{j} belongs to the
\Robject{j}-th category, and \Robject{0} if not.
%
<<cat1>>=
names(x$score) = x$geneAnno$GeneID
sel = !is.na(x$score) & (!is.na(x$geneAnno$go_id))
goids = strsplit(x$geneAnno$go_id[sel], ", ")
genes = rep(x$geneAnno$GeneID[sel], listLen(goids))
categs = cache("categs", cateGOry(genes, unlist(goids, use.names=FALSE)))
@ 
% 
We will selected only those categories that contain at least 3 and no more than 1000 genes.
<<cat2>>=
remGO = which(regexpr("^GO:", nodes(categs)) > 0)
nrMem = listLen(edges(categs)[remGO])
remGO = remGO[nrMem>1000 | nrMem <3]
categs = subGraph(nodes(categs)[-remGO], categs)
@ 
%
As the statistic for the category analysis we use the $z$-score. After
selecting the subset of genes that actually have GO annotation,
%
<<cat3>>=
stats = x$score[ sel & (names(x$score) %in% nodes(categs)) ]
@
%
we are ready to call the category summary functions:
%
<<cat6>>=
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<=1e-3) & (acMean>0.5)
@ 
%
A volcano plot of the $-\log_{10}$ of the $p$-value \Robject{acTtest} versus
the per category mean $z$-score \Robject{acMean} is shown in
Figure~\ref{cellhts-volcano}. For a given category, the $p$-value is calculated from the 
$t$-test against the null hypothesis that there is no difference between the mean
$z$-score of all probes and the mean $z$-score of the probes in that category. 
To select the enriched categories (\Robject{isEnriched}), we considered a significance
level of $0.1\%$ for the $t$-test, and a per category mean $z$-score greater
than $0.5$. This led to the \Sexpr{sum(isEnriched)} categories marked in red in
Figure~\ref{cellhts-volcano} are listed in Table~\ref{tab:enrichedGoCateg}.

 
%
\input{cellhts-enrichedGoCateg} 
%
\myincfig{cellhts-volcano}{0.5\textwidth}{Volcano plot of the $t$-test
  $p$-values and the mean $z$-values of the category analysis for Gene
  Ontology categories. The top categories are shown in red.}
%
<<volcano, fig=TRUE, echo=FALSE, results=hide, include=FALSE, width=3.2, height=3.2>>=
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)
stopifnot(identical(names(acMean), names(acTtest)),  
          identical(names(acMean), names(acNum))) 
@ 
%
<<enrichedGoCateg, echo=FALSE, results=hide>>=
enrichedGOCateg = names(which(isEnriched))

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$"), ]

cellHTS:::dataframeOutput(res, header=TRUE, 
  caption=sprintf("Top %d Gene Ontology categories with respect to $z$-score.", nrow(res)),
  label="enrichedGoCateg", gotable=TRUE)
@ 
%
%------------------------------------------------------------
\section{Appendix: Data transformation}
%------------------------------------------------------------
\myincfig{cellhts-transfplots}{0.95\textwidth}{Comparison between untransformed (left) and 
logarithmically (base 2) transformed (right), normalized data. 
Upper: histogram of intensity values of replicate 1. 
Middle: scatterplots of standard deviation versus mean of the two replicates. 
Bottom: Normal quantile-quantile plots.}

An obvious question is whether to do the statistical analyses on the
original intensity scale or on a transformed scale such as the
logarithmic one.  Many statistical analysis methods, as well as
visualizations work better if (to sufficient approximation)
\begin{itemize}
\item replicate values are normally distributed,
\item the data are evenly distributed along their dynamic range, 
\item the variance is homogeneous along the dynamic range~\cite{Huber2002ismb}.
\end{itemize}

Figure~\ref{cellhts-transfplots} compares these properties for
untransformed and log-transformed normalized data, showing that the difference is small. 
Intuitively, this can be explained by the fact that for
small $x$,
\[
\log(1+x)\approx x
\]
and that indeed the range of the untransformed data is mostly not far
from 1.  Hence, for the data examined here, the choice between
original scale and logarithmic scale is one of taste, rather
than necessity.
%
<<transfplots, fig=TRUE, include=FALSE, width=6.5, height=9>>=
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 = normalizePlateMedian(x, transform=log2)
dvlog = matrix(xlog$xnorm, nrow=prod(dim(xlog$xnorm)[1:2]), ncol=dim(xlog$xnorm)[3])
myPlots(dvlog, main="log2")
@ 


%------------------------------------------------------------
%Bibliography
%------------------------------------------------------------
\bibliography{cellhts}
\bibliographystyle{plain}

\end{document}


