% 
%\VignetteIndexEntry{beadarray image processing} 
%\VignetteKeywords{beadarray image processing} 
%\VignettePackage{beadarray} 
\documentclass{article}



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

\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}}}



\newcommand{\classdef}[1]{
  {\em #1}
}
\begin{document}
\title{Image Analysis with \Rpackage{beadarray}}
\author{Mike Smith}
\maketitle

\section*{Introduction}

<<options, echo=FALSE, eval=TRUE>>==
options(width = 70)
@

<<prelim, echo=FALSE, eval=TRUE>>==
library(beadarray)
@

\noindent From version 2.0 \Rpackage{beadarray} provides more flexibility in the processing of array images and the extraction of bead intensities than its predecessor.  In the past intensity extraction from array images by \Rpackage{beadarray} attmepted to emulated that performed by Illumina, with minimal opportunities for deviation from this.  Whilst the default approach taken in \Rpackage{beadarray} is still to emulated Illumina, we have made each step modular, in order to allow greater flexibility for the user.  This vignette is designed to show how one can read the TIFF images from the BeadArray scanner and implement alternative feature intensity extraction algorithms.

\section{Reading bead-level data into beadarray}

\subsection{Standard Illumina Image Processing}

\noindent The first step in a pipeline for image processing is to read both the TIFF image and the bead-level text file.  The text file contains the identities of each bead, as well as the bead-centre coordinates.  The image processing methods contained within \Rpackage{beadarray} use these coordinates as seed points for intensity extraction, although one can conceive of approaches where bead-centres are calculated seperately, prior to intensity extraction.  However, even with such an approach the Probe ID for each bead will need to be extracted from the \textit{.txt} file.

<<readingTiff, eval=FALSE>>=
tiff <- readTIFF();
data <- readBeadLevelTextFile();
@

\noindent The standard method employed by Illumina's scanner for calculating bead intensity is a four step process, described in Kuhn \textit{et al} \cite{Kuhn2004}.  It can be summarized as:

\begin{itemize}
 \item Calculate background value
 \item Sharpen image
 \item Calculate foreground value
 \item Subtract background for foreground to give final intensity
\end{itemize}

\noindent If the function \Rfunction{readIllumina()} is called with \Rfunarg{useImages = TRUE} then intensities are extracted using code that gives a very close emulation of that used by Illumina.  If one wished to perform this calcuation themselves (outside of \Rfunction{readIllumina()}), it can be done using the following code.

<<standardProcessing, eval=FALSE>>=
bg <- illuminaBackground(tiff, data[,3:4]);
tiffSharp <- illuminaSharpen(tiff);
fg <- illuminaForeground(tiffSharp, data[,3:4]);
finalIntensity <- fg - bg;
@

\noindent  Each of the functions above take a matrix representing the pixel values from the TIFF image as their first argument.  The background and foreground algorithms additionally take a two column matrix containing the coordinates of the bead centres.  If one wished to calculate intensities for only a subset of the beads then supplying only the appropriate bead-centres in this step would achieve this.\\

\noindent  After calculating intensties they need to be inserted into a \Rclass{beadLevelData} object.  The code below shows how to create a new object and insert intensity values.  However, this approach creates an empty \Rclass{beadLevelData} object, which will be lacking any information except that which the user manually inserts.

<<creatingObject, eval=FALSE>>=
BLData <- new(Class = "beadLevelData");
BLData <- insertBeadData(BLData, array = 1, what = "Grn", data = finalIntensity)
@

\noindent An easier alternative to creating your own \Rclass{beadLevelData} object is to use the function \Rfunction{readIllumina()} to read the data as described in the main vignette.  This ensures that any available data (such as sample IDs, scanner metrics, grid sizes etc.) are read in and stored.  One can then choose to overwrite the values generated by \Rfunction{readIllumina()}, or store alternative intensities alongside them. \\

\noindent  The example below first reads the data using the standard arguments to \Rfunction{readIllumina()}, which will extract the intensities from the \textit{.txt} file.  The second step overwrites those intensities with those we calculated previously (which should be very similar).  The final command creates a new entry in the \Rclass{beadLevelData} object (refered to as `GrnLog'), that stores the log transform of the values we calculated earlier.  In this way the user can store a variety of intensity values if they wish to experiment with alternative forms of background subtraction, gradient removal etc.

<<usingReadIllumina, eval=FALSE>>=
BLData <- readIllumina();
BLData <- insertBeadData(BLData, array = 1, what = "Grn", data = finalIntensity)
BLData <- insertBeadData(BLData, array = 1, what = "GrnLog", data = log2(finalIntensity))
@

\subsection{Alternative Methods}

\noindent  The examples above have focused on applying the same intensity extraction algorithms that are employed by Illumina.  However, one may wish to employ an alternative algorithm to test its performance.  The example below implements an alternative method of calculating the background intensity values, as recommended by Smith \textit{et al} \cite{Smith2010}.  

<<alternativeMethods, eval=FALSE>>=
bg <- medianBackground(tiff, data[,3:4]);
@

\noindent We can then use the new background intensities in the same way as previously, before inserting them into the \Rclass{beadLevelData} object.

\section{Parallel Processing}

\noindent  We have included some support for parallel processing in the functions to perform sharpening of the image and the two background calculation methods.  These can offer some increase in throughput when one is using a single computer to analyse a small number of samples.  However if one is dealing with a large number of arrays then there are probably more efficient mechanisms to achieve speedup, such as reading seperate chips on multiple machines (or \texttt{R} sessions) and combining the data after they have been read.\\

\noindent This multicore support is implemented at the \texttt{C} level using the \textbf{OpenMP} library.  Unfortunately adding support for this generates a warning on Bioconductor, so support needs to be added manually and the package build from source.  The procedure is slightly different for users on Linux and Windows machines.\\

\noindent Linux users should create a file called \texttt{Makevars} in the \texttt{beadarray/src} directory and add the following two lines before building the package from source.

\begin{verbatim}
PKG_CFLAGS=-fopenmp
PKG_LIBS=-lgomp
\end{verbatim}

\noindent Windows users should create a file called \texttt{Makevars.win} in the \texttt{beadarray/src} directory and add the following two lines before building the package from source.

\begin{verbatim}
PKG_CFLAGS=-fopenmp
PKG_LIBS=-lgomp -mthreads -lpthreadGC2
\end{verbatim}

\section{Session Info}

<<sessionInfo>>=
sessionInfo()
@

%\section{References}

%\noindent [1] Kuhn K, Baker SC, Chudin E, Lieu MH, Oeser S, Bennett H, Rigault P, Barker D, McDaniel TK, Chee MS, ``A novel, high-performance random array platform for quantitative gene expression profiling'', \textit{Genome Research}, 2004, 14(11):2347-56.\\


%\noindent [2] Smith ML, Dunning MJ, Tavar\'{e} S and Lynch AG, ``Identification and correction of previously unreported spatial phenomena using raw Illumina BeadArray data'', \textit{BMC Bioinformatics}, 2010, 11:208

\bibliographystyle{unsrt}
\bibliography{beadarray}

\end{document}
