% -*- mode: noweb; noweb-default-code-mode: R-mode; -*-
%\VignetteIndexEntry{lumi}
%\VignetteDepends{lumi, affy, Biobase}
%\VignetteKeywords{Illumina}
%\VignetteKeywords{microarray}
%\VignetteKeywords{microarray normalization}
%\VignetteKeywords{lumi}
%\VignettePackage{lumi}
\documentclass[a4paper]{article}

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

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

\author{Pan Du$^\ddagger$\footnote{dupan@northwestern.edu}, Warren A. Kibbe$^\ddagger$\footnote{wakibbe@northwestern.edu}, Simon Lin$^\ddagger$\footnote{s-lin2@northwestern.edu}}
\begin{document}

\setkeys{Gin}{width=1\textwidth} 

\title{Using lumi, a package processing Illumina Microarray }
\maketitle
\begin{center}$^\ddagger$Robert H. Lurie Comprehensive Cancer Center \\ Northwestern University, Chicago, IL, 60611, USA
\end{center}

\tableofcontents
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Overview of lumi} 
Illumina microarray is becoming a popular microarray platform. The BeadArray technology from Illumina makes its preprocessing and quality control different from other microarray technologies. Unfortunately, until now, most analyses have not taken advantage of the unique properties of the BeadArray system. The \Rpackage{lumi} Bioconductor package especially designed to process the Illumina microarray data. The \Rpackage{lumi} package provides an integrated solution for the bead-level Illumina microarray data analysis. The package covers data input, quality control, variance stabilization, normalization and gene annotation. 

The \Rpackage{lumi} package includes a new variance-stabilizing transformation (VST) algorithm that takes advantage of the technical replicates available on every Illumina microarray. A new robust spline normalization (RSN) algorithm, which combines the features of the quantile and loess normalization, is also implemented in this package. Options available in other popular normalization methods are also provided. Multiple quality control plots are provided in the package. To better annotate the Illumina data, a new, vendor independent nucleotide universal identifier (nuID) was devised to identify the probes of Illumina microarray. The nuID indexed Illumina annotation packages is compatible with other Bioconductor annotation packages. Mappings from Illumina Target Id or Probe Id to nuID are also included in the annotation packages. The output of lumi processed results can be easily integrated with  other microarray data analysis, like differentially expressed gene identification, gene ontology analysis or clustering analysis.


% \section{Citation}


\section{Object models of major classes}
The \Rpackage{lumi} package has one major class: {\bf LumiBatch}. {\bf LumiBatch} is inherited from {\bf ExpressionSet} class in Bioconductor for better compatibility. Their relations are shown in Figure \ref{objects-lumi}. {\bf LumiBatch} class includes {\it se.exprs},  {\it beadNum} and {\it detection} in {\bf assayData} slot for additional informations unique to Illumina microarrays.  A controlData slot is used to keep the control probe information, and a QC slot is added for  keeping the quality control information. The S4 function \Rfunction{plot} supports different kinds of plots by specifying the  specific plot type of {\bf LumiBatch} object. See help of \Rfunction{plot-methods} function for details. The {\it history} slot records all the operations made on the {\bf LumiBatch} object. This provides data provenance. Function \Rfunction{getHistory} is to retrieve the {\it history} slot. Please see the help files of {\bf LumiBatch} class for more details. A series of functions: \Rfunction{lumiR},  \Rfunction{lumiB},  \Rfunction{lumiT},  \Rfunction{lumiN} and \Rfunction{lumiQ} were designed for data input, preprocessing and quality control. Function \Rfunction{lumiExpresso} encapsulates the preprocessing methods for easier usability.


\begin{figure}
\includegraphics{objects-lumi}
\caption{Object models in lumi package}
\label{objects-lumi}
\end{figure}


\section{Data preprocessing}

The first thing is to load the  \Rpackage{lumi} package.
<<eval=T>>=
library(lumi)
@

\subsection{Input data with "intelligent-load" technology}
The \Rfunction{lumiR} function supports directly reading the Illumina raw data output of the Illumina Bead Studio toolkit from version 1 to version 3.  It can automatically detect the BeadStudio output version and format and create a new {\bf LumiBatch} object for it. An example of the input data format is shown in in Figure \ref{fig:dataFileSample}. For simplicity, only part of the data of first sample is shown. The data in the highlighted columns are kept in the corresponding slots of {\bf LumiBatch} object, as shown in Figure \ref{fig:dataFileSample}. The \Rfunction{lumiR} function will automatically determine the starting line of the data. The columns with header including \verb+AVG_Signal+ and \verb+BEAD_STD+ are required for the {\bf LumiBatch} object. The sample IDs and sample labels are extracted from the column names of the data file. For example, based on the column name: \verb+AVG_Signal-1304401001_A+, we will extract \verb+"1304401001"+ as the sample ID and \verb+"A"+ as the sample label (The function assumes the separation of the sample ID and the sample label is \verb+"_"+ if it exists in the column name.). The algorithm will check the uniqueness of sample IDs. If the sample ID is not unique, the entire portion after removing \verb+"AVG_Signal"+ will be used as a sample ID. A file format error will be reported if the uniqueness still cannot be satisfied in this way. The  \Rfunction{lumiR} will automatically initialize the QC slot of the {\bf LumiBatch} object by calling \Rfunction{lumiQ}.

If BeadStudio outputted the control probe data, their information will be kept in the controlData slot of the {\bf LumiBatch} object. If BeadStudio outputted the sample summary information, which is called [Samples Table] in the output text file, the information will be kept in BeadStudioSummay within the QC slot of the {\bf LumiBatch} object.

The BeadStudio can output the gene profile or the probe profile. As the probe profile provides unique mapping from the probe Id to the expression profile, outputting probe profile is preferred. When the probe profile is outputted, as show in Figure \ref{fig:dataFileSample}(B), the ProbeId column will be used as the identifier of {\bf LumiBatch} object.

We strongly suggest outputting the header information when using BeadStudio, as shown in Figure \ref{fig:dataFileSample}. 

If a lumi annotation library is provided, the \Rfunction{lumiR} function will automatically mapping the ProbeId or TargetID as nuID (see annotation section for more details), and keep the mapping information in the featureData of the {\bf LumiBatch} object.
 
\begin{figure}
\includegraphics{dataFileSample}
\caption{An example of the input data format}
\label{fig:dataFileSample}
\end{figure}


\begin{Sinput}
> ## specify the file name
> # fileName <- 'Barnes_gene_profile.txt'  	# Not Run
> ## load the data
> # x.lumi <- lumiR(fileName, lib='lumiHuamnV1')		# Not Run
\end{Sinput}

Here, we just load the pre-saved example data, example.lumi, which is a subset of the experiment data package \Rpackage{lumiBarnes} in the Bioconductor. The example data includes four samples "A01", "A02", "B01" and "B02". "A" and "B" represent different Illumina slides (8 microarrays on each slide), and "01" and "02" represent different samples.  That means "A01" and "B01" are technique replicates at different slides, the same for "A02" and "B02". 

\begin{Sinput}
> ## load example data
> data(example.lumi)	
> ## summary of the example data
> example.lumi
\end{Sinput}


<<echo=FALSE>>=
data(example.lumi)		
example.lumi	
@ 

\subsection{Quality control of the raw data}
The quality control of a {\bf LumiBatch} object includes a data summary (the mean and standard deviation, sample correlation,  detectable probe ratio of each sample (microarray)) and different quality control plots.

In the featureData, it also records the presentCount of each probe, which measures the number of detectable measurements (the detection p-value less than the user provided threshold (0.01 by default)) of each probe.

\Rfunction{LumiQ} function will produce the data summary of a {\bf LumiBatch} object and organize the results in a QC slot of {\bf LumiBatch} object. When creating the {\bf LumiBatch} object, the \Rfunction{LumiQ} function will be called to initialize the QC slot of the {\bf LumiBatch} object.

Summary of the quality control information of example.lumi data. If the QC slot of the {\bf LumiBatch} object is empty, function \Rfunction{lumiQ} will be automatically called to estimate the quality control information.
<<echo=T>>=
summary(example.lumi, 'QC')	## summary of the quality control
@

The S4 method \Rfunction{plot} can produce the quality control plots of {\bf LumiBatch} object. The quality control plots includes: the density plot (Figure \ref{fig:density}), box plot (Figure \ref{fig:box}), pairwise correlation between microarrays (Figure \ref{fig:pairs}), pairwise MAplot between microarrays (Figure \ref{fig:MAplot}), density plot of coefficient of varience, (Figure \ref{fig:cv}), and the sample relations (Figure \ref{fig:sampleRelation1}). More details are in the help of \Rfunction{plot,LumiBatch-method} function. Most of these plots can also be plotted by the extended general functions: \Rfunction{hist} (for density plot), \Rfunction{boxplot}, \Rfunction{MAplot}, \Rfunction{pairs}.


% density plot
Plot the density plot of the {\bf LumiBatch} object. See Figure \ref{fig:density}.
\begin{Sinput}
> ## plot the density
> plot(example.lumi, what='density')			
> ## or
> hist(example.lumi)
\end{Sinput}

\begin{figure}
\centering
<<fig=true, width=10, height=8, quiet=T, echo=F>>=
plot(example.lumi, what='density')	## plot the density
@
\caption{Density plot of Illumina microarrays before normalization}
\label{fig:density}
\end{figure}            

% box plot
Plot the box plot of the {\bf LumiBatch} object. See Figure \ref{fig:box}.
\begin{Sinput}
> ## plot the box plot
> plot(example.lumi, what='boxplot')			
> ## or
> boxplot(example.lumi)
\end{Sinput}

\begin{figure}
\centering
<<fig=true, width=10, height=8, quiet=T, echo=F>>=
plot(example.lumi, what='boxplot') 	## box plot
@
\caption{Density plot of Illumina microarrays before normalization}
\label{fig:box}
\end{figure}            

% pair plot with microarray correlation
Plot the pairwise sample correlation of the {\bf LumiBatch} object. See Figure \ref{fig:pairs}.
\begin{Sinput}
> ## plot the pair plot
> plot(example.lumi, what='pair')			
> ## or
> pairs(example.lumi)
\end{Sinput}

\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T, echo=F>>=
plot(example.lumi, what='pair')		## pairwise plots
@
\caption{Pairwise plot with microarray correlation before normalization}
\label{fig:pairs}
\end{figure}            

% MA plot
The MA plot of the {\bf LumiBatch} object. See Figure \ref{fig:MAplot}.
\begin{Sinput}
> ## plot the MAplot
> plot(example.lumi, what='MAplot')			
> ## or
> MAplot(example.lumi)
\end{Sinput}

\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T, echo=F>>=
plot(example.lumi, what='MAplot')	## pairwise MAplot
@
\caption{Pairwise MAplot before normalization}
\label{fig:MAplot}
\end{figure}            


% plot the density plot of coefficient of variance
The density plot of the coefficient of variance of the {\bf LumiBatch} object. See Figure \ref{fig:cv}.
\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T>>=
plot(example.lumi, what='cv')	## density plot of coefficient of varience
@
\caption{Density Plot of Coefficient of Varience}
\label{fig:cv}
\end{figure}            


% sample relations, hierarchical clustering
Plot sample relations using hierarchical clustering, see Figure \ref{fig:sampleRelation1}
\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T>>=
plot(example.lumi, what='sampleRelation')
@
\caption{Sample relations before normalization}
\label{fig:sampleRelation1}
\end{figure}            

% sample relations, MDS
Plot the sampleRelation using MDS, see Figure \ref{fig:sampleRelation2}. The color of the sample is based on the sample type, which is \verb+"01", "02", "01", "02"+ for the sample data. Please see the help of \Rfunction{getSampleRelation} and \Rfunction{plot-methods} for more details.
\begin{Sinput}
> ## plot the sample relations
> plot(example.lumi, what='sampleRelation', method='mds', color=c("01", "02", "01", "02"))	
> ## or
> plotSampleRelation(example.lumi, method='mds', color=c("01", "02", "01", "02"))
\end{Sinput}

\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T, echo=F>>=
plot(example.lumi, what='sampleRelation', method='mds', color=c("01", "02", "01", "02"))
@
\caption{Sample relations before normalization}
\label{fig:sampleRelation2}
\end{figure} 


\subsection{Background correction}
The \Rpackage{lumi} package provides \Rfunction{lumiB} function for background correction. As both \Rfunction{vst} and \Rfunction{log2} transforms require the expression value to be positive. The default background correction method ('forcePositive') will force all expression values to be positive by adding an offset (minus minimum value plus one). It does nothing if all expression values are positive. Other options of \Rfunction{lumiB} include 'none' and 'bg.adjust'. 'bg.adjust' method is based on the \Rfunction{bg.adjust} function in \Rpackage{affy} package. User can also provide their own background correction function with a LumiBatch Object as the first argument and return a LumiBatch Object with background corrected. See \Rfunction{lumiB} help document for more details.


\subsection{Variance stabilizing transform}
Variance stabilization is critical for subsequent statistical inference to identify differential genes from microarray data. We devised a variance-stabilizing transformation (VST) by taking advantages of larger number of technical replicates available on the Illumina microarray. Please see [1] for details of the algorithm.

Function \Rfunction{lumiT} performs variance stabilizing transform with both input and output being {\bf LumiBatch} object.

Do default VST variance stabilizing transform
<<>>=
## Do default VST variance stabilizing transform
lumi.T <- lumiT(example.lumi)
@ 

The \Rfunction{plotVST} can plot the transformation function of VST, see Figure \ref{fig:VST}, which is close to log2 at high expression values, see Figure \ref{fig:VST_log2}. Function \Rfunction{lumiT} also provides options to do \verb+"log2" or "cubicRoot"+ transform. See help of \Rfunction{lumiT} for details.

\begin{Sinput}
> ## plot VST transformation
> trans <- plotVST(lumi.T)
> ## compare the log2 and VST transform
> matplot(log2(trans$untransformed), trans$transformed, main='compare VST and log2 transform', xlab='log2 transformed', ylab='vST transformed')
\end{Sinput}

\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T, echo=F>>=
trans <- plotVST(lumi.T)
@
\caption{VST transformation}
\label{fig:VST}
\end{figure} 

\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T, echo=F>>=
matplot(log2(trans$untransformed), trans$transformed, type='l', main='compare VST and log2 transform', xlab='log2 transformed', ylab='vST transformed')
abline(a=0, b=1, col=2)
@
\caption{Compare VST and log2 transform}
\label{fig:VST_log2}
\end{figure} 


\subsection{Data normalization}
We proposed a robust spline normalization (RSN) algorithm, which combines the features of quanitle and loess nor-malization, is designed to normalize the variance-stabilized data. Please see [1] for details of the algorithm.

Function \Rfunction{lumiN} performs robust spline normalization (RSN) algorithm with both input and output being {\bf LumiBatch} object. \Rfunction{lumiN} also provides options to do \verb+"loess", "quantile", "vsn"+ normalization. See help of \Rfunction{lumiN} for details.

Do default RSN between microarray normaliazation
<<>>=
## Do RSN between microarray normaliazation
lumi.N <- lumiN(lumi.T)
@ 

User can also easily select other normalization method. For example, the following command will run quantile between microarray normaliazation.
<<eval=FALSE>>=
## Do quantile between microarray normaliazation
lumi.N <- lumiN(lumi.T, method='quantile')	## Not Run
@ 


\subsection{Quality control after normalization}
To make sure the data quality meets our requirement, we do a second round of quality control of normalized data with different QC plots. Compare the plots before and after normalization, we can clearly see the improvements.

<<>>=
## Do quality control estimation after normalization
lumi.N.Q <- lumiQ(lumi.N)

## summary of the quality control
summary(lumi.N.Q, 'QC')		## summary of QC
@

% density plot
\begin{figure}
\centering
<<fig=true, width=10, height=8, quiet=T>>=
plot(lumi.N.Q, what='density')		## plot the density
@
\caption{Density plot of Illumina microarrays after normalization}
\label{fig:density.N}
\end{figure}            

% box plot
\begin{figure}
\centering
<<fig=true, width=10, height=8, quiet=T>>=
plot(lumi.N.Q, what='boxplot')		## box plot
# boxplot(lumi.N.Q)
@
\caption{Density plot of Illumina microarrays after normalization}
\label{fig:box.N}
\end{figure}            

% pairwise plot with microarray correlation
\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T>>=
plot(lumi.N.Q, what='pair')			## pairwise plots
@
\caption{Pairwise plot with microarray correlation after normalization}
\label{fig:pairs.N}
\end{figure}            

% MA plot
\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T>>=
plot(lumi.N.Q, what='MAplot')		## plot the pairwise MAplot
@
\caption{Pairwise MAplot after normalization}
\label{fig:MAplot.N}
\end{figure}            

% sample relations, hierarchical clustering
\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T>>=
## plot the sampleRelation using hierarchical clustering
plot(lumi.N.Q, what='sampleRelation')
@
\caption{Sample relations after normalization}
\label{fig:sampleRelation1.N}
\end{figure}            

% sample relations, MDS
\begin{figure}
\centering
<<fig=true, width=10, height=10, quiet=T>>=
## plot the sampleRelation using MDS
plot(lumi.N.Q, what='sampleRelation', method='mds', color=c("01", "02", "01", "02"))
@
\caption{Sample relations after normalization}
\label{fig:sampleRelation2.N}
\end{figure}            


\subsection{Encapsulate the processing steps}
The \Rfunction{lumiExpresso} function is to encapsulate the major functions of Illumina preprocessing. It is organized in a similar way as the \Rfunction{expresso} function in \Rpackage{affy} package. The following code basically did the same processing as the previous multi-steps and produced the same results lumi.N.Q.

<<echo=TRUE>>=
## Do all the default preprocessing in one step
lumi.N.Q <- lumiExpresso(example.lumi)
@ 

Users can easily customize the processing parameters. For example, if the user want to do "quantile" normalization instead of "rsn" normalization, the user can run the following code. For more details, please read the help document of \Rfunction{lumiExpresso} function. 

<<echo=TRUE>>=
## Do all the default preprocessing in one step
lumi.N.Q <- lumiExpresso(example.lumi, normalize.param=list(method='quantile'))
@ 


\section{Gene annotation}
Illumina microarray provides the TargetID or the ProbeID to identify the measurements. The TargetID is used as a public identifier by Illumina and is supposed to be stable. The problem of the TargetID is that it can correspond to several different probes, which are supposed to match the same gene.  Due to the binding affinity difference or alternative splicing, the probes corresponding the the sample TargetID may have quite different expression levels and patterns. If we use TargetID to identify the measurements, then we cannot differentiate the difference between these probes. Another problem of using the TargetID is that the mapping between the TargetID and probes could be changed with our better understanding of the gene. Moreover, the TargetID used by Illumina microarray is not consistent among different versions of arrays. For instance, the same 50mer sequence has two different TargetIDs used by Illumina: \verb+"GI_21070949-S"+ in the \verb+Mouse_Ref-8_V1+ chip and \verb+"scl022190.1_154-S"+ in the \verb+Mouse-6_V1+ chip. This causes difficulties when combining clinical microarray data collected over time using different versions of the chips. 

In order to get unique mapping between microarray measurements and probes, using ProbeID is preferred. However, the ProbeID of Illumina is not stable. It is changing between different versions, even between different batches of Illumina microarrays. To solve these problems, we designed a nucleotide universal identifier (nuID), which encodes the 50mer oligonucleotide sequence and contains error checking and self-identification code. By using nuID, all the problems mentioned above can be easily solved. For details, please read [2]. 


\subsection{Examples of nuID}
\begin{Sinput}
> ## provide an arbitrary nucleotide sequence as an example 
> seq <- 'ACGTAAATTTCAGTTTAAAACCCCCCG'
> ## create a nuID for it
> id <- seq2id(seq)
> print(id)
\end{Sinput}

<<echo=F>>=
## provide an arbitrary nucleotide sequence as an example 
seq <- 'ACGTAAATTTCAGTTTAAAACCCCCCG'
## create a nuID for it
id <- seq2id(seq)
print(id)
@


The original nucleotide sequence can be easily recovered by \Rfunction{id2seq}
<<>>=
id2seq(id)
@

The nuID is self-identifiable. \Rfunction{is.nuID} can check the sequence is nuID or not. 
A real nuID
<<>>=
is.nuID(id)
@

An random sequence
<<>>=
is.nuID('adfqeqe')
@

\subsection{Illumina microarray annotation package}
Because all the Illumina microarrays use 50-mers, by using the nuID universal identifier, we are able to build one annotation database for different versions of the human (or other species) chips. Moreover, the nuID can be directly converted to the probe sequence, and used to get the most updated refSeq matches and annotations.  Annotation packages indexed by nuID for different Illumina expression chips can be downloaded from Bioconductor. 

The Illumina annotation packages are produced by using \Rpackage{AnnBuilder} with small modification. As a result, the format of the package is the same as Affymetrix annotation package, lots of packages designed for Affymetrix can also be used for Illumina annotation package. The mappings between TargetID to nuID and ProbeID to nuID are also included in the Illumina annotation packages. Thus, we can easily mapping between the nuID and TargetID or ProbeID.

Here is some examples:
\begin{Sinput}
> ## load lumi annotation package
> lib <- 'example.lumi'		# Huamn lumi annotation package
> if(require(GO) & require(annotate) & require(lib, character.only=TRUE)) {
+ 	GOId <- 'GO:0004816' 		 # asparagine-tRNA ligase activity
+ 	probe <- lookUp(GOId, lib, 'GO2ALLPROBES')
+ 	# probes under 'GO:0004816'  category
+ 	print(probe)
+ }
\end{Sinput} 

<<echo=F>>=
## load lumi annotation package
lib <- 'lumiHumanV1'		# Huamn lumi annotation package
if(require(GO) & require(annotate) & require(lib, character.only=TRUE)) {
	GOId <- 'GO:0004816' 		 # asparagine-tRNA ligase activity
	probe <- lookUp(GOId, lib, 'GO2ALLPROBES')
	# probes under 'GO:0004816'  category
	print(probe)
}
@


\begin{Sinput}
> # specify a nuID
> nuId <- 'WVUU7XyNw3ucXzwdEk'	
> if (require(annotate) & require(lib, character.only=TRUE)) {
+	# get the gene symbol of nuId
+	getSYMBOL(nuId, lib)
+ }
\end{Sinput} 


<<echo=F>>=
# specify a nuID
nuId <- 'WVUU7XyNw3ucXzwdEk'	
if (require(annotate) & require(lib, character.only=TRUE)) {
	# get the gene symbol of nuId
	getSYMBOL(nuId, lib)
}
@ 


Mapping from nuID to TargetID
<<echo=T>>=
nuId <- "WVUU7XyNw3ucXzwdEk"
nuID2targetID(nuId, lib='lumiHumanV1')
@ 

Mapping from TargetID to nuID
<<echo=T>>=
targetID <- "GI_7262387-S"	
targetID2nuID(targetID, lib='lumiHumanV1')
@ 

\subsection{Transfer Illumina identifier annotated data into nuID annotated}
As the annotation packages include the mappings between TargetID to nuID and ProbeID to nuID. We can easily map the targetID (or Probe Id) to nuID. The function can automatically check whether targetID or Probe Id  was used in the text data file, and convert them as nuID. Function \Rfunction{addNuId2lumi} can transfer a TargetID or Probe Id indexed  {\bf LumiBatch} object as an nuID indexed {\bf LumiBatch} object. And the mapping between the nuID and TargetID is kept in the featureData of the {\bf LumiBatch} object. If a {\bf LumiBatch} object has already been nuID indexed, the function will do nothing.
<<quiet=T>>=
if (require(lumiHumanV1)) {
	lumi.N <- addNuId2lumi(lumi.N, lib='lumiHumanV1')
}
@

The {\bf LumiBatch} object can also be directly transferred as nuID indexed at the very beginning of inputting data using \Rfunction{lumiR}. For example:
\begin{Sinput}
> ## load the data
>  example.lumi <- lumiR(fileName, lib='lumiHumanV1')	 # Not run
\end{Sinput} 


\section{A use case: from raw data to functional analysis}
Figure \ref{useCase} shows the data processing flow chart of the use case. Since the classes in \Rpackage{lumi} package are inherited from class {\bf ExpressionSet}, packages and functions compatible with class {\bf ExpressionSet} or accepting matrix as input   all can be used for \Rpackage{lumi} results. Here we just give two examples: using \Rpackage{limma} to identify differentiated genes and using \Rpackage{GOstats} to annotate the significant genes.

We use the Barnes data set [3] as an example, which has be created as a Bioconductor experiment data package \Rpackage{lumiBarnes}.  The Barnes data set measured a dilution series of two human tis-sues, blood and placenta.  It includes six samples with the titration ratio of blood and placenta as 100:0, 95:5, 75:25, 50:50, 25:75 and 0:100. The samples were hybridized on HumanRef-8 BeadChip (Illumina, Inc) in duplicate. We select samples with titration ratio, 100:0 and 95:5 (each has two technique replicates) in  this data set to evaluate the detection of differential expressions.


\begin{figure}
\includegraphics{useCaseFlowChart}
\caption{Flow chart of the use case}
\label{useCase}
\end{figure}


\subsection{Preprocess the Illumina data}
\begin{Sinput}
> library(lumi)
> ## specify the file name
> # fileName <- 'Barnes_gene_profile.txt'  #  Not run
> ## load the data
> # example.lumi <- lumiR(fileName, lib='lumiHumanV1')	 # Not run
 
> ## load saved data
> load(example.lumi) 
> ## sumary of the daa
> example.lumi
> ## summary of quality control information
> summary(example.lumi, 'QC')

> ## preprocessing and quality control after normalization
> lumi.N.Q <- lumiExpresso(example.lumi, QC.evaluation=TRUE)
> ## summary of quality control information after preprocessing
> summary(lumi.N.Q, 'QC')

\end{Sinput}

\subsection{Identify differentiate genes}
Identify the differentiated genes based on moderated t-test using \Rpackage{limma}.

Retrieve the normalized data
<<>>=
dataMatrix <- exprs(lumi.N)
@

To speed up the processing and reduce false positives, remove the unexpressed genes
<<>>=
presentCount <- pData(featureData(lumi.N))$presentCount
selDataMatrix <- dataMatrix[presentCount > 0,]
selProbe <- rownames(selDataMatrix)
@


\begin{Sinput}
> ## Specify the sample type
> sampleType <- c('100:0', '95:5', '100:0', '95:5')
> if (require(limma)) {
+ 	## compare '95:5' and '100:0'
+ 	design <- model.matrix(~ factor(sampleType))
+ 	colnames(design) <- c('100:0', '95:5-100:0')
+ 	fit <- lmFit(selDataMatrix, design)
+ 	fit <- eBayes(fit)
+       ## Add gene symbols to gene properties
+	if (require(lumiHumanV1) & require(annotate)) {
+             geneSymbol <- getSYMBOL(fit$genes$ID, 'lumiHumanV1')
+              fit$genes <- data.frame(fit$genes, geneSymbol=geneSymbol)
+        }
+ 	## print the top 10 genes
+ 	topTable(fit, coef='95:5-100:0', adjust='fdr', number=10)
+ 	
+ 	## get significant gene list with FDR adjusted p.values less than 0.01
+ 	p.adj <- p.adjust(fit$p.value[,2])		
+ 	sigGene.adj <- selProbe[ p.adj < 0.01]
+ 	## without FDR adjustment
+ 	sigGene <- selProbe[ fit$p.value[,2] < 0.001]
+
+ } 
\end{Sinput}

<<echo=F>>=
## Specify the sample type
sampleType <- c('100:0', '95:5', '100:0', '95:5')
if (require(limma)) {
	##  compare '95:5' and '100:0'
	design <- model.matrix(~ factor(sampleType))
	colnames(design) <- c('100:0', '95:5-100:0')
	fit <- lmFit(selDataMatrix, design)
	fit <- eBayes(fit)
	## Add gene symbols to gene properties
	if (require(lumiHumanV1) & require(annotate)) {
              geneSymbol <- getSYMBOL(fit$genes$ID, 'lumiHumanV1')
              fit$genes <- data.frame(fit$genes, geneSymbol=geneSymbol)
          }
	## print the top 10 genes
	print(topTable(fit, coef='95:5-100:0', adjust='fdr', number=10))
	
	## get significant gene list with FDR adjusted p.values less than 0.01
	p.adj <- p.adjust(fit$p.value[,2])		
	sigGene.adj <- selProbe[ p.adj < 0.01]
	## without FDR adjustment
	sigGene <- selProbe[ fit$p.value[,2] < 0.01]
} 
@

Based on the significant genes identified using \Rpackage{limma} or t-test, we can do further analysis, like GO analysis (\Rpackage{GOstats} package) and machine learning (\Rpackage{MLInterface} package). Next, we will use GO analysis as an example.


\subsection{Gene Ontology analysis}
Based on the significant genes identified using \Rpackage{limma} or t-test, we can further do Gene Ontology annotation. We can use package \Rpackage{GOstats} to do the analysis.

Do Hypergeometric test of Gene Ontology based on the significant gene list (for e. Table \ref{ta:GOggterms} shows the significant GO terms of Molecular Function with p-value less than 0.01. Here only show the significant GO terms of BP (Biological Process). For other GO categories MF(Molecular Function) and CC (Cellular Component), it just follows the same procedure.

\begin{Sinput}
> if (require(GOstats) & require(lumiHumanV1)) {
+ 
+       ## Get the locuslink Id of the gene
+       sigLL <- unique(unlist(mget(sigGene, env=lumiHumanV1LOCUSID, ifnotfound=NA)))
+       sigLL <- as.character(sigLL[!is.na(sigLL)])
+       params <- new("GOHyperGParams",
+               geneIds= sigLL,
+               annotation="lumiHumanV1",
+               ontology="BP",
+               pvalueCutoff= 0.01,
+               conditional=FALSE,
+               testDirection="over")
+          
+       hgOver <- hyperGTest(params)
+          
+       ## Get the p-values of the test
+       gGhyp.pv <- pvalues(hgOver)
+ 	
+       ## select the Go terms with p-value less than 0.001
+       sigGO.ID <- names(gGhyp.pv[gGhyp.pv < 0.001])
+ 	
+       ## Here only show the significant GO terms of BP (Molecular Function)
+       ## 	For other categories, just follow the same procedure.
+       sigGO.Term <- getGOTerm(sigGO.ID)[["BP"]]
+ }
\end{Sinput}


<<echo=F, eval=F>>=
if (require(GOstats) & require(lumiHumanV1)) {

	## Get the locuslink Id of the gene
	sigLL <- unique(unlist(mget(sigGene, env=lumiHumanV1ENTREZID, ifnotfound=NA)))
	sigLL <- as.character(sigLL[!is.na(sigLL)])
	params <- new("GOHyperGParams",
              geneIds= sigLL,
              annotation="lumiHumanV1",
              ontology="BP",
              pvalueCutoff= 0.01,
              conditional=FALSE,
              testDirection="over")
         
         hgOver <- hyperGTest(params)
         	
	## Get the p-values of the test
	gGhyp.pv <- pvalues(hgOver)
	
	## select the Go terms with p-value less than 0.001
	sigGO.ID <- names(gGhyp.pv[gGhyp.pv < 0.001])
	
	## Here only show the significant GO terms of BP (Molecular Function)
	## 	For other categories, just follow the same procedure.
	sigGO.Term <- getGOTerm(sigGO.ID)[["BP"]]
}
@


% Produce the significant GO term table
<<echo=FALSE, results=tex, eval=F>>=
if (require(GOstats) & require(lumiHumanV1) & require(xtable)) {

	##get gene counts at each GO category
	gg.counts <- geneCounts(hgOver)[sigGO.ID]
	total.counts <- universeCounts(hgOver)[sigGO.ID]

	ggt <- unlist(sigGO.Term)
	numCh <- nchar(ggt)
	ggt2 <- substr(ggt, 1, 17)
	ggt3 <- paste(ggt2, ifelse(numCh > 17, "...", ""), sep="")
	
	## output the significant GO categories as a table
	ggMat <- matrix(c(names(sigGO.Term), ggt3, signif(gGhyp.pv[sigGO.ID],5), gg.counts, total.counts),
    		byrow=FALSE, nc=5, dimnames=list(1:length(sigGO.Term), c("GO ID",
   		"Term", "p-value","Significant Genes No.", "Total Genes No.")))
	xtable.matrix(ggMat,
  		caption="GO terms, p-values and counts.", label="ta:GOggterms")
}
@

% Produce the significant GO term table
% latex table generated in R 2.4.0 by xtable 1.4-2 package
% Sat Dec  9 23:14:22 2006
\begin{table}[ht]
\begin{center}
\begin{tabular}{rlllll}
  \hline
 & GO ID & Term & p-value & Significant Genes No. & Total Genes No. \\
  \hline
1 & GO:0009611 & response to wound... & 8.4244e-06 & 42 & 443 \\
  2 & GO:0006955 & immune response & 8.8296e-06 & 68 & 859 \\
  3 & GO:0006952 & defense response & 1.7525e-05 & 72 & 945 \\
  4 & GO:0006950 & response to stres... & 1.9132e-05 & 81 & 1103 \\
  5 & GO:0009607 & response to bioti... & 5.0811e-05 & 72 & 976 \\
  6 & GO:0009613 & response to pest,... & 7.2813e-05 & 45 & 533 \\
  7 & GO:0006954 & inflammatory resp... & 0.00025402 & 25 & 250 \\
  8 & GO:0009605 & response to exter... & 0.00026005 & 46 & 580 \\
  9 & GO:0051707 & response to other... & 0.00040553 & 45 & 575 \\
  10 & GO:0051674 & localization of c... & 0.00082563 & 30 & 348 \\
  11 & GO:0006928 & cell motility & 0.00082563 & 30 & 348 \\
  12 & GO:0040011 & locomotion & 0.00099205 & 30 & 352 \\
   \hline
\end{tabular}
\caption{GO terms, p-values and counts.}
\label{ta:GOggterms}
\end{center}
\end{table}


\section{Reference}

1. Lin, S.M., Du, P., Kibbe, W.A.,  "Model-based Variance-stabilizing Transformation for Illumina Mi-croarray Data", submitted

2. Du, P., Kibbe, W.A. and Lin, S.M., "nuID: A universal naming schema of oligonucleotides for Illumina, Affymetrix, and other microarrays", submitted.

3. Barnes, M., Freudenberg, J., Thompson, S., Aronow, B. and Pav-lidis, P. (2005) "Experimental comparison and cross-validation of the Affymetrix and Illumina gene expression analysis platforms", Nucleic Acids Res, 33, 5914-5923.

%\bibliographystyle{plainnat}
%\bibliography{lumi}


\end{document}
