%\VignetteDepends{ALL, Biobase, genefilter, hgu95av2.db}
%\VignetteIndexEntry{Introduction to Bioconductor Lab}
%\VignettePackage{IWB2011}

\documentclass{article}

\usepackage{BioconductorVignette}
\SweaveOpts{keep.source=TRUE,width=9,height=8,eps=FALSE,include=FALSE}

\title{Working with \Bioconductor{} Objects: Microarray Analysis}
\author{Martin Morgan, Chao-Jen Wong \\ Fred Hutchinson Cancer Research Center}
\date{19-21 January, 2010}

\begin{document}
\maketitle

<<setup, echo=FALSE>>=
options(width = 72)
library("ALL")
library("genefilter")
library("hgu95av2.db")
@ 

%% FIXME: citation
\noindent (Adapted from F.\ Hahne and R.\ Gentleman, `The ALL
Dataset', \emph{in} \Bioconductor{} Case
Studies~\cite{biocCaseStudies})

\section{Structures for genomic data: \Rclass{ExpressionSet}}

Genomic data can be very complex, usually consisting of a number of
different bits and pieces. In \Bioconductor{} we have taken the
approach that these pieces should be stored in a single structure to
easily manage the data. The package \Rpackage{Biobase} contains
standardized data structures to represent genomic data. The
\Rclass{ExpressionSet} class is designed to combine several different
sources of information into a single convenient structure. An
\Rclass{ExpressionSet} can be manipulated (e.g., subsetted, copied),
and is the input to or output of many \Bioconductor{} functions.

The data in an \Rclass{ExpressionSet} consist of 
\begin{itemize}
  \item{\code{assayData}:} {Expression data from microarray
    experiments (\code{assayData} is used to hint at the methods used
    to access different data components, as we show below).}
  \item{metadata:} {A description of the samples in the experiment
    (\code{phenoData}), metadata about the features on the chip or
    technology used for the experiment (\code{featureData}), and
    further annotations for the features, for example gene annotations
    from biomedical databases (\code{annotation}).}
  \item{\code{experimentData}:} {A flexible structure to describe the
    experiment.}
\end{itemize}

The \Rclass{ExpressionSet} class coordinates all of these data, so
that you do not usually have to worry about the details. However, an
\Rclass{ExpressionSet} needs to be created in the first place, because
it will be the starting point for many of the analyses using
\Bioconductor{} software.

\Rclass{ExpressionSet} instances are created in one of two
ways. Often, an \Rclass{ExpressionSet} is the output of an \R{}
function. For instance, \Rfunction{justRMA} in the \Biocpkg{affy}
\Bioconductor{} package reads in manufacturer CEL files and outputs an
\Rclass{ExpressionSet}. Alternatively, an \Rclass{ExpressionSet} can
be assembled from its constituent parts; this is illustrated in the
vignette called ``An introduction to \Biocpkg{Biobase} and
\Rclass{ExpressionSets}''. For future reference, this vignette can be
viewed in your web browser with the command
<<biobase-eset-intro,eval=FALSE>>=
browseVignettes("Biobase")
@ 
%% 
or by visiting the \Biocpkg{Biobase} page on the \Bioconductor{} web
site.

Here we use an \Rclass{ExpressionSet} derived from a study of Acute
Lymphoblastic Leukemia (ALL). The data set is available in the data
package \Biocdatapkg{ALL}. The main object in the \Rpackage{ALL}
package is \Robject{ALL}, an instance of \Robject{ExpressionSet}. It
consist of microarray data (chip series HG-U95Av2) from 128
individuals -- 95 sample with B-cell ALL and 33 with T-cell ALL.  The
expression measures have been preprocessed using the robust multichip
average (RMA) method, implemented in the package \Rpackage{affy}.

Let's first load the \Rpackage{ALL} data, and view the \Robject{ALL}
object.
<<load ALL>>=
library(ALL) # attach the ALL package to the search path
data(ALL)    # load the ALL data into the global work space
ALL          # view the ALL instance -- our first ExpressionSet!
@ 
%% 
This prints a short overview of the \Rclass{ExpressionSet}. The
\Rclass{ExpressionSet} is an example of an \Rcode{S4} class. The
authors of the class have provided \emph{methods} for accessing the
data in the class. For instance, 
\begin{description}
\item [\Rcode{exprs(ALL)}] returns the matrix of expression values
  (probe sets as rows, samples as columns). In the ALL data, the
  expression values are pre-processed and log-transformed.
\item [\Rcode{pData(ALL)}] extracts a data frame describing the
  sample phenotype data.
\item [\Rcode{annotation(ALL)}] reports the type of microarray chip
  used in this experiment.
\end{description}
There are additional methods defined to make common operations
easy. For instance, \Rcode{\$} allows one to access columns of
phenotypic data. This
<<BT>>=
ALL$BT
@ 
%% $
retrieves the vector of values in the \Rcode{BT} column of
\Rcode{pData(ALL)}. 

\section{Subsetting}

Another important method available with the \Rclass{ExpressionSet}
class is subsetting. The underlying expression values in an experiment
are represented as a matrix, and the \Rclass{ExpressionSet} reflects
this data -- it can be subset by rows (features) and / or columns
(samples). For instance, to select the samples that have
\Rcode{mol.biol} phenotype \Rcode{NEG}, we might extract the
\Rcode{mol.biol} phenotype data using \Rcode{\$}, and then compare
this to \Rcode{NEG} to produce a vector of logical values that are
\Rcode{TRUE} whenever the value of \Rcode{mol.biol} is \Rcode{NEG}; to
illustrate, we summarize this vector by tabulating the \Rcode{TRUE}
and \Rcode{FALSE} occurrences
<<mol.biol>>=
table(ALL$mol.biol == "NEG")
@ 
%% 
We can use this to create a subset of the original
\Rclass{ExpressionSet} with just these samples
<<eset-subset>>=
idx <- ALL$mol.biol == "NEG"
ALL[,idx]
@ 
%% 
Note that there are only \Sexpr{sum(idx)} samples in this
\Rclass{ExpressionSet}, and that the phenotype and expression data
have been subset in a coordinated fashion.

Additional methods operating on \Rclass{ExpressionSet} can be found on
the help page accessible with
<<eset-help, eval=FALSE>>=
class?ExpressionSet
@ 

\begin{Exercise}
  The objective of this exercise is to create a subset of
  \Robject{ALL} containing only those samples wtih B-cell ALL tumors
  with either BCR/ABL or NEG abnormalities.
  \begin{itemize}
    \item Which samples originate from B-cell tumors?  The
      \Robject{BT} covariate in the phenotypic data encodes the tissue
      type; use \Rfunction{grep} to select any value of \Robject{BT}
      that begins with the letter \Rcode{B}.

    \item Which samples are associated with either the molecular
      subtype BCR/ABL or NEG? The molecular subtypes are labeled in
      the \Robject{mol.biol} phenotypic covariate. To answer this
      question, convert \Robject{ALL\$mol.biol} from a factor to a
      character vector using \Rfunction{as.character}. Then use the
      \Rcode{\%in\%} function (see \Rcode{?"\%in\%"}) to match either
      element of the character vector \Rcode{c("NEG", "BCR/ABL")}.
      
    \item Use the \code{intersect} function to find the intersection
      between the set of samples from B-cell tumors, and the set of
      samples with either NEG or BCR/ABL molecular biology.
      
    \item Finally, create a subset of the \Robject{ALL} object that
      contains B-cell tumors with NEG or BCR/ABL cytotype.
  \end{itemize}
\end{Exercise}
\begin{Solution}
<<get BT>>=
bcell <- grep("^B", as.character(ALL$BT))
types <- c("NEG", "BCR/ABL")
moltyp <- which(as.character(ALL$mol.biol) %in% types)
idx <- intersect(bcell, moltyp)
ALL_bcrneg <- ALL[, idx]
ALL_bcrneg
@ 
\end{Solution}

\Robject{ALL\_bcrne\$mol.biol} is a \Robject{factor}. The `levels' of
a factor are its set of possible values. Since we have reduced the set
of samples, we need to remove the empty levels .
<<factor>>=
ALL_bcrneg$mol.biol <- factor(ALL_bcrneg$mol.biol)
@ 

\section{Non-specific filtering}

A common next step is to use a non-specific filter to remove probe
sets that we know, \emph{a priori} will not be informative in our
analysis. For instance, probe sets that are not annotated to a
particular ENTREZ gene id will not be informative in a study relies on
ENTREZ ids to understand biology. Here we use the \Rfunction{nsFilter}
functon from the \Biocpkg{genefilter} package to filter on a number of
different criteria.
<<filter>>=
library(genefilter)
library(hgu95av2.db) ## annotation package for the ALL dataset
filt_bcrneg <- nsFilter(ALL_bcrneg, require.entrez=TRUE,
                        remove.dupEntrez=TRUE, feature.exclude="^AFFX")
filt_bcrneg$filter.log
@ 
%% 
This removes probes that do not have ENTREZ gene identifiers,
groups of probes that map to the same ENTREZ gene id, and probes that
are annotated as Affymetrix control probes. An entrance into the
literature on non-specific filtering is~\cite{pmid20460310}.

\begin{Exercise}
  Consult the manual for the function \Rfunction{nsFilter} (use
  \Rcode{?nsFilter} after loading the \Biocpkg{genefilter} library)
  for the different options for filtering features from an
  \Robject{ExpressionSet} and the return value of
  \Robject{filt\_bcrneg}.
\end{Exercise}

<<get the eset>>=
ALLfilt_bcrneg <- filt_bcrneg$eset
@ 

\section{Directions}

The text and exercises introduce you to S4 objects, accessors, and
subsetting. You gained practice using the \R{} help system. The
exercises illustrated key data manipulation steps in a microarray work
flows. There are many opportunities for learning more; The
\Bioconductor{} Case Studies book~\cite{biocCaseStudies} and
\Biocpkg{limma} \Bioconductor{} package vignette (for two-color array
pre-processing) are two excellent next steps.

<<sessionInfo,results=tex,echo=FALSE>>=
toLatex(sessionInfo())
@ 

\bibliographystyle{plain}
\bibliography{IWB2011}

\end{document}
