%\VignetteIndexEntry{Handling probe sequence information}
%\VignetteDepends{Biobase, matchprobes, affy, hgu95av2cdf}
%\VignetteKeywords{Expression Analysis}
%\VignettePackage{matchprobes}

\documentclass{article}
\usepackage{times}
\usepackage{hyperref}

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

\parskip=3mm
\oddsidemargin=3mm
\evensidemargin=3mm
\headheight=-2cm
\parindent0cm

\bibliographystyle{plainnat} 

\begin{document}
\title{Handling probe sequence information: the package 
matchprobes}
\author{Wolfgang Huber and Robert Gentleman}
\maketitle
\tableofcontents

%------------------------------------------------------------
\section{Overview}
%------------------------------------------------------------
This package allows to calculate with the sequences and related
information of the probes on a microarray. The sequences themselves
are stored in separate packages, one for each type of microarray. For
some commonly used types, the packages can be downloaded from
http://www.bioconductor.org. You can also create such packages
yourself. This is described in section~\ref{sec.mpp}.

<<startup,results=hide>>=
library(matchprobes)
library(affy)
library(hgu95av2cdf)
library(hgu95av2probe)
library(hu6800cdf)
library(hu6800probe)
@

%------------------------------------------------------------
\section{Using probe packages}\label{sec.prbpkg}
%------------------------------------------------------------
Help for the probe sequence data packages can be accessed through
\begin{Schunk}
\begin{Sinput}
> ? hgu95av2probe
\end{Sinput}
\end{Schunk}

%------------------------------------------------------------
\subsection{Functions}\label{subsec.fcts}
%------------------------------------------------------------
Complementary sequence
<<complementSeq>>=
example(complementSeq)
@

Reverse sequence
<<reverseSeq>>=
example(reverseSeq)
@

Matching probes
<<matchprobes>>=
pm    <- hgu95av2probe$sequence
query <- hgu95av2probe$sequence[21:23]
m     <- matchprobes(query, pm)
unlist(m)
@

%------------------------------------------------------------
\subsection{Mismatch sequences for Affymetrix GeneChips}
\label{subsec.mmseq}
%------------------------------------------------------------
The mismatch sequences are not explicitely stored in the 
probe packages. They are implied by the match sequences, by 
flipping the middle base. This can be done with the 
\Rfunction{complementSeq}. 
<<mismatchSeq>>=
mm <- complementSeq(hgu95av2probe$sequence, start=13, stop=13)
cat(pm[1], mm[1], sep="\n")
@
For Affymetrix GeneChips the 
length of the probe sequences is 25, and since we start 
counting at 1, the middle position is 13.

The function \Rfunction{matchprobes} also responds to this 
type of implied mismatch sequences, by multipying the matching 
indices with -1:

<<matchprobes-mismatch>>=
m     <- matchprobes(query, c(pm, mm))
unlist(m)
@

{\bf Note}(wh): This may need to be set as an option, rather than 
as a standard behavior - I don't particularly like that the 
position 13 is hardcoded into \Rfunction{matchprobes}. 
Best may in fact be to drop the mismatch matching altogether, 
since we could simply use \Rfunction{complementSeq(query, start=13, stop=13)}
on the query string(s) and then call \Rfunction{matchprobes} 
with standard exact matching.

%------------------------------------------------------------
\subsection{Base content}
\label{subsec.basecontent}
%------------------------------------------------------------
The base content of the probes on an array can be obtained
through the function \Rfunction{basecontent}.

<<basecontent>>=
bcpm <- basecontent(hgu95av2probe$sequence)
as.matrix(bcpm[1:5,])
bcmm <- basecontent(complementSeq(hgu95av2probe$sequence, start=13, stop=13))
as.matrix(bcmm[1:5,])
@

%------------------------------------------------------------
\section{CombineAffyBatch}\label{sec.combine}
%------------------------------------------------------------
Load the data. 
<<loadcel>>=
f1 <- system.file("extdata", "118T1.cel.gz", package="matchprobes")
f2 <- system.file("extdata", "CL2001032020AA.cel.gz", package="matchprobes")
pd1 <- new("phenoData", pData=data.frame(id="pi"), varLabels=list("phenovar"))
pd2 <- new("phenoData", pData=data.frame(id="bh"), varLabels=list("phenovar"))
x1 <- read.affybatch(filenames=f1, compress=TRUE, phenoData=pd1)
x2 <- read.affybatch(filenames=f2, compress=TRUE, phenoData=pd2)
@
\Robject{f1} and \Robject{f2} contain the filenames of the CEL files. 
In practice, they will be vectors with many filenames per array type, here, 
for demonstration, we only use one single CEL file per array type, 
which are provided in the subdirectory \Rfile{extdata} of the package.

Combine the data. For this to work it is required that the packages 
\Rpackage{hu6800probe} and \Rpackage{hgu95av2probe} are installed.
<<docombine>>=
res <- combineAffyBatch(list(x1, x2), c("hu6800probe","hgu95av2probe"), newcdf="comb")
## hu6800probe
@

The function returns a list \Robject{res} with two elements:
an \Robject{AffyBatch} \Robject{data} and a
\Robject{CDF enviroment} \Robject{newcdf}.
<<dorma>>=
comb <- res$cdf
z    <- rma(res$dat)
@

View the results.
<<viewres,fig=TRUE>>=
par(mfrow=c(2,2))
plot(exprs(res$dat), main="after combine", pch=".", xlab="f1", ylab="f2", log="xy")
plot(exprs(z),       main="after RMA",     pch=".", xlab="f1", ylab="f2")
## the number of probes in each "new" probe set
prs <- multiget(ls(comb), comb)
nrprobes <- sapply(prs, function(x) nrow(x))
barplot(as.vector(table(nrprobes)), 
   xlab="number of probes per probeset", 
   ylab="frequency")
@

%------------------------------------------------------------
\section{Creating probe packages}\label{sec.mpp}
%------------------------------------------------------------
Probe packages are a convenient way for distributing and storing the 
probe sequences and related information. 

%------------------------------------------------------------
\subsection{Affymetrix genechips}\label{subsec.affy}
%------------------------------------------------------------
In this section we see how a probe package can be created for
Affymetrix genechips from the tabulator-separated sequence files that
can be obtained from the vendor (at
\url{http://www.affymetrix.com/support}). As an example, the file
\Rfile{HG-U95Av2\_probe\_tab.gz} is provided in the \Rfile{extdata}
subdirectory of the \Rpackage{matchprobes} package.

<<makeprobepackage>>=
filename <- system.file("extdata", "HG-U95Av2_probe_tab.gz", package="matchprobes")
outdir   <- tempdir()
me       <- "Wolfgang Huber <w.huber@dkfz.de>"
makeProbePackage("HG-U95Av2",
                 datafile   = gzfile(filename, open="r"),
                 outdir     = outdir,
                 maintainer = me, 
                 version    = "0.0.1", 
                 force      = TRUE)
dir(outdir)
@

%------------------------------------------------------------
\subsection{Other chiptypes}\label{subsec.otherarraytypes}
%------------------------------------------------------------ 
To deal with different file formats and additional types of probe
annotation data from public or in-house databases, the function
\Rfunction{makeProbePackage} offers a great deal of flexibility. The
user can specify her own import function through the argument
\Robject{importfun}. By default, its value is
\Rfunction{getProbeDataAffy}, a function that reads tabular Affymetrix
genechip sequence files. Import functions for other types of arrays
can be adapted from this prototype.

The help pages and R code contained in the produced packages are
derived from a template directory that obeys the usual R package
conventions~\cite{RExt}. A prototype for such a directory is provided
within the package \Rpackage{matchprobes}.  To facilitate the
automated production of large numbers of similar packages, we provide
a text substitution mechanism similar to the one used in the GNU
\Rfunction{configure} system.

The input parameters of an import function are
\begin{itemize}
\item a character string naming the array type
\item the input data files
\item a character string with the directory name of a package 
  template directory, containing at least a file \Rfile{DESCRIPTION}, 
  a directory \Rfile{man} with a file @PKGNAME@.Rd, 
  a directory \Rfile{data}, and possible other directories and files,
  conforming to the usual R package conventions.
\item ... any sort of further parameters, as necessary.
\end{itemize}

The output of an import function is a named list with elements
\begin{itemize}
\item \texttt{pkgname}: a character string with the package name.
\item \texttt{dataEnv}: an environment containing an arbitrary number 
  of data objects; these make the core of the package. Among these,
  there should be a data frame whose name is the value of
  \texttt{pkgname} with a column \Robject{sequence}.
  This data frame may have other columns such as the $x$- and $y$-position of the probes
  on the array, identifiers linking a probe to genomic databases, or its
  length and relative position within the gene it represents.  
  The objects in the environment will be saved as individual 
  \Rfile{.rda} files of the same name into the data directory of the 
  produced package. Documentation needs to be provided for the 
  columns of the data frame, as well as  for the other objects.
\item \texttt{symVals}: a named list, containing the symbol-value 
  substitutions that are used in the text processing. It must at least 
  contain the elements
  \begin{itemize}
  \item \texttt{ARRAYTYPE}: name of the array
  \item \texttt{DATASOURCE}: a textual description of how the data
     were obtained. It should contain the URL, or the name of the
     company / person.
  \end{itemize}
\end{itemize}

For more details, please refer to the help files for the functions 
\Rfunction{makeProbePackage} and \Rfunction{getProbeDataAffy}.
For an example, refer to the source code of \Rfunction{getProbeDataAffy}.

\begin{thebibliography}{10}
\bibitem{RExt}
R Foundation (1999).
\newblock {\em Writing R Extensions}.
\newblock \url{http://www.r-project.org}.
\end{thebibliography}

\end{document}

