% -*- mode:Latex -*-
%\VignetteIndexEntry{Vignette for the documentation by value package (docval)}
%\VignetteDepends{}
%\VignetteKeywords{add-on preprocessing microarray svn rma}
%\VignettePackage{docval}

%% check whether we are running pdflatex
%%--------------------------------------
\newif\ifpdf
  \ifx\pdfoutput\undefined
  \pdffalse                 % we are not running pdflatex 
\else
  \pdfoutput=1              % we are running pdflatex 
  \pdfcompresslevel=9       % compression level for text and image;
  \pdftrue
\fi

%% choose global options depending on whether 
%% we are running under pdflatex
%%-------------------------------------------        
\ifpdf
  \documentclass[a4paper, 12pt, fleqn]{article}  
  \usepackage[pdftex]{graphicx}
  \usepackage[pdftex]{color}
  \usepackage[pdftex,linktocpage]{hyperref}
  \hypersetup{pdfstartview={FitH},
              pdfpagemode={UseOutlines},
              pdftitle={},
              pdfauthor={Dennis Kostka},
              pdfsubject={}}
  \usepackage[pdftex]{graphicx}
  \usepackage{thumbpdf}
\else
  \documentclass[a4paper, 12pt, fleqn]{article}
  \usepackage[dvips]{graphicx}
  \usepackage[dvips]{color}
  \usepackage{hyperref}
  \usepackage{psfig}
\fi

%% packages
%%---------
\usepackage[a4paper,left=3cm,right=3cm,top=3.5cm,bottom=3.5cm,nohead]{geometry}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{yfonts}
\usepackage{amssymb}
\usepackage{latexsym}  
\usepackage{array}
\usepackage[numbers]{natbib}
\usepackage{color}
\usepackage{sectsty}
\usepackage{colortbl}
%% page style:
%%------------
\setlength{\parskip}{0ex}%{1.5ex}
\setlength{\parindent}{0cm}
%\allsectionsfont{\sffamily}


\definecolor{lg}{rgb}{0.9,0.9,0.9}
\definecolor{dg}{rgb}{0.6,0.6,0.6}


%% some new commands
%%------------------
%% bold symbols (vectors)
\newcommand{\bbet}{\boldsymbol{\beta}}
\newcommand{\bbett}{\boldsymbol{\beta}^{\prime}}
\newcommand{\balph}{\boldsymbol{\alpha}}
\newcommand{\bsig}{\boldsymbol{\sigma}}
\newcommand{\bx}{\boldsymbol{x}}
\newcommand{\by}{\boldsymbol{y}}
\newcommand{\bX}{\boldsymbol{X}}
\newcommand{\bY}{\boldsymbol{Y}}
\newcommand{\bK}{\boldsymbol{K}}
\newcommand{\bmu}{\boldsymbol{\mu}}
\newcommand{\bt}{\boldsymbol{\tau}}
\newcommand{\bvth}{\boldsymbol{\vartheta}}
%% bold symbols (matrices)
\newcommand{\BA}{\mathbf{A}}
\newcommand{\BX}{\mathbf{X}}
\newcommand{\BK}{\mathbf{K}}
\newcommand{\BB}{\mathbf{B}}
\newcommand{\BW}{\mathbf{W}}
\newcommand{\BSig}{\mathbf{\Sigma}}
%% indexed vectors
\newcommand{\bxi}{\boldsymbol{x}^{\,(i)}}
\newcommand{\byi}{\boldsymbol{y}^{\,(i)}}
\newcommand{\bXa}{\boldsymbol{X}}
%% sets and spaces
\newcommand{\mcD}{\mathcal{D}}
\newcommand{\mcK}{\mathcal{K}}
\newcommand{\mcF}{\mathcal{F}}
\newcommand{\mcH}{\mathcal{H}}
\newcommand{\mcL}{\mathcal{L}}
\newcommand{\mcN}{\mathcal{N}}
%% colors
\newcommand{\red}[1]{\textcolor{red}{#1}}
\newcommand{\blue}[1]{\textcolor{blue}{#1}}
\newcommand{\green}[1]{\textcolor{green}{#1}}
\newcommand{\gray}[1]{\textcolor{greay}{#1}}
%% misc
\newcommand{\phead}[1]{{ \flushleft \sf\small\textbf{#1}\quad}}
\newcommand{\mbf}{\mathbf}
\newcommand{\arsinh}{\mbox{arsinh}} 
\newcommand{\Ys}{Y^\star}
\newcommand{\ys}{y^\star}

%% the math operators
%%-------------------
\DeclareMathOperator*{\argmax}{argmax}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator*{\diag}{diag}
\DeclareMathOperator*{\sgn}{sgn}

\begin{document}

%\VignetteIndexEntry{docval}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TITLE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Vignette for the documentation by value package (docval)}

\author{ Dennis Kostka and Rainer Spang \\ 
  \vspace{.01cm}\\
   \small Max Planck Institute for Molecular Genetics \\
   \small Ihnestrasse 63-73, 14195 Berlin, Germany }
   
  \date{} 

\maketitle

\tableofcontents

\newpage



\section{Using the software}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection{Preprocessing}
%%%%%%%%%%%%%%%%%%%%%%%%%%

We proceed in three steps:
\begin{itemize}

\item[$i$)] Read in the data. We have prepared an \texttt{AffyBatch} 
object  called \texttt{abo.w}, featuring ten \texttt{.cel} files:

<<results=hide>>=
library(docval)
data(willenbrock)
@
<<>>=
abo.w
@
\item[$ii$)] Preprocessing using \texttt{rma}
<<results=hide>>=
exs.rma  = wrap.val(abo.w[,-10],method="rma")
scl.rma  = preproc(description(exs.rma))$val
@
\item[$iii$)] Preprocessing using \texttt{vsn}
<<results=hide>>=
exs.vsn  = wrap.val(abo.w[,-10],method="vsn")
scl.vsn  = preproc(description(exs.vsn))$val
@
\end{itemize}

\subsection{Building the classification rule}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A classification rule is derived using \texttt{pamr}. 
The function \texttt{pamr.fil} is a simple wrapper, performing
model selection based on crossvalidation-error.

<<results=hide>>=
 labs    = ((as.numeric(pData(abo.w)$IMMUN == "T")-1/2)*2)[-10]
 sig.vsn =  pamr.fil(exs.vsn,labs,fil=FALSE) 
 sig.rma =  pamr.fil(exs.rma,labs,fil=FALSE)
@

\subsection{Documentation by value}

We document the classification rules by value: The signature, together
with the data-dependent scale information, is stored in a binary format
for later application to external samples:

<<results=hide>>=
sig.byval.rma = list(sig=sig.rma, scl=scl.rma)
sig.byval.vsn = list(sig=sig.vsn, scl=scl.vsn)
@

When publishing the signature, the file \texttt{"sig\_byval\_rma.rdat"} or
\texttt{"sig\_byval\_vsn.rdat"} need to be made available on supplemental
web-pages.


\subsection{Diagnosing an external patient}

We diagnose an external patient (\texttt{external\_patient.CEL.gz}).
Again, we proceed in three steps:
 
\begin{itemize}
\item[$i$)] Read in the data. For convenience, we take the left out patient.
<<results=hide>>=
abo.extrnl = abo.w[,10]

@

\item[$ii$)] Add-on preprocess the data, transforming it to a
study-consistent scale.  The function \texttt{wrap.pag.add} utilizes the
results from the theory section.  In order to do that, data-dependent
information stored with the signature has to be retrieved.
<<results=hide>>=
exs.extrnl.rma = wrap.val.add(abo.extrnl,sig.byval.rma$scl,method="rma")
exs.extrnl.vsn = wrap.val.add(abo.extrnl,sig.byval.vsn$scl,method="vsn")
@
\item[$iii$)] Predict the class labels of the external patient, using
the classifier derived beforehand:
<<>>=
diag.rma = sig.byval.rma$sig(exprs(exs.extrnl.rma))
diag.vsn = sig.byval.vsn$sig(exprs(exs.extrnl.vsn))
@
\end{itemize}


\section{Concepts}

\subsection{Documenting quantile normalization for add-on preprocessing (rma)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Assume we have $p$ probes and $n$ arrays. Let $\bX$ be the $p\times n$
background corrected probe-level expression matrix on log scale. Let
$\mathcal{P}$ be a the permutation sorting the columns of $\bX$ and
$\mathcal{P}^{-1}$ its inverse.  Then the  quantile normalized version
$\tilde{\bX}$ of $\bX$ is obtained via:
\[
\tilde{\bX} = \mathcal{P}^{-1}((\mathcal{P}\bX) \mathbf{1})\quad ,
\]
where $\mathbf{1}$ is a $n\times p$ matrix  with all elements equal
to $1/n$.  Further on, let $\bmu$ be equal to the first  column of
$(\mathcal{P}\bX) \mathbf{1}$ and let $\bx \in \mathbb{R}^p$ be an
external array.  If $\mathcal{P}_{\bx}$ is the permutation sorting
the entries of $\bx$, the add-on-quantile-normalized version of $\bx$
consistent with the study is given by via
\[
\tilde{\bx} = \mathcal{P}^{-1}_{\bx} (\bmu) \quad .
\]
Since $\mathcal{P}_x$ depends on $\bx$ only, quantile normalization is
fully documented by $\bmu$.

 
\subsection{Documenting the variance stabilizing transformation for add-on preprocessing (vsn)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Let $\bX$ be a raw probe-level $p \times n $ expression
matrix   of $p$ probes and $n$ samples. The model of Huber et
al.\cite{huber2003parameter,newton2001differential} relates a random
variable $X_{ki}$ ($k=1\ldots p$ and $i=1\ldots n$) stochastically
with the true abundance $\mu_k$ for probe $k$, given probe $k$ is not
differentially expressed:
\begin{equation}
  \label{eq:model}
  \arsinh (a_i + X_{ki}b_i) =: h_i(X_{ki}) =  \mu_k +\epsilon_{ki},\quad \epsilon_{k_i} \sim N(0,c^2) \quad .
\end{equation}
Here $a_i \in \mathbf{R}$, $b_i \in \mathbf{R_{+}}$ and $c\in
\mathbf{R_+}$ are unknown parameters. \citet{huber2002variance} explain
how to estimate these parameters, together with $\mu_k$ and $c$  from
the data at hand. In the main paper, this method is referred to as the
\texttt{vsn} method.

Assume  \texttt{vsn} normalized internal data is at hand.  That is, for
$n$ arrays we have normalized expression values $\{\hat{h}_i(x_{ki})\}$,
with $i=1 \ldots n$ and $k = 1\ldots p$,  and  corresponding  parameter
estimates $\{(\hat{a}_{i},\hat{b}_i)\}$. Further on, a set $\mathcal{K}$
of not differentially expressed genes has been specified.  We also have
the estimates of the $\hat{\mu}_k$ for each gene  $k \in \mathcal{K}$,
as well as estimates of the variance of the residuals in Equation
(\ref{eq:model}):
\begin{equation*}
  \label{eq:plugin}
\begin{aligned}
  &\hat{\mu}_k &=&\quad  \frac{1}{n}\sum_{i=1}^n \hat{h}_i (x_{ki}) \quad \mbox{and} \quad
  \hat{c}^2   &=&\quad  \frac{1}{n|\mathcal{K}|}\sum_{k \in \mathcal{K}}\sum_{i=1}^n (\hat{h}_i(x_{ki} )- \hat{\mu}_k)^2 \quad .
\end{aligned}
\end{equation*} 

Given an external  $(n+1)$-th sample $\{x_k^\star\}_{k=1}^p$, we want to
transform it to the scale determined by the $n$ core arrays.  By employing
the same stochastic model as for the core data (Equation (\ref{eq:model}))
and plugging in available estimates, we get the following model for the
new sample:
\begin{equation*}
  \label{eq:newmodel} 
\arsinh (a_{n+1}+ b_{n+1}X^\star_k) = \hat{\mu}_k +\epsilon_{k}, \quad \epsilon_k \sim N(0,\hat{c}^2) \;\;\; \mbox{for} \;\;\; k \in \mathcal{K}. 
\end{equation*}

Maximum likelihood estimators for the parameters $a_{n+1}$ and $b_{n+1}$
are available as
\begin{equation*}
  \label{eq:maxlikeest}
\begin{aligned}
& (\hat{a}_{n+1}, \hat{b}_{n+1})   &=&\quad  \argmin_{(a,b)}  &&\sum_{k \in \mathcal{K}} \frac{( h(x^\star_k)-\hat{\mu}_k )^2}{2\hat{c}^2}  - \sum_{k\in \mathcal{K}}\log ( \partial_{x^\star_k}h(x^\star_k)) \quad 
\end{aligned}
\end{equation*}
and can be calculated numerically. The estimates are completely determined
by the $\hat{\mu}_k$, $\hat{c}^2$ and measurement values on the external
array. So is the variance stabilizing transformation $h(x^\star_k)
= h(\hat{a}_{n+1} + \hat{b}_{n+1}x^\star_k)$ bringing the $(n+1)$-th
sample to the core scale.  Therefore \texttt{vsn} normalization is fully
documented by $\{\hat{\mu}_k\}_{k=1}^{|\mathcal{K}|}$ and $\hat{c}$.



\subsection{Documenting an additive model of probeset summary for add-on preprocessing (rma,vsn)}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Let $\bX$ be the $p\times n$ background corrected and normalized
probe-level expression matrix on log scale. Let $\bX ^{(k)}$ be the
submatrix indexed by the probes belonging to  probe set $k$ across
all arrays.  Then an additive model assumes
\[
\bX ^{(k)}_{ij}  = p_i + g_j + \epsilon \quad,
\]
where $p_i$ is a probe specific effect and $g_j$ represents the
abundance of \texttt{mRNA} of gene $k$ in array $j$. $\hat{p}_i$
and $\hat{g}_j$ can be estimated  by a median polish procedure
\cite{tukey1977exploratory}.  In that case, $\mbox{\small \texttt{median}}
(\bx^{(k)} - \boldsymbol{\hat{p}}^{(k)})$ denotes a study-consistent
estimate of the expression of gene $k$, given the external array $\bx$.
$\boldsymbol{\hat{p}}^{(k)}$ is the vector of estimated probe effects
associated with  probeset $k$, and $\bx^{(k)}$ is the vector of normalized
expression values of the same probeset.  That is, the additive model
is fully documented by keeping track of the  probe effects $\hat{p}_i$
for all probes $i$ on the array.



\begin{thebibliography}{1}

\bibitem{huber2003parameter}
  W Huber, A von Heydebreck, H Sültmann, A Poustka and M Vingron. ``Parameter estimation for the calibration and variance stabilization of microarray data.'', \emph{Statistical Applications in Genetics and Molecular Biology}, 2(1):Art 3, 2003. ISSN 1544--6115.

\bibitem{newton2001differential}
  MA Newton, CM Kendziorski, CS Richmond, FR Blattner and KW Tsui. ``On differential variability of expression ratios: improving statistical inference about gene expression changes from microarray data.'', \emph{J Comput Biol}, 8(1):37-52, 2001. URL \texttt{http://dx.oi.org/10.1089/106652701300099074}.

\bibitem{huber2002variance}
  W Huber, A von Heydebreck, H Sültmann, A Poustka and M Vingron. ``Variance stabilization applied to microarray data calibration and to the quantification of differential expression.'', \emph{Bioinformatics}, 18 Suppl 1:S96-104, 2002.

\bibitem{tukey1977exploratory}
  JW Tukey. \emph{Exploratory Data Analysis}. Addison-Wesley, 1977.



\end{thebibliography}


%\bibliographystyle{plainnat}
%\bibliography{/home/kostka/tex/bib/pag_norm}




\end{document}



