%\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{\file}[1]{{\texttt{#1}}}

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

\bibliographystyle{plainnat} 

\begin{document}
\title{Handling probe sequence information: the package 
matchprobes (draft)}
\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)
@

%------------------------------------------------------------
\section{Using probe packages}\label{sec.prbpkg}
%------------------------------------------------------------
The probe sequence data for the array type HG-U95Av2 
(an Affymetrix GeneChip) can be loaded with the command
\begin{Schunk}
\begin{Sinput}
> library(hgu95av2probe)
> help(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:} 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 be to drop the option altogether, since we can
simply use \Rfunction{complementSeq(query, start=13, stop=13)}
on the query string(s) and then call \Rfunction{matchprobes} 
with standard exact matching.

%------------------------------------------------------------
\subsection{C-, G-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{Combine}\label{sec.combine}
%------------------------------------------------------------
Load the data
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% FIXME: provide two CEL files in the data directory and have 
% this example work
% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
\begin{Schunk}
\begin{Sinput}
x1 <- read.affybatch("118T1.cel")
x2 <- read.affybatch("CL2001032020AA.CEL")
\end{Sinput}
\end{Schunk}

Combine the data. For this to work it is required that the packages 
\Rpackage{hu6800probe} and \Rpackage{hgu95av2probe} are installed.
\begin{Schunk}
\begin{Sinput}
res = combine(list(x1, x2), c("hu6800probe","hgu95av2probe"))
\end{Sinput}
\end{Schunk}

The function returns a list \Robject{res} with two elements:
an \Robject{AffyBatch} \Robject{data} and a
\Robject{CDF enviroment} \Robject{newcdf}.
\begin{Schunk}
\begin{Sinput}
newcdf <- res$newcdf
z <- rma(res$data)
\end{Sinput}
\end{Schunk}

View the results.
\begin{Schunk}
\begin{Sinput}
x11(width=7, height=7)
par(mfrow=c(2,2))

plot(exprs(res$data), main="after combine", pch=".", log="xy")
plot(exprs(z),        main="after RMA", pch=".")

## the number of probes in each "new" probe set
prs = multiget(ls(newcdf), newcdf)
nrprobes = sapply(prs, function(x) nrow(x))
barplot(table(nrprobes), xlab="no probes in probeset", ylab="frequency")
\end{Sinput}
\end{Schunk}

%------------------------------------------------------------
\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 data files that can be obtained 
from the vendor (at http://www.netaffx.com). As an example, the file
\file{HG-U95Av2\_probe\_tab.gz} is provided in the \file{data} subdirectory
of the \Rpackage{matchprobes} package.

<<makeprobepackage>>=
filename <- file.path(.path.package("matchprobes"), "data", "HG-U95Av2_probe_tab.gz")
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}
%------------------------------------------------------------
This requires writing your own \Rfunction{importfun}. 

\emph{Input of an importfun:}
\begin{itemize}
\item directory name of a package template directory, containing at 
least 
a file \file{DESCRIPTION}, 
a directory \file{man} with a file @PKGNAME@.Rd, 
a directory \file{data}, and possible other directories and files,
conforming to the usual R package conventions.
\end{itemize}

Output of an importfun:
\begin{itemize}
\item an environment with a dataframe named @PKGNAME@ and a column
`sequence'. The dataframe may have other columns, and the environment 
may contain other data objects. These will be saved as individual 
.rda files of the same name into the data directory of the produced 
package. Documentation needs to be provided for the columns, and for 
the other objects.
\item Symbols:
  \begin{itemize}
  \item ARRAYTYPE: name of the array
  \item DATASOURCE: a textual description of how the data were obtained. 
     It should contain the URL, or the name of the company / person.
  \item FORMAT
  \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}.


\end{document}

