%\VignetteIndexEntry{vsn}
%\VignetteDepends{Biobase,vsn,affy,marrayNorm,multtest}
%\VignetteKeywords{Expression Analysis}
%\VignettePackage{vsn}

\documentclass[11pt]{article}
\usepackage[margin=2.5cm,noheadfoot]{geometry}

%------------------------------------------------------------
% newcommand
%------------------------------------------------------------
\newcommand{\code}[1]{{\tt #1}}
\newcommand{\Rfunc}[1]{{\tt #1}}

\newcommand{\myincfig}[3]{%
  \begin{figure}[htbp]
    \begin{center}
      \includegraphics[width=#2]{#1}
      \caption{\label{#1}#3}
    \end{center}
  \end{figure}
}

%------------------------------------------------------------
\begin{document}
%------------------------------------------------------------
\title{Robust calibration and variance stabilization with VSN}
\author{Wolfgang Huber}
\maketitle

<<loadotherpackages,echo=FALSE,results=hide>>=
library(affy)
library(marrayNorm)
library(multtest)
@

%------------------------------------------------------------
\section*{Introduction}
%------------------------------------------------------------
\code{vsn} is a method to preprocess DNA microarray intensity data. As
input, the function \code{vsn} takes the raw intensity measurements
from the DNA probes on a series of microarrays. The intensities from
each array are calibrated by a suitable affine transformation, then
transformed by a variance-stabilizing transformation.  After this,
systematic array- or dye-biases should be removed, and the variance
should be approximately independent of the mean intensity.  This is
useful for subsequent analyses such as hypothesis tests, ANOVA
modeling, clustering, or classification that assume that the variance
is the same for all observations\footnote{% 
  Note that \code{vsn} only addresses the dependence of the variance
  on the mean intensity. There may be other factors influencing the
  variance, such as gene-inherent properties, or changes of the
  tightness of transcriptional control in different conditions.  If
  necessary, these need to be addressed by other methods.}.  
Differences between the transformed values are the
so-called "generalized log-ratios".  If both numerator and denominator
are well above background, generalized log-ratios coincide with the
usual log-ratios:
$h(x_i)-h(x_j)\approx\log(x_i)-\log(x_j)=\log(x_i/x_j)$ if $x_i,
x_j\gg0$~\cite{HuberISMB2002,HSG2002}. In contrast to log-ratios,
they remain well-defined and statistically meaningful if $x_i$ or
$x_j$ are close to zero.

%------------------------------------------------------------
\section{The data} 
%------------------------------------------------------------
The prefered input type for \code{vsn} are objects of class
\code{exprSet}. It is also possible to pass matrices, data frames with
numeric columns only, and objects of class \code{marrayRaw}.
Furthermore, \code{vsn} can be used as a normalization method in the
function \code{expresso} in the package \code{affy}
(cf.~Section~\ref{affy}).

To load intensity data from your own experiments, you can use the
function \code{read.table}, the \code{read}-functions from the package
\code{marrayInput}, or \code{ReadAffy} from the package \code{affy}.

%------------------------------------------------------------
\section{Running vsn on the data from a single cDNA array} 
%------------------------------------------------------------
The package includes example data from a cDNA array on which two
biologically highly similar samples, one labeled in green (Cy3), one in red
(Cy5), were hybridized.
% begin code ----------------------------------------
<<loadpackage,results=hide>>= 
library(vsn) 
data(kidney) 
@ 
% end code -----------------------------------------
The two columns of the matrix \code{exprs(kidney)} contain the
green and red intensities, respectively. Let's try out \code{vsn} on
these example data. In Fig.~\ref{vsn-nkid-scp} you can see the
scatterplot of the calibrated and transformed data. For comparison,
the scatterplot of the log-transformed raw intensities is also shown.
% begin code ------------------------------------------------------
<<nkid-scp,include=FALSE,fig=TRUE,width=8,height=4>>=
par(mfrow=c(1,2))
nkid <- vsn(kidney)
plot(exprs(nkid), main = "vsn", pch = ".")
plot(log.na(exprs(kidney)), main = "raw", pch = ".")
@
% end code ------------------------------------------------------------
\myincfig{vsn-nkid-scp}{\textwidth}{Scatterplots of the kidney 
example data}
\code{vsn} returns the transformed intensities in an object of class
\code{exprSet}. Its slot \code{exprs} is a matrix of the same size as
the input data.  The plot in Fig.~\ref{vsn-nkid-scp} shows the complete
set of $n=9216$ red and green intensities, without any thresholding or
masking of data points.  To verify the variance stabilization, there
is the function \code{meanSdPlot}. For each probe $k=1,\ldots,n$
it shows the estimated standard deviation $\hat{\sigma}_k$ on the
$y$-axis versus the rank of the average $\hat{\mu}_k$ on the $x$-axis,
\begin{equation}
\hat{\mu}_k     =\frac{1}{d}  \sum_{i=1}^d h_{ki}\quad\quad
\hat{\sigma}_k^2=\frac{1}{d-1}\sum_{i=1}^d (h_{ki}-\hat{\mu}_k)^2.
\end{equation}
% begin code ------------------------------------------------------------
<<nkid-sdmean,include=FALSE,fig=TRUE,width=8,height=4,results=hide>>=
par(mfrow=c(1,2))
meanSdPlot(nkid, ranks=TRUE)
meanSdPlot(nkid, ranks=FALSE)
@
\myincfig{vsn-nkid-sdmean}{\textwidth}{Standard deviation versus 
rank of the mean, and the mean, respectively}
% end code ------------------------------------------------------------
Such a plot is shown in Fig.~\ref{vsn-nkid-sdmean}.
The red dots, connected by lines, show the running median of the standard 
deviation\footnote{Window width: 10\%, window midpoints 5\%, 10\%, 15\%, \ldots.}.
Within each window, the median may be considered a pooled estimator of the standard 
deviation, and the curve given by the red line is an estimate
of the systematic dependence of the standard deviation on the mean. After 
variance stabilization, this should be approximately a horizontal line.
It may have some random fluctuations, but should not show an overall trend. 
If this is not the case, that could indicate a data quality 
problem, see Section~\ref{sec.qc}. The rank ordering distributes the data 
evenly along the $x$-axis, which is in many cases useful for the visualization. 
A plot in which the $x$-axis shows the average intensities themselves is 
obtained by calling the \code{plot} command with the argument 
\code{ranks=FALSE}.

The parameter estimation in \code{vsn} works in an iterative manner.
To verify that the iterations have converged, you can call the
function \code{vsnPlotPar}.
% begin code -------------------------------------------------------
<<nkid-iter,include=FALSE,fig=TRUE,width=8,height=4,results=hide>>=
par(mfrow=c(1,2))
vsnPlotPar(nkid, "offsets")
vsnPlotPar(nkid, "factors")
@
% end code ---------------------------------------------------------
\myincfig{vsn-nkid-iter}{\textwidth}{Iteration trajectory of the 
calibration and transformation parameters}
The plots in Fig.~\ref{vsn-nkid-iter} show the values of the estimated
calibration and variance stabilization parameters on the $y$-axis as a
function of the iteration index. All curves should reach a plateau
well before the last iteration. If this is not the case, the number of
iterations may be increased through the parameter \code{iter}, and/or
\code{lts.quantile} may be decreased. It could also indicate a
data quality problem, see Section~\ref{sec.qc}.

The "generalized log-ratios"~\cite{HuberISMB2002} for this experiment
may be obtained for further processing through
% begin code -------------------------------------------------------
<<nkid-histM,include=FALSE,fig=TRUE,results=hide>>=
M <- exprs(nkid)[,2] - exprs(nkid)[,1] 
hist(M, breaks=50, col="#d95f0e")
@
% end code -----------------------------------------------------------
\myincfig{vsn-nkid-histM}{0.4\textwidth}{Histogram of generalized 
log-ratios for the kidney example data}
The histogram is shown in Fig.~\ref{vsn-nkid-iter}.

%---------------------------------------------------------------------
\section{Calibration} 
%---------------------------------------------------------------------
We can access the transformation and calibration parameters through
% begin code ---------------------------------------------------------
<<opts,echo=FALSE,results=hide>>=
oopt <- options(digits=3)
on.exit(options(oopt))
@
<<nkid-calib1>>=
prep <- preproc(description(nkid))
names(prep)
prep$vsnParams
@
% end code -----------------------------------------------------------
The \code{description} slot of an \code{exprSet} is an object of class
\code{MIAME}, and may contain annotation (or ``metadata'') pertinent
to the experiment represented by the object. For an exprSet with $d$
columns, \code{prep\$vsnParams} is a numeric vector of length $2d$. Its
elements $1,\ldots,d$ contain the additive calibration and
transformation parameters $a_i$, its elements $d+1,\ldots,2d$ the
multiplicative ones. Compare the above numbers with the final values
in Fig.~\ref{vsn-nkid-iter}.

If $y_{ki}$ is the matrix of uncalibrated data, 
with $k$ indexing the rows and $i$ the columns, then the calibrated 
data $y_{ki}'$ is obtained through scaling by $b_i$ and shifting 
by $a_i$:
\begin{equation}\nonumber
y_{ki}' = a_i+b_i y_{ki}
\end{equation}
Now suppose the kidney example data were not that well measured, and 
the red channel had a baseline that was shifted by 500 and a scale 
that differed by a factor of $0.25$:
<<nkid-calib2>>=
bkid <- kidney
exprs(bkid)[,"red"] <- 0.25 * (500+exprs(bkid)[,"red"])
@
We can again call \code{vsn} on this data
<<nkid-calib-warn1,echo=FALSE,results=hide>>=
options(warn=-1)
@
<<nkid-calib3,results=hide>>=
nbkid <- vsn(bkid)
<<nkid-calib4,fig=TRUE,include=FALSE,width=8,height=4>>=
par(mfrow=c(1,2))
plot(exprs(bkid),  main = "raw", pch = ".", log="xy")
plot(exprs(nbkid), main = "vsn", pch = ".")
preproc(description(nbkid))$vsnParams
@
\myincfig{vsn-nkid-calib4}{\textwidth}{Scatterplots for badly biased 
data. Left hand side: raw data on log-log scale, right hand side: after 
calibration and transformation with vsn.}
<<nkid-calib-warn2,echo=FALSE,results=hide>>=
options(warn=0)
@
The factor for the red channel is now about four times as large as
before.  The result is shown in Fig.~\ref{vsn-nkid-calib4}.

%------------------------------------------------------------
\section{Running vsn on the data from multiple cDNA arrays} 
%------------------------------------------------------------
The package includes example data from a series of 8 cDNA arrays 
on which different lymphoma were hybridized together with a reference 
cDNA~\cite{Alizadeh}.
<<lymphoma>>=
data(lymphoma)
pData(lymphoma)
@
The 16 columns of the \code{lymphoma} object contain the red 
and green intensities, respectively, from the $8$ slides, as shown 
in the table. Thus, the CH1 intensities are in columns 
$1, 3, \ldots, 15$, the CH2 intensities in columns $2, 4, \ldots, 16$.

There are now two modes of operation for \code{vsn}: First,
\code{vsn} can be called for each slide in turn. Second,
if there is reason to believe that the slides are 
of similar quality, it can also be called on all of them at once:
%------------------------------------------------------------
<<lym-sdmean1,results=hide>>=
lym <- vsn(lymphoma, verbose=FALSE)
@
<<lym-sdmean2,include=FALSE,fig=TRUE>>=
meanSdPlot(lym)
@
\myincfig{vsn-lym-sdmean2}{0.5\textwidth}{Standard deviation 
versus rank of the mean for the lymphoma example data}
%------------------------------------------------------------
This calculation may take a while. Again, Fig.~\ref{vsn-lym-sdmean2} 
helps to visually verify that the variance stabilization worked.
As above, we can obtain the "generalized log-ratios" for each slide,
by subtracting the common reference intensities from those for the 
8 samples:
%------------------------------------------------------------
<<lym-M,include=FALSE,fig=TRUE,width=8,height=4,results=hide>>=
refrs <- (1:8)*2-1
samps <- (1:8)*2
M <- exprs(lym)[,samps] - exprs(lym)[,refrs] 
colnames(M) <-  pData(lymphoma)[samps, "sample"]
A <- rowMeans(exprs(lym))
par(mfrow=c(1,2))
plot(A, M[,"CLL-13"], pch=".")
abline(h=0, col="red")
plot(A, M[,"DLCL-0032"], pch=".")
abline(h=0, col="red")
@
\myincfig{vsn-lym-M}{\textwidth}{Mean-difference plots for two slides 
from the lymphoma example data}
%------------------------------------------------------------
Fig.~\ref{vsn-lym-M} shows the analagon to the $M$-vs-$A$-plots
as described in reference~\cite{Dudoit578}. Note that in the left
scatterplot, there is a cloud of points at low intensities that is
concentrated slightly off the line $M=0$.  In the right scatterplot, a
similar cloud sits right on the $M=0$ line. This could be related to a
quality problem with the left slide (e.\,g. related to the PCR
amplification or the printing, see Section~\ref{sec.qc}).

%-------------------------------------------------------------------
\section{Comparing calibration and data transformation methods}
%-------------------------------------------------------------------
To compare different microarray calibration and data transformation
methods one needs to specify a measure of goodness. One approach is to
compare the obtained values against a known truth. This can be done
in controlled spike-in experiments and in dilution series, which
allow to systematically assess the performance of the methods at 
different biologically relevant spike-in concentrations. Like
any statistical method, the methods may make different choices with 
respect to the trade-off between bias and variance.

Here, we just consider a much coarser and simpler criterion:
the overall sensitivity and specificity in detecting differential 
transcription of the combined procedures of data preprocessing and  
statistical hypothesis testing. Such an analysis can generally be 
applied to any data set that contains replicated measurements of 
samples from biologically clearly distinct, known 
groups (e.\,g.~\cite{HuberISMB2002}).

We use the $t$-test to find genes that are differentially transcribed 
between CLL and DLCL. The comparison relies on the idea that for 
a given selection of genes we can estimate the {\em False Discovery 
Rate} through a permutation method~\cite{StoreyTibshirani}.
Thus, if a method consistently produces a smaller false discovery 
rate for the same size gene list than another method, it may be 
considered more specific.

In the following piece of code, data is transformed to log-ratios and 
normalized by the global median. The result, together with that of 
the previous section, is stored in the 3-dimensional array \code{Ms}. 
Fig.~\ref{vsn-Ms-plot} shows, for one of the arrays, a comparison 
of the resulting log-ratios.
<<Ms-plot,results=hide,fig=TRUE,include=FALSE>>=
library(marrayNorm)
Ms <- array(0, dim=c(dim(M), 2))
Ms[,,1] <- M
mr <- new('marrayRaw', maGf=exprs(lymphoma)[,refrs], maRf=exprs(lymphoma)[,samps])
mn <- maNorm(mr, norm="median", echo=T)
Ms[,,2] <- mn@maM
plot(Ms[,1,], xlab="vsn", ylab="global median", main="M (array 1)", pch=".")
@
\myincfig{vsn-Ms-plot}{0.5\textwidth}{Generalized log-ratio from
\code{vsn} ($x$-axis) versus the log-ratio after global median
normalization ($y$-axis)}

The function \code{calc.fdr} provides a very simple-minded 
implementation of the estimation of the false discovery rate from 
the permutation distribution of the sorted test statistics. 
The result is shown in Fig.~\ref{vsn-fdr3}.
\myincfig{vsn-fdr3}{0.5\textwidth}{False Discovery Rate}
<<fdr1,results=hide>>=
library(multtest)
calc.t <- function(classlabel, M) {
  t <- mt.teststat(M, classlabel) 
  t[abs(t)>1e30] <- NA
  return(t)
}
@
<<fdr2,results=hide>>=
calc.fdr <- function(M, classlabel, nrgeneselect) {
  t <- calc.t(classlabel, M)
  st <- sort(abs(t), decreasing=TRUE)
  threshold <- st[nrgeneselect]

  n <- length(classlabel)
  grp2 <- which(classlabel)
  nck <- nchoosek(n-1, length(grp2)-1)
  permclasslabel <- matrix(0, nrow=n, ncol=ncol(nck))
  for(p in 1:ncol(nck)) {
   permclasslabel[      n, p] <- 1
   permclasslabel[nck[,p], p] <- 1
  }     

  permt <- apply(permclasslabel, 2, calc.t, M)
  fdr <- numeric(length(nrgeneselect))
  for(j in 1:length(threshold)) {
    pnrsel <- apply(permt, 2, function(pt) length(which(abs(pt)>=threshold[j])))
    fdr[j] <- mean(pnrsel) / nrgeneselect[j]
  }
  return(fdr)
}
@
<<fdr3,fig=TRUE,include=FALSE,results=hide>>=
classlabel <- regexpr("CLL",colnames(M)) > 0
nrgeneselect <- c(10, 20, 50, 100, 200, 500, 1000)  
fdr <- apply(Ms, 3, calc.fdr, classlabel, nrgeneselect) 
plot(range(nrgeneselect), range(fdr), type='n', log='x', 
xlab='No. of genes selected', ylab='Estimated FDR')
for(j in 1:ncol(fdr)) 
  lines(nrgeneselect, fdr[,j], type='b', pch=19, lty=j)
legend(10, max(fdr), c("vsn", "global median"), lty=1:ncol(fdr))
@

%---------------------------------------------------------
\section{Running vsn on Affymetrix data} \label{affy}
%---------------------------------------------------------
The package \code{affy} provides excellent functionality for reading and 
processing Affymetrix genechip data. To use \code{vsn} for the 
calibration and transformation of the probe intensities, a wrapper 
for \code{vsn} is provided that can be used within the 
data processing routines of \code{affy}. See the documentation for 
the package \code{affy} for more information about data structures 
and other available methods. There are different ways to treat the
mismatch (MM) intensities, to adjust for background, and summarize 
the probe sets. 

The following code is taken from the example in the man page for 
\code{normalize.AffyBatch.vsn}. The resulting plot is shown in 
Fig.~\ref{vsn-affy}.
<<affy,include=FALSE,results=hide,fig=TRUE,width=7,height=7>>=
library(affy)
library(affydata)
data(Dilution)

## let affy know about vsn
normalize.AffyBatch.methods <- c(normalize.AffyBatch.methods, "vsn")

es1 = expresso(Dilution[1:2], 
               bg.correct       = FALSE,   ## bg correction is done by vsn
               normalize.method = "vsn",
               pmcorrect.method = "pmonly", 
               summary.method   = "medianpolish")

es2 = expresso(Dilution[1:2], 
               bgcorrect.method = "rma",
               normalize.method = "quantiles", 
               pmcorrect.method = "pmonly",
               summary.method   = "medianpolish")

## extract expression values
x1 = exprs(es1)
x2 = exprs(es2) 
 
## graphics output
oldpar = par(mfrow=c(2,2), pch=".")

## scatter plot
plot(x1, main="vsn: chip 3 vs 4")
plot(x2, main="rma: chip 3 vs 4")

## rank(mean) - difference plot
ylim = c(-0.7, 0.7)
plot(rank(rowSums(x1)), diff(t(x1)), ylim=ylim, main="rank(mean) vs differences")
abline(h=0, col="red")

plot(rank(rowSums(x2)), diff(t(x2)), ylim=ylim, main="rank(mean) vs differences")
abline(h=0, col="red")

## reset old plotting parameters 
par(oldpar)
@
\myincfig{vsn-affy}{\textwidth}{\code{normalize.AffyBatch.vsn} example}

%------------------------------------------------------------------------
\section{Verifying and assessing the performance of vsn with simulated
data} \label{sagmb}
%------------------------------------------------------------------------
There are two functions \Rfunc{sagmbSimulateData} and
\Rfunc{sagmbAssess} that can be used to generate simulated data and
assess the difference between the 'true' and 'estimated' data
calibration and transformation by vsn. An example is shown in the code
chunk below. Reference~\cite{HuberSAGMB2003} describes in more detail
(i) the simulation model, (ii) the assessment strategy, and (iii) a
comprehensive suite of assessments with respect to the number of
probes $n$, the number of arrays $d$, the fraction of differentially
expressed genes $de$, and the fraction of up-regulated genes $up$.
<<sagmb-calc,results=hide>>=
n     <- c(500, 1000, 2000, 4000, 8000)
d     <- 2
de    <- c(0, 0.2)
up    <- 0.5
nrrep <- 8
res <- array(NA, dim=c(length(n), nrrep, length(de)))
for (i in seq(along=de)) {
  for (k in seq(along=n)) {
    for (r in 1:nrrep) {
      sim          <- sagmbSimulateData(n[k], d, de=de[i], up=0.5)
      ny           <- vsn(sim$y)  
      res[k, r, i] <- sagmbAssess(exprs(ny), sim)
    }
  }
}
@
<<sagmb-show,include=FALSE,fig=TRUE,width=7,height=4>>=
par(mfrow=c(1,2))
for (i in seq(along=de)) {
  matplot(n, res[,,i], pch=20, log="xy", col="#909090", main=paste("de=", de[i]))
  lines(n, rowMeans(res[,,i]), col="blue")
}
@
\myincfig{vsn-sagmb-show}{\textwidth}{Estimation error for the transformation:
the root mean squared difference between true and estimated transformed data,
as a function of the number of genes $n$. If \Rfunc{vsn} works correctly, the
estimation error should decrease roughly as $n^{-1/2}$.
}

%---------------------------------------------------------
\section{Quality control}\label{sec.qc}
%---------------------------------------------------------
\code{vsn} makes some assumptions about your data that need to hold 
if it is to produce meaningful results. We have found them appropriate for 
many microarray experiments, but it is your responsibility to make sure 
that they hold for your data.  

First, \code{vsn} assumes that the measured signal $y_{ik}$ 
increases, to sufficient approximation, proportionally to the
mRNA abundance $c_{ik}$ of gene $k$ on the $i$-th array, or on
the $i$-th color channel:
\begin{equation}\label{assumption} 
y_{ik}\approx a_i + b_i b_k c_{ik}.  
\end{equation} 
For a series of $d$ single-color arrays such as Affymetrix arrays
or cDNA nylon membranes, $i=1,\ldots,d$, and the different factors 
$b_i$ reflect the different initial amounts of sample mRNA, or different
overall reverse transcription, hybridization and detection efficiencies. 
The probe affinity $b_k$ contains factors that affect
all measurements with probe $k$ in the same manner, such as
sequence-specific labelling 
efficiency. The $b_k$ are assumed to be the same across all arrays.
There can be a non-zero overall offset $a_i$ for each color channel.
For a two-color cDNA array, $i=1,2$, and the $b_i$ take into account
the different overall efficiencies of the two dyes\footnote{% 
It has been reported that for some genes the dye bias is different from 
gene to gene, such that the proportionality factor does not simply factorize 
as in~(\ref{assumption}). As long as this only occurs sporadically, this 
should not have much effect on the estimation of the calibration and 
variance stabilization parameters. Further, by using an appropriate 
experimental design such as color-swap or reference design, the effects of 
gene-specific dye-biases to subsequent analyses can also be reduced.}.

Situations in which the assumption~(\ref{assumption}) is violated include:

{\em Saturation.} The biochemical reactions and/or the photodetection 
can be run in such a manner that saturation effects occur.
It may be possible to rescue such data by using non-linear transformations. 
Alternatively, it is recommended that the experimental parameters are 
chosen to avoid saturation.

{\em Print-tip and PCR effects.} In cDNA microarray data, systematic 
patterns have been observed that are associated with the print-tips 
used for the spotting of DNA, and with the microtitre plates 
used to store and amplify the DNA. Possibly, these effects could be included
in~(\ref{assumption}) by replacing $a_i$ and $b_i$ by print-tip and/or plate 
specific coefficients, but presently \code{vsn} does not include this option.

{\em Batch effects.} The probe affinities $b_k$ may differ between different 
manufacturing batches of arrays due, e.g., to different qualities of DNA 
amplification or printing. \code{vsn} cannot be used to simultaneously 
calibrate and transform data from different batches.

How to reliably diagnose and deal with such violations is beyond the scope 
of this vignette; see the references for more~\cite{Dudoit578,HSG2002}. 

\paragraph{Variance.} A further assumption that \code{vsn} makes is that 
the measurement error (more exactly: the variance) is the sum of two 
contributions: an additive component that has roughly the same size for all 
probes on an array, and a multiplicative component that is roughly proportional 
in size to the signal's true value, with a proportionality factor (called
the {\em coefficient of variation}) that is the same for all 
genes~\cite{RockeDurbin2001}.

\paragraph{Most genes unchanged assumption.} 
\code{vsn} assumes that only a minority of genes on the arrays is
detectably differentially transcribed across the experiments. The
allowed size of that minority is controlled by the parameter
\code{lts.quantile} and must be less than or equal to $1/2$.  
By default, \code{lts.quantile}=0.5, that is, the maximum.

\paragraph{Processing biases.} 
Image analysis software for cDNA arrays typically estimates a {\em
local background} associated with each probe intensity. For Affymetrix
arrays, the intensities from {mismatch} probes are thought to
represent the level of non-specific signal. In both cases, the raw
probe intensities may be {\em adjusted} by subtracting these
background estimates. Some software packages, however, bias the
adjustment through rules based on the data values. For example,
Affymetrix' MAS 5.0 software uses the mismatch intensity only if it is
smaller than the probe's intensity, and otherwise employs a heuristic to
make sure that the net intensities always remain positive. As a
consequence, the intensities are systematically over-estimated, and
cannot be used with \code{vsn}. For Affymetrix data, we recommend to
use \code{vsn} on the probe intensities from the "CEL file". For cDNA
data, we recommend to use only background adjustment procedures that
estimate the background independent of the observed foreground
intensity.

\begin{thebibliography}{10}

\bibitem{HuberISMB2002}
W. Huber, A. von Heydebreck, H. {S\"ultmann}, A. Poustka, and M. Vingron.
\newblock Variance stablization applied to microarray data calibration and to
  quantification of differential expression.
\newblock {\em Bioinformatics}, 18:S96--S104, 2002.

\bibitem{HuberSAGMB2003}
W. Huber, A. von Heydebreck, H. {S\"ultmann}, A. Poustka, and M. Vingron.
\newblock Parameter estimation for the calibration and variance stabilization 
of microarray data.
\newblock {\em Statistical Applications in Genetics and Molecular Biology}, 
Vol. 2: No. 1, Article 3, 2003. 
http://www.bepress.com/sagmb/vol2/iss1/art3

\bibitem{HSG2002}
W. Huber, A. von Heydebreck, and M. Vingron.
\newblock Analysis of microarray gene expression data.
\newblock To appear in the {\em Handbook of Statistical Genetics}, 2003.
Eds.: D. J. Balding, M. Bishop, C. Cannings. 
John Wiley \& Sons, Inc., 

\bibitem{RockeDurbin2001}
David~M. Rocke and Blythe Durbin.
\newblock A model for measurement error for gene expression analysis.
\newblock {\em Journal of Computational Biology}, 8:557--569, 2001.

\bibitem{Dudoit578}
S. Dudoit, Y.~H. Yang, T.~P. Speed, and M.~J. Callow.
\newblock Statistical methods for identifying differentially expressed genes in
  replicated {cDNA} microarray experiments.
\newblock {\em Statistica Sinica}, 12:111--139, 2002.

\bibitem{StoreyTibshirani}
J.~D. Storey and R.~Tibshirani.
\newblock SAM thresholding and false discovery rates for detecting differential
  gene expression in {DNA} microarrays.
\newblock To appear in: G. Parmigiani, E.~S. Garrett, R.~A. Irizarry and 
  S.~L. Zeger (eds.), The analysis of gene expression data: Methods and 
  software. New York, Springer 2003.

\bibitem{Alizadeh}
A.~A. Alizadeh et al.
\newblock Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling.
\newblock {\em Nature}, 403:503--511, 2000.

\end{thebibliography}
\end{document}





