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

\usepackage{/homes/madman/R-1.9.0/share/texmf/Sweave}
\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}.

\begin{Schunk}
\begin{Sinput}
> library(matchprobes)
> library(affy)
> library(hgu95av2cdf)
> library(hgu95av2probe)
> library(hu6800cdf)
> library(hu6800probe)
\end{Sinput}
\end{Schunk}

%------------------------------------------------------------
\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
\begin{Schunk}
\begin{Sinput}
> example(complementSeq)
\end{Sinput}
\begin{Soutput}
cmplmS> seq <- c("AAACT", "GGGTT")

cmplmS> complementSeq(seq)
[1] "TTTGA" "CCCAA"

cmplmS> seq <- c("CGACTGAGACCAAGACCTACAACAG", "CCCGCATCATCTTTCCTGTGCTCTT")

cmplmS> complementSeq(seq, start = 13, stop = 13)
[1] "CGACTGAGACCATGACCTACAACAG" "CCCGCATCATCTATCCTGTGCTCTT"
\end{Soutput}
\end{Schunk}

Reverse sequence
\begin{Schunk}
\begin{Sinput}
> example(reverseSeq)
\end{Sinput}
\begin{Soutput}
rvrsSq> w <- c("hey there", "you silly fool")

rvrsSq> reverseSeq(w)
[1] "ereht yeh"      "loof yllis uoy"

rvrsSq> w <- "able was I ere I saw Elba"

rvrsSq> reverseSeq(w)
[1] "ablE was I ere I saw elba"
\end{Soutput}
\end{Schunk}

Matching probes
\begin{Schunk}
\begin{Sinput}
> pm <- hgu95av2probe$sequence
> query <- hgu95av2probe$sequence[21:23]
> m <- matchprobes(query, pm)
> unlist(m)
\end{Sinput}
\begin{Soutput}
match1 match2 match3 
    21     22     23 
\end{Soutput}
\end{Schunk}

%------------------------------------------------------------
\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}. 
\begin{Schunk}
\begin{Sinput}
> mm <- complementSeq(hgu95av2probe$sequence, start = 13, stop = 13)
> cat(pm[1], mm[1], sep = "\n")
\end{Sinput}
\begin{Soutput}
TCTCCTTTGCTGAGGCCTCCAGCTT
TCTCCTTTGCTGTGGCCTCCAGCTT
\end{Soutput}
\end{Schunk}
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:

\begin{Schunk}
\begin{Sinput}
> m <- matchprobes(query, c(pm, mm))
> unlist(m)
\end{Sinput}
\begin{Soutput}
 match1  match2  match3  match4  match5  match6 
     21 -199105      22 -199106      23 -199107 
\end{Soutput}
\end{Schunk}

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

\begin{Schunk}
\begin{Sinput}
> bcpm <- basecontent(hgu95av2probe$sequence)
> as.matrix(bcpm[1:5, ])
\end{Sinput}
\begin{Soutput}
 A T C  G
 2 9 9  5
 6 3 9  7
 5 4 9  7
 5 6 8  6
 5 7 3 10
\end{Soutput}
\begin{Sinput}
> bcmm <- basecontent(complementSeq(hgu95av2probe$sequence, start = 13, 
+     stop = 13))
> as.matrix(bcmm[1:5, ])
\end{Sinput}
\begin{Soutput}
 A  T C G
 1 10 9 5
 7  2 9 7
 4  5 9 7
 5  6 7 7
 5  7 4 9
\end{Soutput}
\end{Schunk}

%------------------------------------------------------------
\section{CombineAffyBatch}\label{sec.combine}
%------------------------------------------------------------
Load the data. 
\begin{Schunk}
\begin{Sinput}
> 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)
\end{Sinput}
\end{Schunk}
\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.
\begin{Schunk}
\begin{Sinput}
> res <- combineAffyBatch(list(x1, x2), c("hu6800probe", "hgu95av2probe"), 
+     newcdf = "comb")
\end{Sinput}
\begin{Soutput}
package:hu6800probe	hu6800probe	
package:hgu95av2probe	hgu95av2probe	
34428 unique probes in common
\end{Soutput}
\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}
> comb <- res$cdf
> z <- rma(res$dat)