% -*- mode:R -*-
%\VignetteIndexEntry{Estimation of local false discovery rate}
%\VignetteDepends{splines,stats,golubEsets,vsn}
%\VignetteKeywords{Gene expression analysis}
%\VignettePackage{twilight}

\documentclass[11pt,a4paper,fleqn]{report}

\usepackage{compdiag}
\usepackage{amsmath}
\usepackage[bf]{caption}
\setlength{\captionmargin}{30pt}
\fboxsep=2mm
%\parindent0mm
%\parskip1ex
\SweaveOpts{eps=false}

\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}
\newcommand{\Rmethod}[1]{{\textit{#1}}}
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}

  

\title{Estimation of local false discovery rate \bigskip \\ User's guide to the Bioconductor package \Rpackage{twilight} 1.1.0}%'
\author{Stefanie Scheid\footnote{Corresponding author: 
        \texttt{stefanie.scheid@molgen.mpg.de}}
        and Rainer Spang \bigskip \\
        \small Max Planck Institute for Molecular Genetics\\ 
         Computational Diagnostics\\ Department for Computation Molecular Biology \\
         Ihnestr. 63-73, D-14195 Berlin, Germany}
\reportnr{01}
\year{2004}
\abstract{This is the vignette of the Bioconductor compliant package \Rpackage{twilight}. We describe our implementation of a stochastic search algorithm to estimate the local false discovery rate. In addition, the package provides functions to test for differential gene expression in the common two-condition setting.}
\date{}
  


\begin{document}
\maketitle

<<eval=FALSE,echo=FALSE,results=hide>>=
q(save="no")

<<echo=FALSE,results=hide>>=
library(splines)
library(stats)
library(golubEsets)
library(vsn)
oldopt <- options(digits=3)
on.exit( {options(oldopt)} )
options(width=70)

@
\chapter{Introduction}

In a typical microarray setting with gene expression data observed under two conditions, the local false discovery rate describes the probability that a gene is not differentially expressed between the two conditions given its corrresponding observed score or $p$-value level. The resulting curve of $p$-values versus local false discovery rate offers an insight into the \textbf{twilight zone} between clear differential and clear non-differential gene expression. The Bioconductor compliant package \Rpackage{twilight} contains two main functions: Function \Rfunction{twilight.pval} performs a two-condition test on differences in means for a given input matrix or expression set (\Rclass{exprSet}) and computes permutation based $p$-values. Function \Rfunction{twilight} performs the successive exclusion procedure described in Scheid and Spang (2004) \cite{scheid04} to estimate local false discovery rates and effect size distributions.




@
\chapter{Implemented methods}
\section{Function \texttt{twilight.pval}: Testing effect sizes}

\fbox{
\begin{minipage}{0.95\textwidth}
\Rfunction{twilight.pval(xin, yin, method="fc", paired=FALSE, B=1000, yperm=NULL, balance=FALSE, quant.ci=0.95, s0=NULL, verbose=TRUE)}
\end{minipage}
}
\bigskip       
       
The input object \Rfunarg{xin} is either a pre-processed gene expression set of class \Rclass{exprSet} or any data matrix where rows correspond to genes and columns to samples. Each sample was taken under one of two distinct conditions, for example under treatment A or treatment B. The functions in package \Rpackage{twilight} are not limited to microarray data only but can be applied to any two-sample data matrix. However, it is necessary for both expression set or numerical matrix that values are on \textbf{additive scale} like log or arsinh scale. The function does not check or transform the data to additive scale. The input vector \Rfunarg{yin} contains condition labels of the samples. Vector \Rfunarg{yin} has to be numeric and dichotomous. Note that in terms of \textit{under}- and \textit{over}-expression, the samples of the higher labeled condition are compared to the samples of the lower labeled condition.

We are given a pre-processed matrix for samples belonging to two distinct conditions A and B, and gene expression values on additive scale. For gene $i$ in the experiment ($i=1,\dots,N$), $\bar \alpha_i$ is the mean expression under condition A and $\bar \beta_i$ is the mean expression under condition B. To test the null hypothesis of no differential gene expression, function \Rfunction{twilight.pval} compares the mean expression levels $\bar \alpha_i$ and $\bar \beta_i$. The current version offers three test variants: The classical $t$-test uses score $T_i$ with
\begin{equation}
  T_i=\frac{\bar \alpha_i - \bar \beta_i }{s_i},
\end{equation}
where $s_i$ denotes the root pooled variance. The $t$-test is called with \Rfunarg{method="t"}.

The $t$-test score can be misleadingly high if $s_i$ is very small. To overcome this problem, the $Z$-test enlarges the denominator by a fudge factor $s_0$ \cite{tusher01}, \cite{efron01}:
\begin{equation}
  Z_i=\frac{\bar \alpha_i - \bar \beta_i}{s_i + s_0}.
\end{equation}
The $Z$-test is called with \Rfunarg{method="z"}. Fudge factor $s_0$ is set to \Rfunarg{s0=NULL} by default and is only evaluated if \Rfunarg{method="z"}. In that case, it is the median of all root pooled variances $s_1,\dots,s_N$. However, the fudge factor can be chosen manually. Note that if \Rfunarg{method="z"} is chosen with \Rfunarg{s0=0}, the test call is altered to \Rfunarg{method="t"}, the $t$-test as described above.

The third variant is based on log ratios only with score
\begin{equation}
  F_i=\bar \alpha_i - \bar \beta_i.
\end{equation}
The distribution of scores $F_i$ under the alternative is called \textit{effect size distribution}. With expression values on log or arsinh scale, $\exp (|F_i|)$ is an estimator for the fold change. We call $\exp (|F_i|)$ the \textit{fold change equivalent score} \cite{scheid04}. Note that the package contains a function for plotting the effect size distribution which is only available if function \Rfunction{twilight.pval} was run with \Rfunarg{method="fc"}, the fold change test.

Function \Rfunction{twilight.pval} handles paired and unpaired data. In the unpaired case (\Rfunarg{paired=FALSE}), only one microarray was hybridized for each patient, like in a treatment and control group setting. In the paired case (\Rfunarg{paired=TRUE}), we observed expression values of the same patient under both conditions. The typical example are before and after treatment experiments, where each patient's expression was measured twice. The input arguments \Rfunarg{xin} and \Rfunarg{yin} do not need to be ordered in a specific manner. It is only necessary that samples within each group have the same order, such that the first samples of the two groups represent the first pair and so on.%'
As an example, we apply function \Rfunction{twilight.pval} on the leukemia data set of Golub et al.~(1999) \cite{golub99} as given in \Rfunction{library(golubEsets)}. For normalization, apply the variance-stabilizing method \Rfunction{vsn} in \Rfunction{library(vsn)} \cite{huber02}.
<<results=hide>>=
data(golubMerge)
golubNorm <- vsn(exprs(golubMerge))
id <- as.numeric(golubMerge$ALL.AML)

@
There are 72 samples either expressing acute lymphoblastic leukemia (ALL) or acute myeloid leukemia (AML). As the AML patients are labeled with ``2'' and ALL with ``1'', we compare AML to ALL expression.
<<>>=
golubMerge$ALL.AML
id

@
Additionally to computation of scores, empirical $p$-values are calculated. Argument \Rfunarg{B} specifies the number of permutations with default set to \Rfunarg{B=1000}. The distribution of scores under the null hypothesis is estimated by computing test scores from the same input matrix with randomly permuted class labels. These permutations are either balanced or unbalanced, with default \Rfunarg{balance=FALSE}. The permutation options are described in detail in section \ref{twilight.combi}. For computing empirical $p$-values, we count for each gene how many of \textit{all} absolute permutation scores exceed the absolute observed score, and divide by \Rfunarg{B}$\cdot$(number of genes).

Permutation scores are also used to compute expected scores as described in Tusher et al.~(2001) \cite{tusher01}. In addition, we compute confidence bounds for the maximum absolute difference of each set of permutation scores to expected scores. The width of the confidence bound is chosen with \Rfunarg{quant.ci}. With default \Rfunarg{quant.ci=0.95}, the maximum absolute difference of permutation to expected scores exceeded the confidence bound in only 5\% of all permutations.

Using the optional argument \Rfunarg{yperm}, a user-specified permutation matrix can be passed to the function. In that case, \Rfunarg{yperm} has to be a \textit{binary} matrix where each row is one vector of permuted class labels. The label "1" in \Rfunarg{yperm} corresponds to the higher labeled original class. If the permutation matrix is specified, no other permutation is done and argument \Rfunarg{B} will be ignored. Besides \Rfunction{set.seed}, argument \Rfunarg{yperm} can be used to reproduce results by fixing the matrix of random permutations.

Continuing the example above, we perform a fold change test on the expression data in \Robject{golubNorm} which was transformed to arsinh scale by normalization with \Rfunction{vsn}.
<<>>=
library(twilight)
pval <- twilight.pval(golubNorm, id)

@
The function checks whether complete enumeration of all permutations is possible. Complete enumeration is performed as long as the number of permutations does not exceed 10\,000. The admissible sample sizes depend on whether the experiment was paired or unpaired and whether balanced or unbalanced permutations are used. Details are given in section \ref{twilight.combi}.
The values in the accompanying data set \Robject{expval} were computed in the same manner as in the example above.
<<>>=
data(expval)
expval

@
The output object of function \Rfunction{twilight.pval} is of class \Rclass{twilight} with several elements stored in a list.
<<>>=
class(expval)
names(expval)

@
The element \Robject{quant.ci} contains the corresponding input value which is passed to the plotting function. Element \Robject{ci.line} is used for plotting confidence bounds and contains the computed quantile of maximum absolute differences.  The output dataframe \Robject{result} contains a matrix with several columns.
<<>>=
names(expval$result)

@
The dataframe stores observed and expected scores and corresponding empirical $p$-values. The genes are ordered by absolute observed test scores. Genes with observed score exceeding the confidence bounds are marked as ``1'' in the binary vector \Robject{result\$candidate}. The output object is passed to function \Rfunction{plot.twilight} to produce a plot as in Tusher et al.~(2001) \cite{tusher01} with additional confidence lines and genes marked as candidates, see Figure \ref{fig_scores}.
<<>>=
expval$result[1:7,1:5]

@
\begin{figure}[t]
\begin{center}
<<fig=TRUE,echo=FALSE,width=6,height=4.5>>=
plot(expval,which="scores",grayscale=F,legend=F)
@
\caption{Expected versus observed test scores. Deviation from the diagonal line gives evidence for differential expression. The red lines mark the 95\% confidence interval on the absolute difference between oberved and expected scores. The plotting call is \texttt{plot(expval,which="scores",grayscale=F,legend=F)}.}\label{fig_scores}
\end{center}
\end{figure}

@
In addition, $q$-values and the estimated percentage of non-induced genes $\pi_0$ are computed as described in Remark B of Storey and Tibshirani (2003) \cite{storey03}. These are stored in \Robject{result\$qvalue} (see above) and \Robject{pi0}. The remaing output elements of \Robject{expval} are left free to be filled by function \Rfunction{twilight}. With \Rfunarg{"qvalues"}, Figure \ref{fig_qvalues} shows the plot of $q$-values against the corresponding number of rejected hypotheses.
<<>>=
expval$pi0
@

@
\begin{figure}[t]       
\begin{center}
<<fig=TRUE,echo=FALSE,width=6,height=4.5>>=
plot(expval,which="qvalues")
@
\caption{Stairplot of $q$-values against the resulting size of list of significant genes. A list containing all genes with $q \leq q_0$ has an estimated global false discovery rate of $q_0$. The plotting call is \texttt{plot(expval,which="qvalues")}.}\label{fig_qvalues}
\end{center}
\end{figure}

Column \Robject{result\$index} contains the original gene ordering of the input object. With these numbers, resorting of the \Robject{result} table is possible without knowing the original order of the row names.





@
\clearpage
\section{Function \texttt{twilight.pval}: Testing correlation}

\fbox{
\begin{minipage}{0.95\textwidth}
\Rfunction{twilight.pval(xin, yin, method="fc", B=1000, quant.ci=0.95, verbose=TRUE)}
\end{minipage}
}
\bigskip       

From version 1.1.0 on, function \Rfunction{twilight.pval} offers the computation of correlation scores instead of effect size scores. Now, vector \Rfunarg{yin} can be any clinical parameter consisting of numerical values and having length equal to the number of samples. With \Rfunarg{method="pearson"}, Pearson's coefficient of correlation to \Rfunarg{yin} is computed for every gene in \Rfunarg{xin}. With \Rfunarg{method="spearman"}, \Rfunarg{yin} and the rows of \Rfunarg{xin} are converted into ranks and Spearman's rank correlation is computed.

Note that most input arguments of \Rfunction{twilight.pval} will be ignored. Only \Rfunarg{B} takes effect and causes the computation of correlation scores based on \Rfunarg{B} random permutations of \Rfunarg{yin}. All successive analyses like expected scores, $p$- and $q$-values are kept as before. As an illustration, we search for genes with high correlation to the highest scoring gene found in the effect size test. Figure \ref{fig_corr} displays the resulting scores.

<<>>=
gene <- exprs(golubNorm)[expval$result$index[1],]
corr <- twilight.pval(golubNorm,gene,method="spearman",quant.ci=0.99)
corr
@
Note that the overall percentage of non-induced genes $\pi_0$ is now interpreted as the overall percentage of genes not correlated to the clinical parameter under the null hypothesis.

<<>>=
corr$result[1:10,1:5]
@

@
\begin{figure}[t]
\begin{center}
<<fig=TRUE,echo=FALSE,width=6,height=4.5>>=
plot(corr,which="scores",grayscale=F,legend=F)
@
\caption{Expected versus observed Spearman correlation scores. Deviation from the diagonal line gives evidence for significant correlation. The red lines mark the 99\% confidence interval on the absolute difference between oberved and expected scores. The plotting call is \texttt{plot(corr,which="scores",grayscale=F,legend=F)}.}\label{fig_corr}
\end{center}
\end{figure}




@
\clearpage
\section{Function \Rfunction{twilight}}

\fbox{
\begin{minipage}{0.95\textwidth}
\Rfunction{twilight(xin, lambda=NULL, B=0, boot.ci=0.95, clus=NULL, verbose=TRUE)}
\end{minipage}
}
\bigskip       

Local false discovery rates (fdr) are estimated from a simple mixture model given the density $f(t)$ of observed scores $T=t$:
\begin{equation}
f(t) = \pi_0 \, f_0(t) + (1-\pi_0) \, f_1(t) \quad \Rightarrow \quad \mbox{fdr}(t)=\pi_0 \, \frac{f_0(t)}{f(t)},
\end{equation}
where $\pi_0 \in [0,1]$ is the overall percentage of non-induced genes. Terms $f_0$ and $f_1$ are score densities under no induction and under induction respectively. Assume that there exists a transformation $W$ such that $U=W(T)$ is uniformly distributed in $[0,1]$ for all genes not differentially expressed. In a multiple testing scenario these $u$-values are $p$-values corresponding to the set of observed scores. However, we do not regard the local false discovery rate as a multiple error rate but as an exploratory tool to describe a microarray experiment over the whole range of significance, thus speaking of $u$-values instead of $p$-values.

Mapping scores to $u$-values allows to assume $f_0(u)$ to be the uniform density instead of specifying the null density $f_0(t)$ with respect to a chosen scoring method. The implemented successive exclusion procedure (SEP) splits any vector of $u$-values into a uniformly distributed null part and an alternative part. The uniform part represents genes that are not differentially expressed. Based on histogram counts of the null and the alternative part, we estimate local false discovery rates. These denote the conditional probability of being non-induced given we observe a certain $u$-value level. The proportion of the uniform part to the total number of genes in the experiment is a natural estimator for percentage $\pi_0$. The successive exclusion procedure is described in detail in Scheid and Spang (2004) \cite{scheid04}. The functionality of \Rfunction{twilight} is not limited to microarray experiments. In principle, any vector of $u$-values can be passed to \Rfunction{twilight} as long as the assumption of uniformity under the null hypothesis is valid.


The objective function in \Rfunction{twilight} includes a penalty term that is controlled by regularization parameter $\lambda \geq 0$. The regularization ensures that we find a separation such that the uniform part contains as many $u$-values as possible. As percentage $\pi_0$ is often underestimated, the inclusion of a penalty term results in a more ``conservative'' estimate that is usually less biased. If not specified (\Rfunarg{lambda=NULL}), function \Rfunction{twilight.getlambda} finds a suitable $\lambda$.

The estimates for probability $\pi_0$ and local false discovery rate are averaged over 10 runs of SEP. In addition, bootstrapping can be performed to give bootstrap estimates and bootstrap percentile confidence intervals on both $\pi_0$ and local false discovery rate. The number of bootstrap samples is set by argument \Rfunarg{B}, and the width of the bootstrap confidence interval is set by argument \Rfunarg{boot.ci}.

Function \Rfunction{twilight} takes \Rclass{twilight} objects or any vector of $u$-values as input and returns a \Rclass{twilight} object. If the input is of class \Rclass{twilight}, the function works on the set of empirical $p$-values and fills in the remaining output elements. Note that the estimate for $\pi_0$ is replaced, and $q$-values are recalculated with the new estimate $\pi_0$.
As an example, we run SEP with 1000 bootstrap samples and $95\%$ boostrap confidence intervals: \Rfunction{twilight(xin=expval, B=1000, boot.ci=0.95)}, as was done for data set \Robject{exfdr}.
<<>>=
data(exfdr)
exfdr

<<>>=
exfdr$result[1:5,6:9]

@
The output elements \Robject{result\$fdr}, \Robject{result\$mean.fdr}, \Robject{result\$lower.fdr} and \linebreak[5] \Robject{result\$upper.fdr} contain the estimated local false discovery rate, the bootstrap average and upper and lower bootstrap confidence bounds. These values are used to produce the following plots which are only available after application of function \Rfunction{twilight}. First, we plot $u$-values against the corresponding conditional probabilities of being induced given the $u$-value level, that is $1-\mbox{fdr}$, see Figure \ref{fig_fdr}. Going back to observed scores, we produce a \textit{volcano plot}, that is observed scores versus local false discovery rate, see Figure \ref{fig_volcano}.

@
\begin{figure}[t]
\begin{center}
<<fig=TRUE,echo=FALSE,width=6,height=4.5>>=
plot(exfdr,which="fdr",grayscale=F,legend=T)
@
\caption{Curve of estimated local false discovery over $u$-values. The red lines denote the bootstrap mean (solid line) and the 95\% bootstrap confidence interval on the local false discovery rate (dashed lines). The bottom ticks are 1\% quantiles of $u$-values. The plotting call is  \texttt{plot(exfdr,which="fdr",grayscale=F,legend=T)}.}\label{fig_fdr}
\end{center}
\end{figure}


@
\begin{figure}[t]
\begin{center}
<<fig=TRUE,echo=FALSE,width=6,height=4.5>>=
plot(exfdr,which="volcano")
@
\caption{Volcano plot of observed test scores versus local false discovery rate. The bottom ticks are 1\% quantiles of observed scores. The plotting call is \texttt{plot(exfdr,which="volcano")}.}\label{fig_volcano}
\end{center}
\end{figure}

@
Output element \Robject{effect} contains histogram information about the effect size distribution, that is log ratio under the alternative. One run of the successive exclusion procedure results in a split of the input $u$-value vector into a null and an alternative part. We estimate the effect size distribution from the distribution of log ratio scores corresponding to $u$-values in the alternative part. Again, this estimate is averaged over 10 runs of the procedure. Argument \Rfunarg{which="effectsize"} produces the histogram of all observed log ratios overlaid with the averaged histogram of log ratios in the alternative, see Figure \ref{fig_effectsize}. The x-axis is changed to fold change equivalent scores or rather to increase in effect size. Given an observed log ratio $F$, the increase in effect size is $(\exp (|F|)-1) \cdot \mbox{sign} (F) \cdot 100\%$. A value of 0\% corresponds to no change (fold change of 1), a value of 50\% to fold change 1.5 and so on. A value of -100\% corresponds to a 2-fold down-regulation, that is fold change of 0.5.

@
\begin{figure}[t]
\begin{center}
<<fig=TRUE,echo=FALSE,width=7,height=4.5>>=
plot(exfdr,which="effectsize",legend=T)
@
\caption{Observed effect size distribution (gray histogram) overlaid with the estimated effect size distribution under the null hypothesis (black histogram). The plotting call is \texttt{plot(exfdr,which="effectsize",legend=T)}.}\label{fig_effectsize}
\end{center}
\end{figure}

@
The last plotting argument \Rfunarg{which="table"} tabulates the histogram information in terms of fold change equivalent scores and log ratios.
<<>>=
tab <- plot(exfdr,which="table")
tab[1:8,]
@

       
@
The input argument \Rfunarg{clus} of function \Rfunction{twilight} is used to perform parallel computation within \Rfunction{twilight}. Parallelizing saves computation time which is especially useful if the number of bootstrap samples \Rfunarg{B} is large. With default \Rfunarg{clus=NULL}, no parallelizing is done. If specified, \Rfunarg{clus} is passed as input argument to \Rfunction{makeCluster} in \Rfunction{library(snow)}. Please make sure that \Rfunction{makeCluster(clus)} works properly in your environment.






@
\clearpage
\section{Function \Rfunction{twilight.combi}}\label{twilight.combi}

\fbox{
\begin{minipage}{0.95\textwidth}
\Rfunction{twilight.combi(xin, pin, bin)}
\end{minipage}
}
\bigskip       

Function \Rfunction{twilight.combi} is used within \Rfunction{twilight.pval} to completely enumerate all permutations of a \textit{binary} input vector \Rfunarg{xin}. Argument \Rfunarg{pin} specifies whether the input vector corresponds to paired or unpaired data. Argument \Rfunarg{bin} specifies whether permutations are balanced or unbalanced. Note that the resulting permutations are always ``as balanced as possible'': The balancing is done for the smaller subsample. If its sample size is odd, say 7, \Rfunction{twilight.combi} computes all permutations with 3 and 4 samples unchanged.

As first example, compute all unbalanced permutations of an unpaired binary vector of length 5 with two zeros and three ones. The number of rows are 
\begin{equation}
m = \frac{5!}{2! \cdot 3!} = 10.
\end{equation}
<<>>=
x <- c(rep(0,2),rep(1,3))
x

<<>>=
twilight.combi(x,pin=FALSE,bin=FALSE)

@
Each row contains one permutation. The first row contains the input vector. In balanced permutations, we omit those rows where both original zeros have been shifted to the last three columns. The number of balanced rows is
\begin{equation}
m = {2 \choose 1} \cdot \frac{3!}{1! \cdot 2!} = 6.
\end{equation}
<<>>=
twilight.combi(x,pin=FALSE,bin=TRUE)
@
Note that the function returns six balanced rows \textit{and} the original input vector although it is not balanced.
       
Next, consider a paired input vector with four pairs. The first zero and the first one are the first pair and so on. In paired settings, values are flipped only within a pair. The number of rows is
\begin{equation}
m = \frac{1}{2} \cdot 2^4 = 2^3 = 8.
\end{equation}
<<>>=
y <- c(rep(0,4),rep(1,4))
y

<<>>=
twilight.combi(y,pin=TRUE,bin=FALSE)
@
       
The matrix above contains only half of all possible $2^4=16$ permutations. The reversed case \Rfunction{1 - twilight.combi(y, pin=TRUE, bin=FALSE)} is omitted as this will lead to the same absolute test scores as \Rfunction{twilight.combi(y, pin=TRUE, bin=FALSE)}. The same concept applies to balanced paired permutations. Now, two pairs are kept fixed and two pairs are flipped in each row. The number of balanced rows is
\begin{equation}
m = \frac{1}{2} \cdot {4 \choose 2} = 3.
\end{equation}
<<>>=
twilight.combi(y,pin=TRUE,bin=TRUE)
@
Again, the input vector is part of the output.

The complete enumeration of \Rfunction{twilight.combi} is limited by the sample sizes. The function returns \Rfunarg{NULL} if the resulting number of rows exceeds 10\,000. If \Rfunarg{NULL} is returned, function \Rfunction{twilight.pval} uses the functions \Rfunction{twilight.permute.unpair} and \Rfunction{twilight.permute.pair} which return a matrix of random permutations. For example, use the latter function to compute 7 balanced permutations of the paired vector \Rfunarg{y}. Similar to \Rfunction{twilight.combi}, these two functions return the input vector in the first row of their output matrices.
<<>>=
twilight.permute.pair(y,7,bal=TRUE)
@






@
\section*{Acknowledgements}

This work was done within the context of the Berlin Center for Genome Based Bioinformatics (BCB), part of the German National Genome Network (NGFN), and supported by BMBF grants 031U109C and 03U117 of the German Federal Ministry of Education.








\chapter{Differences to earlier versions}


\section{Changes in version 1.1.0}

The computation of $p$-values in \Rfunction{twilight.pval} changed from gene-wise to pooled $p$-values. For the computation of a gene-wise $p$-value for gene $i$, only the permutation scores of gene $i$ are taken into account. For pooled $p$-values, all permutation scores of all genes are taken as null distribution. This change has several advantages: First, gene-wise $p$-values were not monotonically increasing with scores because each gene had its own null distribution. Thus, two genes with almost equal scores might get quite different $p$-values. Now, the null distribution is the same for all genes, that is the union set of all permutation scores. Second, pooled $p$-values are less granular than gene-wise $p$-values. Gene-wise $p$-values are computed from \Rfunarg{B} permutation scores whereas pooled $p$-values are computed from \Rfunarg{B} $\cdot$ (number of genes) scores.

These two important features gave rise to further changes: The ordering of the \Robject{result} table is now more intuitive because the most significant genes on top have the highest scores, the lowest $p$- and $q$-values and are candidates (if there are any). In addition, the default value of the number of permutations \Rfunarg{B} is lowered to 1000 permutations. Computation of pooled $p$-values is slower than for gene-wise $p$-values. On the other hand, changing to pooled $p$-values increases the number of values in the null distribution by the factor of number of genes. Hence, even with less permutations, the number of null values is larger than before.

We integrated Pearson and Spearman correlation coefficients into \Rfunction{twilight.pval}. Each gene is correlated to an numerical input vector. Expected scores are computed from random permutations of the input vector.
       
The \Robject{result} table contains an additional \Robject{index} column with genes indices which comes in handy for sorting back to original ordering.

All output matrices of the permutation functions \Rfunction{twilight.combi}, \Rfunction{twilight.permute.pair} and \Rfunction{twilight.permute.unpair} have the original labeling vector as first row. This is also the case if balanced permutations are wanted, although the input vector is not balanced. Hence, the permutation matrix within \Rfunction{twilight.pval} now includes the original labeling even for balanced permutations implying that the smallest \textit{possible} $p$-value is 1/(number of permutations).


\section{Changes in version 1.0.3}

A \Rfunction{print.twilight} function was added which produces a short information about the contents stored in the \Robject{twilight} object.


\section{Changes in version 1.0.2}

The \Rfunarg{which} argument of the plot command changed from \Rfunction{plot1} style to more intuitive labels like \Rfunarg{scores} or \Rfunarg{fdr}.

       




       







       
\begin{thebibliography}{1}

\bibitem{efron01}
B. Efron, R. Tibshirani, J.D. Storey and V.G. Tusher, "Empirical Bayes Analysis of a Microarray Experiment", \emph{J. Am. Stat. Assoc.}, vol. 96, no. 456, pp. 1151-1160, 2001.

\bibitem{huber02}
W. Huber, A. von Heydebreck, H. S{\"u}ltmann, A. Poustka and M. Vingron, ``Variance stabilization applied to microarray data calibration and to the quantification of differential expression'', \emph{Bioinformatics}, vol. 18, suppl. 1, pp. S96-S104, 2002.

\bibitem{golub99}
T.R. Golub, D.K. Slonim, P. Tamayo, C. Huard, M. Gaasenbeek, J.P. Mesirov, H. Coller, M.L. Loh, J.R. Downing, M.A. Caligiuri, C.D. Bloomfield and E.S. Lander, "Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring", \emph{Science}, vol. 286, pp. 531-537, 1999.

\bibitem{scheid04}
S. Scheid and R. Spang, "A stochastic downhill search algorithm for estimating the local false discovery rate", \emph{IEEE Transactions on Computational Biology and Bioinformatics}, vol. 1, no. 3, pp. 98-108, 2004.

\bibitem{storey03}
J.D. Storey and R. Tibshirani, ``Statistical significance for genomewide studies'', \textit{Proc. Natl. Acad. Sci.}, vol. 100, no. 16, pp. 9440-9445, 2003.

\bibitem{tusher01}
V. Tusher, R. Tibshirani and C. Chu, ``Significance analysis of microarrays applied to ionizing radiation response'', \textit{Proc. Natl. Acad. Sci.}, vol. 98, no. 9, pp. 5116-5121, 2001.

\end{thebibliography}




\end{document}
