%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
% \VignetteIndexEntry{eSet metadata structures}
%\VignetteDepends{Biobase}
%\VignetteKeywords{Expression Analysis}
%\VignettePackage{Biobase}
\documentclass[12pt]{article}

\usepackage{amsmath,pstricks}
\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}


\textwidth=6.2in
\textheight=8.5in
%\parskip=.3cm
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.3in

\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}

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



\bibliographystyle{plainnat}

\begin{document}
\title{eSet structures, August 2005}
\author{VJ Carey $<$\url{stvjc@channing.harvard.edu}$>$}
\maketitle
\tableofcontents

\section{Preliminaries}

There is a need to expand the capabilities of \Rclass{exprSet}, the primary
multiassay container class in Biobase:
\begin{itemize}
\item we should have the ability to have multiple matrices of assay results,
where rows correspond to reporters, columns correspond to samples
\item we should have the  ability to have fairly rich metadata content
for both samples and reporters
\item we should have the ability to make use of pass-by-reference performance
with large assay data structures
\end{itemize}
We are open to other concerns about container design and would like developers
to consider the following questions:
\begin{itemize}
\item what sorts of reporter data structures and functionalities would
be of interest?  Use cases would be helpful
\item are there compelling concepts of auditability of these structures 
that you would like to address?  I added a character slot `history' in
which one can place a vector of strings, capturing for example
results of match.call()
\item how should we handle geneNames, reporterNames, sampleNames data?  Currently
I have independent slots for reporterNames and sampleNames, but this information
could be stored as rownames/colnames of the reporterInfo and phenoData structures.
\end{itemize}
Please discuss openly on the bioc-devel list.


\section{Introduction}

The high-throughput information that we wish to work with
has the following general structure.
\begin{itemize}
\item There are $N$ samples that have been assayed.  Sample-level
metadata are collected roughly uniformly for all samples, and
take the form of $q$-vectors of mixed continuous and categorical
variables for each sample.
\item There are $G$ reporters on which assay
values are obtained.  $G$ can be as large as $10^7$.
Reporter names come from a vocabulary that can be translated
into information about biological systems.  
\item There are various forms of experiment-level metadata
that are collected.
\end{itemize}

In Bioconductor 1.6,
the \Rclass{exprSet} class managed information of the type
described above, with some restrictions.  In an \Rclass{exprSet}
instance, the \texttt{exprs} slot is intended to hold a $G \times N$
matrix of numerical assay results.  An \texttt{se.exprs} slot
has the same dimensions and structure.  The column names and
row names of the \texttt{exprs} slot provide sample identification
and reporter identification.  Experiment-level metadata can be
contained in slots called \texttt{description}, \texttt{annotation}
and \texttt{notes}.  Sample level data are collected in a structure
of class \texttt{phenoData}.  This has a slot called \texttt{pData}
containing an $N \times q$ data frame, a list of \texttt{varLabels}
that map attribute names to longer descriptive labels, and
a data frame called \texttt{varMetadata} which can capture
variable-level information such as units of measurement.


The \Rclass{eSet} class is a provisional implementation
of an extended container concept.  The current definition
is:

<<gb,results=hide>>=
library(Biobase)
<<echo=FALSE>>=
options(width=60)
<<lkes>>=
getClass("eSet")
package.version("Biobase")
@

The \Rclass{listOrEnv} virtual class confers flexibility
on the assay data container component.
<<lkle>>=
getClass("listOrEnv")
@
The \texttt{assayData} slot can be occupied by either a list or
an environment.  This allows accommodation of multiple $G \times N$
matrices representing assay data.
Environments can be used to obtain pass-by-reference
semantics.  Both list and environment representations
for assay data relax structural
restrictions that exist with an \Rclass{exprSet} instance.

\section{Annotated dataset class}

One superclass from which \Rclass{eSet} inherits is
\Rclass{annotatedDataset}.
<<lkad>>=
getClass("annotatedDataset")
@
A variety of operators such as \verb+[[+ and \verb+$+ are
defined for objects of classes that extend \Rclass{annotatedDataset}.
See the help page on this class.

A key formal container conferred by \Rclass{annotatedDataset} inheritance is
\Rclass{phenoData}.
<<lkpd>>=
getClass("phenoData")
@
The \Rclass{phenoData} slots are \texttt{pData}, the $N \times q$ data frame
holding sample level non-assay information, \texttt{varLabels}, a named
list of $q$ variable labels, and \texttt{varMetadata}, a data frame of
arbitrary dimensions that provides additional information on phenoData variables.
See section \ref{varmet} below for an example.

In addition to the \Rclass{phenoData} container, we now allow a
\texttt{reporterInfo} slot to hold a data.frame of $G$ rows with
metadata about reporters.  The intention is that columns in
\texttt{reporterInfo} will likely hold factors to facilitate
splitting up the reporter set, into control and active reporters,
for example.


\section{Example}

We include a coercion tool, and can thus use the Golub data
to illustrate.  We will add contrived information on
standard errors and presence calls to illustrate the
flexibility.

<<dogol>>=
data(golubMergeSub)
gmes <- as(golubMergeSub, "eSet")
assayData(gmes)[["se.exprs"]] <- runif(length(exprs(gmes)),1.5,2.5)*assayData(gmes)[["exprs"]]
ratpred <- (assayData(gmes)[["exprs"]]/assayData(gmes)[["se.exprs"]] > 2)
assayData(gmes)[["presence"]] <- ratpred
gmes
@
Note that the show method is slightly elaborated to cope with
the added potential content of the assayData component.  The
\Rfunction{exprs} method still works, provided that an element
of the \texttt{assayData} slot has name \texttt{exprs}.
<<lke>>=
dim(exprs(gmes))
@
The use of \texttt{exprs} as an accessor is permitted but does
a bit more than the analogous accessor for \Rclass{exprSet}
objects.  Assay data are, in this example, contained in a
named list of matrices.
<<lkas>>=
dim(assayData(gmes)[["exprs"]])
@
Note the interpretation of \texttt{dim}:
<<lkd>>=
dim(gmes)
@
Subsetting operations perform as expected:
<<lksub>>=
gmes[1:10,1:10]
@
\section{Validity conditions}
<<lkval>>=
getValidity(getClass("eSet"))
@
I would like to have a constraint on reporterNames length
but more discussion is needed.

Let's illustrate the use of the reporterInfo construct.  Suppose
we wish to make it convenient to bind the gene symbols with the
data.
<<dobi>>=
data(bbsym)
gs <- unlist(mget(reporterNames(gmes), bbsym))
gs <- data.frame(name=gs)
rownames(gs) <- reporterNames(gmes)
reporterInfo(gmes) <- gs
@
Now let's restrict our eSet to those reporters for which
symbols contain the substring \texttt{NFK}.
<<donfk>>=
gmes[ grep("NFK", as.character(reporterInfo(gmes)$name)) , ]
@

\section{Environments for assay results}

It is permissible to employ an environment containing matrices
as the assay data container.  Here we construct the \Rclass{eSet}
instance explicitly, notice the direct containment of
sample names and reporter names in slots.
<<doenv>>=
ee <- new.env()
assign("exprs", exprs(golubMergeSub), ee)
gme <- new("eSet", assayData=ee, phenoData=phenoData(golubMergeSub),
 sampleNames=sampleNames(golubMergeSub), reporterNames=geneNames(golubMergeSub))
gme
@
Subsetting works as expected, and the \Rfunction{exprs} or \Rfunction{assayData} extractors
can be used.
<<lki>>=
gme[1:3,1:4]
dim(exprs(gme))
@

\section{phenoData metadata}
\label{varmet}

The \Robject{varMetadata} slot of \Rclass{phenoData} helps to
encode units or other important features of phenoData variables.
<<dovm>>=
pes <- phenoData(gme)
@
It is necessary to initialize the metadata-enabled \Rclass{phenoData}
structure.
<<dopes>>=
pes <- convertVarLabels(pes)
@
Then metadata can be added field by field.
<<addf>>=
pes <- addVarMetadataEntry(pes, "PS", "units", "pct/100")
phenoData(gme) <- pes
getUnits(phenoData(gme), "PS")
@
It would be possible to have a reporter metadata structure as well.
This could be used for preserving information on, e.g., the version
of the mapping between reporter IDs and gene symbols.  Before adding
this feature I want to hear reactions from the community.

%--------------------------------------------------
\section{A suggestion for a basic history mechanism}
%--------------------------------------------------
It is desirable that an \Rclass{eSet} object contains some basic information
about itself. Note that the goal of such a history mechanism is less than
complete reproducibility: for this we have the vignettes (Sweave) technology.

Currently, preprocessing functions in the \Rpackage{affy} and \Rpackage{vsn}
packages already do add some information on what they did to their returned
\Rclass{exprSet} objects. For example, see the last dozen lines in the function
\Rfunction{vsn} in the file vsn/R/vsn.R of the package \Rpackage{vsn}.
Here is a toy example for an \Rclass{eSet} object:
%
<<proposedBasicHist1,echo=TRUE, results=hide>>=
setGeneric("wonderfulPreproc", 
  function(x, lambda) standardGeneric("wonderfulPreproc"))
setMethod("wonderfulPreproc", c("eSet", "numeric"),
  function(x, lambda) {
    y = assayData(x)$exprs
    gamma = mean(log(diag(t(y)%*%y)))
    assayData(x)$exprs = 1/(1+(y^lambda)/gamma)

    cmt = list(
     comment = paste("Processed by wonderful processing method on", date()),
     lambda = lambda,
     gamma = gamma)

    x@description@preprocessing = append(x@description@preprocessing, 
       cmt)
    return(x)
})
Pgme = wonderfulPreproc(gme, pi/2)
@ 
The resulting \Robject{description} slot then looks as follows.
<<proposedBasicHist2>>=
description(Pgme)
description(Pgme)@preprocessing
@

%%%%
\subsection{Open question}
The mechanism proposed above allows to put arbitrary objects into the 
list \Robject{preprocessing}. In particular, these can then be interrogated
and used by subsequent processing functions. For an example, please see the 
first dozen lines in the function
\Rfunction{meanSdPlot} in the file vsn/R/meanSdPlot.R of 
the package \Rpackage{vsn}.

But how does this relate to the \Robject{history} slot in the \Rclass{eSet},
which is simply of class character?

%%%%
\subsubsection{Seth's comments}
I think the \Robject{history} slot should be a list.  I can see two
possible ways to go with this.

Method I: Use automation.  Modify every replacement function for
\Rclass{eSet} to optionally accept a \Robject{history} argument.  The
default is to grab the current or parent call (as a \Rclass{call}
object) and append this to the \Robject{history} list.  The benefit is
we get history tracking without any work.  It remains to be seen
whether this is useful.  

Method II: Require intervention by implementors to record history.
Here we could define an \Rclass{CallLoggingObject} class that has a
\Robject{history} slot (list).  The \Rclass{eSet} class would extend
\Rclass{CallLoggingObject} as would any class wanting to log history.
Usage would look like:

\begin{verbatim}
## record history
logCall(object) <- aCall

## access history
callHistory(object)
\end{verbatim}

The big disadvantage, is that users/developers have to call
\Robject{logCall} or else history is not recorded.

In recent discussions, Robert prefers method I because it doesn't
require changes to preprocessing or other code.  It should be fairly
easy to implement.  We can see if the history log it produces is
useful and iterate.

 
%--------------------------------------------------
\section{Session Information}
%--------------------------------------------------

The version number of R and packages loaded for generating the vignette were:

\begin{verbatim}
<<echo=FALSE,results=tex>>=
toLatex(sessionInfo())
@
\end{verbatim}


\end{document}
