% \VignetteIndexEntry{GlobalAncova}
% \VignetteDepends{GlobalAncova, multtest, globaltest, golubEsets, hu6800, vsn}
% \VignetteKeywords{Expression Analysis}
% \VignettePackage{GlobalAncova}

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

\documentclass[a4paper]{article}

\title{Global testing of differential gene expression}

\author{Manuela Hummel \and Reinhard Meister \and Ulrich Mansmann} %\and Ramona Scheufele }


\begin{document}

\maketitle \tableofcontents \newpage

\section{Abstract}

In studies about differential gene expression between different clinical diagnoses
the main interest may often not be
in single genes but rather in groups of genes that are associated with a pathway
or have a common location in the genome.
In such cases it may be better to perform a global test because the problems of multiple testing
can be avoided. 
The approach presented here is an ANCOVA global test on phenotype effects
and gene--phenotype interaction.

Testing many pathways simultaneously is also possible. This, of course, causes again need for correction
for multiple testing. Besides the standard approaches for correction we introduce a closed
testing procedure in which the experiment--wise error rate equals the required level of 
confidence of the overall test.  

<<echo=F, results=hide>>=
#library(Biobase)
library(GlobalAncova)
library(globaltest)
sI <- sessionInfo()
@

This document was created using R version \Sexpr{paste(R.version$major,".",R.version$minor,sep="")}
and versions \Sexpr{sI$otherPkgs$GlobalAncova$Version} 
and \Sexpr{sI$otherPkgs$globaltest$Version} of the packages \Rpackage{GlobalAncova} and \Rpackage{globaltest}, respectively.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Changes to Previous Versions} 

\noindent {\bf Version 3.3.3}

\begin{itemize}
\item The permutation approach is now implemented in \emph{C} and therefore faster.
\item If the number of possible phenotype permutations is smaller than the number specified in \Rfunarg{perm} (i.e. in very small sample sizes), 
all possible permutations are considered for the permutation test.   
\item Some more error messages are included.
\item In \Rfunction{Plot.genes} and \Rfunction{Plot.subjects} bar labels can be manipulated with the argument \Rfunarg{bar.names}.
\end{itemize}

\vspace{0.2cm}
\noindent {\bf Version 3.x.x}

\begin{itemize}
\item Besides the permutation-based p-values also asymptotic p-values based on an approximation of the distribution of the test statistic are provided. Theoretical F-test p-values are no longer displayed since they are not valid in case of correlations or non-normality. 
\item The focus level procedure for finding interesting Gene Ontology subgraphs from Goeman and Mansmann (2007) \cite{GoeMan:07} was adapted for the use with \Rfunction{GlobalAncova}.
\item Sequential and type III decompositions of the residual sum of squares, adjustment for global covariates and pair-wise comparisons of different levels of a categorical factor are implemented. These functionalities are described in the additional vignette \emph{GlobalAncovaDecomp.pdf}.
\item Now the parameter \Rfunarg{test.genes} allows for specifying the gene group which a graph shall be based on in the plotting functions. 
\end{itemize}

\vspace{0.2cm}
\noindent {\bf Version 2.5.1}

\begin{itemize}
\item Testing several groups of genes is now more efficient and less time consuming.
\item In the gene and subjects plots bar heights (gene-wise reduction in sum of squares and subject-wise reduction in sum of squares, respectively) can be returned. %This eases the identification of genes with most influence regarding differential expression or of subjects behaving like outliers.
\item Plots are more flexible regarding graphical parameters like specification of colors, titles, axis labels and axis limits.
\end{itemize}

\vspace{0.2cm}
\noindent {\bf Version 2.x.x}

\begin{itemize}
\item The major modification in the new version is the transfer from simple two group comparisons to a general linear model framework where arbitrary clinical variables (in especially with more groups or also continuous ones), time trends, gene--gene interactions, co--expression and so forth can be analysed. 
\item According to the new framework also the diagnostic plots are more flexible. The variable defining the coloring of bars can now be specified by the user, see section \ref{plots} for details.
\item A bug was fixed concerning testing only a single gene for differential expression with the global ANCOVA F--test.
\item Within the closed testing procedure a bug was fixed concerning testing non--disjunct groups of genes.
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

%The ANCOVA global test is a simultaneous test on phenotype main effect and gene--phenotype
%interaction in a two--way layout linear model. 
The ANCOVA global test is a test for the association between expression values and clinical entities.
The test is carried out by comparison of linear models via the extra sum of squares principle.
If the mean expression level for at least one
gene differs between corresponding models the global null hypothesis, which is the intersection of all
single gene null hypotheses, is violated.
As our test is based on the sum of gene-wise reduction in sum of squares due to phenotype, all systematic differences in gene expression between phenotypes equally contribute to the power of the test.

Single genes are not, in general, the primary focus of gene expression experiments. The 
researcher might be more interested in relevant pathways, functional sets or genomic regions
consisting of several genes.  
Most of the current methods for studying pathways analyse differential
expression of single genes. In these methods pathways where many genes show minor changes in their expression
values may not be identified.
Goeman's global test and the ANCOVA global test were designed to address this issue.

Applying global tests for differential expression in pathways substantially reduces the number of tests compared to gene-wise multiple 
testing. The amount of correction for multiple testing decreases.
Function (KEGG, GO) or location (chromosome, cytoband) could be used as grouping criteria, for example.

We want to compare our method with the global test of Goeman et al., 2004 \cite{Goe:04}.
Therefore text and examples in this document follow to a certain extent the vignette presented in the R-package 
\Rpackage{globaltest}.
Our function \Rfunction{GlobalAncova} tests whether the expectation of expression levels differs between biological entities for a 
given group of genes.
This vignette has its focus on the practical use of the test. 
For more details about the mathematical background and the interpretation of results,
we refer to the paper by Mansmann and Meister, 2005 \cite{MaMei:05}
and the technical report of Hummel, Meister and Mansmann \cite{Hummel:07}.

This document shows the functionality of the R-package \Rpackage{GlobalAncova}. 
The datasets, all necessary R--packages and our package \Rpackage{GlobalAncova}
are available from the Bioconductor website ({\it www.bioconductor.org}).

First we load the packages and data we will use.

<<initialize, results=hide>>=
library(GlobalAncova)
library(globaltest)
library(golubEsets)
library(hu6800)
library(vsn)
library(multtest)

data(Golub_Merge)
golubX <- justvsn(Golub_Merge)
@

This creates a dataset \Robject{golubX}, which is of the format
\Rclass{ExpressionSet},
the standard format for gene expression data in
BioConductor. It consists of \Sexpr{length(featureNames(golubX))} genes and \Sexpr{length(sampleNames(golubX))}
samples (the data are from \cite{Gol:99}). We used \Rpackage{vsn}
to normalize the data. Other appropriate normalization
methods may be used as well. 
From several phenotype variables we use ``ALL.AML'' as the clinical diagnoses of interest.
ALL and AML are two types of acute leukemia. 
There are \Sexpr{sum(golubX$ALL.AML=="ALL")} patients
with ALL and \Sexpr{sum(golubX$ALL.AML=="AML")} with AML.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Global Testing of a Single Pathway} \label{onepathway}

\subsection{Golub Data and Cell Cycle Pathway}

Suppose we are interested in testing whether AML and ALL have
different gene expression patterns for certain pathways, for example from the
KEGG database. 
With \Rpackage{globaltest} we answer the question whether the expression profile has prognostic power with respect
to diagnosis of AML or ALL.
\Rpackage{GlobalAncova} asks for differences in mean expression between the two clinical groups.

\subsubsection{Testing all Genes}

We start by applying our test to all genes in the Golub dataset so that differences in the overall gene-expression pattern can be 
demonstrated. 

<<echo=F>>=
options(width=70)
@

<<all>>=
gr <- as.numeric(golubX$ALL.AML=="ALL")
ga.all <- GlobalAncova(xx=exprs(golubX), group=gr, covars=NULL, perm=100)
@ 

The first input \Rfunarg{xx} is a \Sexpr{dim(exprs(golubX))[[1]]} $\times$ 
\Sexpr{dim(exprs(golubX))[[2]]} matrix that contains the expression values of all genes and samples. 
The second input \Rfunarg{group} is a vector that
defines the clinical diagnosis for the \Sexpr{dim(exprs(golubX))[[2]]} patients. 
 
Note that \Rpackage{GlobalAncova} is not restricted to the analysis of dichotomous phenotype groups. More complex tasks like variables 
with more groups or also continuous ones, time trends, gene--gene interactions and co--expression can be performed as well. Some examples 
will be given in section \ref{vantVeer}. The realization of such tasks is done by definition of two linear models that shall be compared 
via the extra sum of squares principle. Hence model formulas for the full model containing all parameters and the reduced model, where 
the terms of interest are omitted, have to be given. An alternative is to provide the formula for the full model and a character vector 
naming the terms of interest. 
Those names can be chosen by previous output of the \Rfunction{GlobalAncova} function. 
Consequently we could run the same analysis as above with two possible further function calls shown below (output is omitted). In both cases a data frame with information about all variables for each sample is required. In the case of microarray data this can be the corresponding \Robject{pData} object. 
%However, one has to step gently with variable names including dots or colons. Those could be interpreted as contrasts or interaction terms and therefore cause problems in the function. %Thus here we build a new model data frame.

<<all2, eval=F>>= 
GlobalAncova(xx=exprs(golubX), formula.full=~ALL.AML, formula.red=~1, model.dat=pData(golubX), perm=100)
GlobalAncova(xx=exprs(golubX), formula.full=~ALL.AML, test.terms="ALL.AMLAML", model.dat=pData(golubX), perm=100)
@  

To avoid alpha--inflation due to correlated data and effects of non--normality of the data tests for significance of the resulting F--ratios are performed using a 
permutation test approach. 
We apply permutation of samples which is equivalent to permuting rows of the full design matrix. Note that permutation is only 
conducted for such columns of the design matrix that correspond to the variables of interest. Values of additional covariates 
remain in the original order. This prevents us for destroying covariate effects. Still the permutation approach is not optimal since 
residuals may be correlated. However, this does not seem to be a severe problem.  
The argument \Rfunarg{perm} defines the number of permutations, which is 10,000 for default. Here we set \Rfunarg{perm} to just 100 or 1000 so that creating this vignette 
will not last too long. For getting more reliable results one should recompute the examples with more permutations. \\
As an alternative to the permutation approach an approximation of the F--statistic nominator according to \cite{Rob:49} yields 
asymptotic p--values. 
%To use the approximation specify the parameter \Rfunarg{method = "approx"} (\Rfunarg{method = "permutation"} is the default) or \Rfunarg{method = "both"} to get both p--values. 
Note that the approximation is not feasible for very large gene groups since the huge gene expression covariance matrix has to be estimated, which is not possible for too many genes. 
The default value for group size (\Rfunarg{max.group.size}) is 2500, groups above this size are treated by the permutation approach. 
When using work stations with good working memory this number may be increased.
The estimation of the covariance matrix is carried out with the R package \Rpackage{corpcor} from \cite{corpcor}.\\
Whether the permutation--based or the asymptotic p--values or both should be calculated is controlled by the argument \Rfunarg{method}. 

The result of the \Rfunction{GlobalAncova} function is a typical ANOVA table with information about sums of squares, degrees of freedom 
and mean sums of squares for the effect and error term, respectively. 
Besides F--statistics %and classical F--test p--values,
there are given either p--values from the permutation test or the asymptotic p--values or both. 
%Which of them are calculated is controled by the argument \Rfunarg{method}.  
The names of all involved parameters are displayed as well as the name(s) of the tested effect(s).

<<all.display>>= 
ga.all 
@ 

From this result we conclude that 
the overall gene expression profile for all \Sexpr{length(featureNames(golubX))} genes is associated 
with the clinical outcome. 
This means that samples with different AML/ALL status tend to have different expression profiles.
We expect most pathways (especially the ones containing many genes) also to be associated with the phenotype groups.

If we apply Goeman's global test we get

<<gt.all, results=hide>>=
gt.all <- globaltest(golubX, "ALL.AML")
gt.all
@
<<echo=F>>=
gt.all
@

Both tests show that the data contain overwhelming evidence for differential gene expression
between AML and ALL.

\subsubsection{Testing the Cell Cycle Pathway}

Now we ask the more specific question of whether there is evidence for differential gene expression
between both diagnoses restricted to genes belonging to the cell cycle pathway.
First we load all KEGG pathways. 

<<loadKEGG>>=
kegg <- as.list(hu6800PATH2PROBE)
@

<<echo=F>>=
cellcycle <- kegg[["04110"]]
@

The list \Robject{kegg} consists of \Sexpr{length(kegg)} pathways. Each pathway is represented
by a
vector of gene names. We are mainly interested in the cell cycle pathway 
which has the identifier ``04110'' 
in the KEGG database. It corresponds to \Sexpr{length(cellcycle)} probe sets on the 
hu6800 chip. 

<<cc, eval=F>>=
cellcycle <- kegg[["04110"]]
@

We apply the global test to this pathway using the option \Rfunarg{test.genes}.

<<cellcycle>>=
ga.cc <- GlobalAncova(xx=exprs(golubX), group=gr, test.genes=cellcycle, method="both", perm=1000) 
ga.cc
@

Also with \Rpackage{globaltest} we get a very small p--value

<<gt.cellcycle, results=hide>>=
gt.cc <- globaltest(X=golubX, Y="ALL.AML", genesets=cellcycle)
gt.cc
@
<<echo=F>>=
gt.cc
@

The test results clearly indicate that the expression pattern of the cell cycle pathway 
is different between the two clinical groups. 

\subsubsection{Adjusting for Covariates}

Covariate information can be incorporated by specifying the \Rfunarg{covars} option.

For example if we want to adjust for whether samples were taken from bone marrow or from peripheral blood (\Robject{BM.PB}), we can do this by

<<adjust>>=
ga.cc.BMPB <- GlobalAncova(xx=exprs(golubX), group=gr, covars=golubX$BM.PB, 
test.genes=cellcycle, method="both", perm=1000) 
ga.cc.BMPB
@


With the more general function call we would simply adjust the definitions of model formulas, namely 
\Rfunarg{formula.full =} $\sim$ \Rfunarg{ALL.AML + BM.PB} and \Rfunarg{formula.red =} $\sim$ \Rfunarg{BM.PB}.

%The source of the samples apparently has some explanatory effect on the outcome resulting in a smaller F--ratio than in the model without adjusting. 
%But the influence of the genes is still highly significant.
The source of the samples does not seem to have an explanatory effect on the outcome since F--statistics and p--values are very similar to the model without adjustment. 

With the \Rpackage{globaltest} we get a higher p--value. 

<<gt.adjust, results=hide>>=
gt.cc.BMPB <- globaltest(X=golubX, Y=ALL.AML ~ BM.PB, genesets=cellcycle)
gt.cc.BMPB
@
<<echo=F>>=
gt.cc.BMPB
@

Permutation based p--values can also be obtained with Goeman's test, however only when covariates are absent. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\subsection{van t'Veer Data and p53--Signalling Pathway}
\subsection{van't Veer Data and p53--Signalling Pathway} \label{vantVeer}

<<echo=F>>=
data(vantVeer)
data(phenodata)	
data(pathways)	
@

We present another example from a study on breast cancer from van't Veer et al. (2002) \cite{vV:02}. 
This example illustrates how more complex tasks than comparing %expression levels between 
just two clinical groups can be performed with \Rpackage{GlobalAncova}.
%The authors present a set of 70 genes as a signature to predict cancer recurrence.
A subset of the data consisting of the expression values for \Sexpr{dim(phenodata)[1]} patients without {\it BRCA1} or {\it BRCA2} 
mutations is available with the package. The dataset (\Robject{vantVeer}) is restricted to \Sexpr{dim(vantVeer)[1]} genes associated 
with \Sexpr{length(pathways)} cancer related pathways that are provided as a list named (\Robject{pathways}), too. 
We take one gene from the original data additionally to the expression set, namely
'AL137718'. This gene is part of the original van't Veer prognosis signature. We will later use it to demonstrate
how signature genes can be related to pathways.
Information about 
some of the originally surveyed covariates is stored in \Robject{phenodata}. %Further we provide a list of \Sexpr{length(pathways)} cancer related pathways (\Robject{pathways}).
%This example illustrates the role of covariates in the context of global testing in more detail.
The tumour suppressor protein p53 contributes as a transcription factor to cell cycle arrest and apoptosis induction.
%Oncogene. 2004 Apr 22;23(19):3376-84. 	
%Identification of Tcf-4 as a transcriptional target of p53 signalling.
%Rother K, Johne C, Spiesbach K, Haugwitz U, Tschop K, Wasner M, 
%Klein-Hitpass L, Moroy T, Mossner J, Engeland K.
Therefore, first the p53-signalling pathway is selected as a candidate, where differential expression
between relevant prognostic groups, defined by the development of distant metastases within five years, was expected. 

<<vV>>=
data(vantVeer)
data(phenodata)	
data(pathways)
metastases <- phenodata$metastases
p53 <- pathways$p53_signalling	
@

%First we compute the Global Ancova without adjusting for covariates and get a 
We get a significant result with the global ANCOVA. 
 
<<p53test>>=	
set.seed(111)
vV.1 <- GlobalAncova(xx=vantVeer, group=metastases, test.genes=p53, method="both", perm=1000)
vV.1
@


\subsubsection{Analysis of Various Clinical Groups}

In the new version of the package also clinical variables with more than two groups can be considered. For demonstration we 
investigate differential expression for the three ordered levels of tumour grade.

<<grade>>=
set.seed(111)
vV.3 <- GlobalAncova(xx=vantVeer, formula.full=~grade, formula.red=~1, model.dat=phenodata, test.genes=p53, method="both", perm=1000)
vV.3
@ 

\subsubsection{Gene--Gene Interaction}

Now we want to go into the matter of other interesting biological questions. For example one might ask if there exists interaction 
between the expression of genes which the authors in \cite{vV:02} presented as signature for prediction of cancer recurrence and the 
expression of genes in a certain pathway. This question can be answered by viewing the expression values of the signature genes as linear 
regressors and to test their effects on the expression pattern of the pathway genes. For demonstration we pick the signature gene "AL137718",
which is not part of any of the pathways,
and test its effect on the p53--signalling pathway. 
Assume that we also want to adjust for the Estrogen receptor status. The analysis can be carried out 
in the following way.

<<interact>>=
signature.gene <- "AL137718"
model <- data.frame(phenodata, signature.gene=vantVeer[signature.gene, ])
set.seed(111)
vV.4 <- GlobalAncova(xx=vantVeer, formula.full=~signature.gene + ERstatus, formula.red=~ERstatus, model.dat=model, test.genes=p53, method="both", perm=1000)
vV.4
@

Assuming a significance level of 0.05 we get a significant effect of the signature gene on the p53--signalling pathway.

\subsubsection{Co--Expression}

Next we want to analyse co--expression regarding the clinical outcome of building distant metastases within five years. This can be done by simply adding the variable \Robject{metastases} to the full and reduced model, respectively. Such layout corresponds to testing the linear effect of the signature gene stratified not only by Estrogen receptor status but also by \Robject{metastases}. 

<<coexpr>>=
set.seed(111)
vV.5 <- GlobalAncova(xx=vantVeer, formula.full=~metastases + signature.gene + ERstatus, formula.red=~metastases + ERstatus, model.dat=model, test.genes=p53, method="both", perm=1000)
vV.5
@

Again we get a significant result.

Supposably the most interesting question in this case concerns differential co--expression. Differential co--expression is on hand if the effect of the signature gene behaves different in both \Robject{metastases} groups. In a one dimensional context this would become manifest by different slopes of the regression lines. Hence what we have to test is the interaction between \Robject{metastases} and \Robject{signature.gene}.

<<diffcoexpr>>=
set.seed(111)
vV.6 <- GlobalAncova(xx=vantVeer, formula.full=~metastases * signature.gene + ERstatus, formula.red=~metastases + signature.gene + ERstatus, model.dat=model, test.genes=p53, method="both", perm=1000)
vV.6
@

We observe a significant differential co-expression between the chosen signature gene and the p53--signalling pathway.

With \Rpackage{globaltest} we can also test gene-gene interaction, also adjusted for phenotype groups. But it is
not possible to test for differential co-expression or the influence of more than one signature gene on a
pathway. On the other hand globaltest is able to deal with survival times as clinical outcome.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Testing Several Pathways Simultaneously}

Systems biology involves the study of mechanisms underlying complex biological processes as integrated systems of many diverse interacting components, often referred to as pathways. 

We regard the possibility to investigate differential gene expression simultaneously for several of those pathways as a contribution towards understanding biological relevant relations. 

The user can apply \Rfunction{GlobalAncova} to compute p--values for a couple of pathways with one call by specifying the \Rfunarg{test.genes} option.  The members of each pathway to be tested must belong to genes in the expression--matrix.
Afterwards a suitable correction for multiple testing has to be applied. An alternative based on the closed testing approach is described later.

%Suppose for example that we want to test the \Sexpr{length(pathways)} cancer related pathways with the van t'Veer data. 
Suppose for example that for sake of simplicity we want to test the first four of the cancer related pathways with the van't Veer data. 
We proceed as follows.

<<pathways>>=
metastases <- phenodata$metastases
ga.pw <- GlobalAncova(xx=vantVeer, group=metastases, test.genes=pathways[1:4], method="both", perm=1000) 
ga.pw
@

%Suppose for example that we want to test the first five KEGG pathways with the Golub data. We proceed as follows.
%
%<<kegg>>=
%ga.kegg <- GlobalAncova(exprs(golubX), gr, test.genes=kegg[1:5], 
%perm=1000) 
%ga.kegg
%@

The result is a matrix whose rows correspond to the different pathways. 

With the \Rfunction{globaltest} we get a similar matrix.

<<gt.pathways>>=
gt.pw <- globaltest(X=vantVeer, Y=metastases, genesets=pathways[1:4])
gt.pw
@


\subsection{Simultaneous Adjustment of p--values}

Next we show how to extract p--values for correction for multiple testing.
Note however that due to the extremely high correlations between
these tests, many procedures that correct for multiple testing
here are inappropriate. 
An appropriate way of adjusting would be for example the method of Holm, 1979 \cite{Holm:79}.
An alternative to such adjustments that is not affected by correlations between tests is a closed
testing procedure. 
For this approach you need a family of null hypotheses that is closed
under intersection. Then a single hypothesis can be rejected at level $\alpha$ if it
is rejected along with all hypotheses included in it (\cite{Mar:76}). %(Marcus et al. 1976). 

For the adjustment according to Bonferroni and Holm we build a vector of
the raw p--values. The function \Rfunction{mt.rawp2adjp} provides several adjusting
methods. We here display only the raw and ``Holm'' adjusted p--values. 
To obtain the original order of the pathways we order the result of \Rfunction{mt.rawp2adjp} according to \Robject{index}.

<<>>=
ga.pw.raw <- ga.pw[ ,"p.perm"]
ga.pw.adj <- mt.rawp2adjp(ga.pw.raw)
ga.result <- ga.pw.adj$adjp[order(ga.pw.adj$index), c("rawp", "Holm")]
rownames(ga.result) <- names(pathways)[1:4]
ga.result             
gt.pw.raw <- p.value(gt.pw)   
gt.pw.adj <- mt.rawp2adjp(gt.pw.raw)
gt.result <- gt.pw.adj$adjp[order(gt.pw.adj$index), c("rawp", "Holm")]
rownames(gt.result) <- names(pathways)[1:4]             
gt.result
@

Allowing a family--wise error rate of 0.05 %with \Rfunction{GlobalAncova} all but two pathways remain significant.
%The \Rfunction{globaltest} appears to be more conservative in this case resulting in only three significant pathways.
all but one pathways remain significant for both methods.


\subsection{Closed Testing Procedure} \label{closed}

Closed testing procedures (\cite{Mar:76}) %(Marcus et al., 1976 \cite{Mar:76}) 
offer a versatile and powerful approach to the multiple testing problem. Implementation is non--trivial, therefore, the program given in this version should be regarded as a prototype.  

In order to apply the closed testing procedure we first have to create the required family of
hypotheses by building all intersections between the ``natural'' hypotheses
tested above and all intersections of those new hypotheses and so on. 

The resulting family of hypotheses can be illustrated in a directed graph. 
If we just for the sake of illustration assume that we have only four hypotheses named ``1'', $\ldots$ ``4'' then
the node ``1-2-3-4'' for example stands for the global hypothesis that the genes
of all four pathways are not differentially expressed.
Now the interesting hypothesis ``1'' for example can be rejected if also the hypotheses
``1-2-3-4'', ``1-2-3'', ``1-2-4'', ``1-3-4'', ``2-3-4'', ``1-2'', $\ldots$, ``1-4'' are rejected.
These relationships are represented by the edges of the graph.
%Note that the required package \Rpackage{Rgraphviz} for plotting the graph
%currently only is available for Unix/Linux.

<<hypoth, echo=F, fig=T>>=
if(require(Rgraphviz))
{
res <- GlobalAncova:::Hnull.family(pathways[1:4])
graph <- new("graphNEL", nodes=names(res), edgemode="directed")
graph <- addEdge(from=c(rep(names(res)[15],4),rep(names(res)[10],3),rep(names(res)[11],3),
  rep(names(res)[13],3),rep(names(res)[14],3),rep(names(res)[5],2),rep(names(res)[6],2),
  rep(names(res)[7],2),rep(names(res)[8],2),rep(names(res)[9],2),rep(names(res)[12],2)), 
  to=c(names(res)[10],names(res)[11],names(res)[13],names(res)[14],names(res)[5],names(res)[6],
  names(res)[8],names(res)[5],names(res)[7],names(res)[9],names(res)[6],names(res)[7],names(res)[12],
  names(res)[8],names(res)[9],names(res)[12],names(res)[1],names(res)[2],names(res)[1],names(res)[3],
  names(res)[1],names(res)[4],names(res)[2],names(res)[3],names(res)[2],names(res)[4],names(res)[3],
  names(res)[4]), graph, weights=rep(1, 28))

att <- list()
lab <- c("1","2","3","4","1-2","1-3","1-4","2-3","2-4","1-2-3","1-2-4","3-4","1-3-4","2-3-4","1-2-3-4")
names(lab) <- names(res)
att$label <- lab

plot(graph, nodeAttrs=att)
} else {
plot(1, 1, type="n", main="This plot requires Rgraphviz", 
     axes=FALSE, xlab="", ylab="")
}
@

We can compute the closed testing procedure using the function

<<closed.test>>=
ga.closed <- GlobalAncova.closed(xx=vantVeer, group=metastases, 
test.genes=pathways[1:4], previous.test=ga.pw, level=0.05, method="approx")
@

%<<closed.test>>=
%ga.closed <- GlobalAncova.closed(xx=exprs(golubX), group=gr, 
%test.genes=kegg[1:5], previous.test=ga.kegg, level=0.05, perm=100)
%@

where \Rfunarg{test.genes} is again a list of pathways. In order to shorten computing time we can provide
the results of the previous application of \Rfunction{GlobalAncova} for the pathways of interest.
The option \Rfunarg{level} allows to manipulate the level of significance. 
There is again the possibility to choose between permutation and asymptotic p-values via the option \Rfunarg{method}. %\Rfunarg{perm} again gives the desired number of permutations used in the permutation test. 
Note that if you provide results of previous computation, the type of p-values has to correspond, i.e. if we now want to use \Rfunarg{method = "approx"} in the previous test we should have used \Rfunarg{method = "approx"} or \Rfunarg{method = "both"} such that asymptotic p-values are available.

Also for \Rfunction{GlobalAncova.closed} all three different function calls as for \Rfunction{GlobalAncova} itself are possible. 

The function \Rfunction{GlobalAncova.closed} provides the formed null hypotheses 
(this means lists of genes to be tested simultaneously), the test results for each pathway 
of interest and the names of significant and non significant pathways.
Names for the intersections of hypotheses are built by simply coercing the names of the respective pathways.
If for a pathway one single hypothesis can not be rejected there is no need to test all the 
remaining hypotheses. That is why in test results of non significant pathways lines 
are filled with NA's after a p--value $> \alpha$ occured. Here only test results for the first pathway are displayed.

<<names.closed.test>>=
names(ga.closed)
rownames(ga.closed$test.results[[1]])
rownames(ga.closed$test.results[[1]]) <- NULL
ga.closed$test.results[1]
ga.closed$significant
ga.closed$not.significant
@

We get the same significant and non significant pathways as before.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Multiple testing on the Gene Ontology graph}

When testing gene sets defined by the Gene Ontology it is of special interest to incorporate the hierarchical structure of the GO graph. Goeman and Mansmann (2007)\cite{GoeMan:07} developed the \emph{focus level} method, a multiple testing approach on the GO that combines the correction method of Holm (1979) \cite{Holm:79} and the closed testing procedure from Marcus et al. (1976) \cite{Mar:76} (also used in section \ref{closed}). The method is originally implemented in package \Rpackage{globaltest}. We adapted the corresponding functions such that the procedure now is available also with \Rfunction{GlobalAncova}.
For details see the vignette of \Rpackage{globaltest}.  

First, the GO graph must be prepared for the data set at hand. We use again the data of Golub et al.

<<makeGOstructure>>=
bp <- makeGOstructure(golubX, "hu6800")
@

This creates a \Rclass{GOstructure} object (class from \Rpackage{globaltest}) for the Golub data, using
the annotation package \Rpackage{hu6800}.

The focus level procedure requires definition of a
\emph{focus level}. This is a set of GO terms reflecting the level of specificity that is of most interest. The focus level can be specified as an arbitrary vector of GO identifiers. However, there is also an automatic way for defining the focus level via the function \Rfunction{getFocus}.

<<getFocus>>=
focusBP <- getFocus(bp, maxatoms=5)
str(focusBP)
@

The argument \Rfunarg{maxatoms} determines the maximum
complexity of the subgraphs of all descendants of focus level nodes.
We choose a lower value (default 10) which leads to more specific focus level terms and less time consuming processing.

The significant subgraph of predefined size (\Rfunarg{stopafter = 30}) can be calculated with

<<gtGO, results=hide>>=
go30 <- GAGO(exprs(golubX), group=pData(golubX)$ALL.AML, focus = focusBP, GO = bp, stopafter = 30)
@
<<>>=
str(go30)
@

All arguments for specifying the linear model used in \Rfunction{GlobalAncova} can be given here. Only the parameter 
\Rfunarg{method} is not available because the focus level procedure does only work with the asymptotic test. Note however, that 
still a number of permutations can be specified (\Rfunarg{perm}, default 10,000) since very large GO terms (with more annotated genes 
than defined by parameter \Rfunarg{max.group.size}, default 2500) are tested permutation based. 

See the \Rpackage{globaltest} vignette for examples how to visualize the results of the focus level procedure.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\setkeys{Gin}{width=0.6\textwidth}

\section{Diagnostic Plots} \label{plots}

There are two types of diagnostic plots available supporting communication and interpretation of results of the global ANCOVA.
The \Rfunction{Plot.genes} visualizes the 
influence of individual genes on the test result while the
\Rfunction{Plot.subjects} visualizes the influence of individual 
samples. Both plots are based on the decomposition of sums of squares.

%We use again the Golub data for demonstration of the plot functions.
We use again the van't Veer data constricted to the genes of the p53--signalling pathway for demonstration of the plot functions.


\subsection{Gene Plot}

The influence of each gene on the outcome of the test can be assessed and visualized with a diagnostic plot generated by our 
function \Rfunction{Plot.genes}.
It corresponds to the function \Rfunction{geneplot} in the \Rpackage{globaltest} package.
The function \Rfunction{Plot.genes} gives a graphical display of single gene-wise analysis for all genes. Bars are always positive as 
a reduction of sum of squares is always achieved in this case.
The bar height indicates the influence of the respective gene on the test statistic. 
The added reference line is the residual mean square error per gene and corresponds to the expected height of the 
bars under the null hypothesis which says that the gene is not associated with the clinical outcome. 
The actual and expected bar heights also correspond to the nominator and denominator of gene-wise F-statistics. Hence the ratio of the 
two values is a measure for the association of the respective gene with the phenotype.
Bar heights for all genes can be returned by setting the option \Rfunarg{returnValues} to \Robject{TRUE}. This helps to detect genes 
with most influence on the global statistics.
Note however that comparisons between different gene groups can not easily be done by means of these values directly since 
different group sizes have an impact on global significance.  
 
%Covariate information can be included in the same way as in the \Rfunction{GlobalAncova} function with the \Rfunarg{covars} option.
The bars can be colored according to a variable of interest with the option \Rfunarg{Colorgroup} in order to show in which of the groups 
a gene has the highest expression values. 
The automatically chosen bar labels can be manipulated with the parameter \Rfunarg{bar.names}.

The commands for creating gene plots in the \Rpackage{GlobalAncova} and
the \Rpackage{globaltest} are as follows. Note that for the former one again three alternatives for function calls are provided, see 
section \ref{onepathway} for details.   
%For facility of inspection it is useful not to plot the bars for all genes at one time
%but only for a few, for example 40. 

The two approaches show almost the same results (figures \ref{ga.genepl1} and \ref{gt.genepl1}). We prefer plotting horizontal bars 
rather than vertical because we think it is easier to read off the bar heights this way.
%In the group variable \Rfunarg{gr} 0 represents AML and 1 the ALL patients.

<<geneplot>>=
Plot.genes(xx=vantVeer, group=metastases, test.genes=p53)
gt.vV <- globaltest(X=vantVeer, Y=metastases, genesets=p53)
geneplot(gt.vV)
@

\begin{figure}[htb!]
\begin{center}
<<genes, fig=TRUE, echo=F, width=6, height=8>>=
Plot.genes(xx=vantVeer, group=metastases, test.genes=p53)
@
\vspace{-0.4cm}
\caption{{\small \label{ga.genepl1} Gene Plot for the van't Veer data with \Rpackage{GlobalAncova}. 
Shown are the genes of the p53--signalling pathway. The bar height indicates the influence of the respective gene on the test 
statistic. The color shows in which of the phenotype groups
the gene has higher expression values. The reference line is the residual mean square 
error per gene.}}
\end{center}
\end{figure}

\begin{figure}[htb!]
\begin{center}
<<gt_genes, fig=TRUE, echo=F>>=
geneplot(gt.vV)
@
\vspace{-0.4cm}
\caption{{\small \label{gt.genepl1} Gene Plot for the van't Veer data with \Rpackage{globaltest}. 
Shown are the genes of the p53--signalling pathway. The bar height indicates the influence of the respective gene on the test 
statistic. The colour shows in which of the phenotype groups
the gene has higher expression values. The reference line gives the expected height of the 
bar under the null hypothesis. Marks indicate with how many standard deviations the bar 
exceeds the reference line.}}
\end{center}
\end{figure}

%<<geneplot>>=
%Plot.genes(exprs(golubX)[cellcycle[1:40], ], gr)
%gp.cc <- geneplot(gt.cc)
%plot(gp.cc[1:40])
%@

%\begin{figure}[htb!]
%\begin{center}
%<<genes, fig=TRUE, echo=F>>=
%Plot.genes(exprs(golubX)[cellcycle[1:40], ], gr)
%@
%\vspace{-0.4cm}
%\caption{{\small Gene Plot for the Golub data with \Rpackage{GlobalAncova}. 
%Shown are the first 40 genes of the cell cycle pathway. The bar height indicates the influence of the respective gene on the test 
%statistic. The colour shows in which of the phenotype groups
%the gene has higher expression values. The reference line is the residual mean square 
%error per gene.}}
%\end{center}
%\end{figure}
%
%\begin{figure}[htb!]
%\begin{center}
%<<gt_genes, fig=TRUE, echo=F>>=
%plot(gp.cc[1:40])
%@
%\vspace{-0.4cm}
%\caption{{\small Gene Plot for the Golub data with \Rpackage{globaltest}. 
%Shown are the first 40 genes of the cell cycle pathway. The bar height indicates the influence of the respective gene on the test 
%statistic. The colour shows in which of the phenotype groups
%the gene has higher expression values. The reference line gives the expected height of the 
%bar under the null hypothesis. Marks indicate with how many standard deviations the bar 
%exceeds the reference line.}}
%\end{center}
%\end{figure}

In this case where only the influence of one variable is of interest (and therefore the easiest version of possible function calls is 
chosen), the same variable is assumed to be relevant for coloring. However one is free to specify another coloring. For example for 
the same plot we could ask which genes are higher expressed in samples with either positive or negative Estrogen receptor status, see 
figure \ref{ga.genepl2}.

<<geneplot>>=
Plot.genes(xx=vantVeer, formula.full=~metastases, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="ERstatus")
@

\begin{figure}[htb!]
\begin{center}
<<genes2, fig=TRUE, echo=F, width=6, height=8>>=
Plot.genes(xx=vantVeer, formula.full=~metastases, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="ERstatus")
@
\vspace{-0.4cm}
\caption{{\small \label{ga.genepl2} Gene Plot for the van't Veer data with \Rpackage{GlobalAncova}. 
Shown are the genes of the p53--signalling pathway. The bar height indicates the influence of the respective gene on the test 
statistic. The color shows in which of the specified phenotype groups, in this case Estrogen receptor status,
the gene has higher expression values. The reference line is the residual mean square 
error per gene.}}
\end{center}
\end{figure}
   


\subsection{Subjects Plot}

The function \Rfunction{Plot.subjects}
visualizes the influence of the individual samples on the test result
and corresponds to the \Rfunction{sampleplot} of Goeman.
The function \Rfunction{Plot.subjects} gives information on the reduction of sum of squares per subject. Here we sum over genes. Large reduction demonstrates a good approximation of a subject's gene expressions by the corresponding group means. If an individual does not fit into the pattern of its phenotype, negative values can occur.   
A small p--value will 
therefore generally coincide with many positive bars. If there 
are still tall negative bars, these indicate deviating samples: 
removing a sample with a negative bar would result in a lower 
p-value. 
%Again we can use \Rfunarg{covars} for covariate adjustment.
The bars are colored to distinguish samples 
of different clinical entities that can again be specified by the user through the option \Rfunarg{Colorgroup}. With the option 
\Rfunarg{sort} it is also possible to sort the bars with respect to the phenotype groups. 
Bar labels can be changed with the argument \Rfunarg{bar.names}.
Also in the subjects plot bar heights can be returned by setting the option \Rfunarg{returnValues} to \Robject{TRUE}. 
That may help to detect, not only visually, samples which do not fit into their respective clinical groups.
%Before plotting we add the sample names to the expression matrix. 
%Otherwise with \Rfunction{Plot.subjects} samples would just
%be enumerated from 1 to \Sexpr{length(pData(golubX))}.

We compare again the different approaches (figures \ref{ga.subjects1} and \ref{gt.subjects1}):

<<subjectsplot>>=
#colnames(exprs(golubX)) <- pData(golubX)[ ,1]
Plot.subjects(xx=vantVeer, group=metastases, test.genes=p53, legendpos="bottomright")
sampleplot(gt.vV)
@

\begin{figure}[htb!]
\begin{center}
<<subjects, fig=TRUE, echo=F, width=6, height=14>>=
Plot.subjects(xx=vantVeer, group=metastases, test.genes=p53, legendpos="bottomright")
@
\vspace{-0.4cm}
\caption{{\small \label{ga.subjects1} Subjects Plot for the van't Veer data with \Rpackage{GlobalAncova}. The bar height 
indicates the influence of the respective sample on the test result. If an individual does 
not fit into the pattern of its phenotype, negative values can occur. Bars are colored 
corresponding to phenotype groups.}}
\end{center}
\end{figure}

\begin{figure}[htb!]
\begin{center}
<<gt_subjects, fig=TRUE, echo=F>>=
sampleplot(gt.vV)
@
\vspace{-0.4cm}
\caption{{\small \label{gt.subjects1} Subjects Plot for the van't Veer data with \Rpackage{globaltest}. The bar height 
indicates the influence of the respective sample on the test result. If an individual does 
not fit into the pattern of its phenotype, negative values can occur. Bars are colored 
corresponding to groups. The reference line shows the expected influence of the samples 
under the null hypothesis. Marks on the bars indicate the standard deviation of the influence 
of the sample under the null hypothesis.}}
\end{center}   
\end{figure}

%<<subjectsplot>>=
%colnames(exprs(golubX)) <- pData(golubX)[ ,1]
%Plot.subjects(exprs(golubX), gr)
%sampleplot(gt.cc)
%@

%\begin{figure}[htb!]
%\begin{center}
%<<subjects, fig=TRUE, echo=F>>=
%Plot.subjects(exprs(golubX), gr)
%@
%\vspace{-0.4cm}
%\caption{{\small Subjects Plot for the Golub data with \Rpackage{GlobalAncova}. The bar height 
%indicates the influence of the respective sample on the test result. If an individual does 
%not fit into the pattern of its phenotype, negative values can occur. Bars are coloured 
%corresponding to groups.}}
%\end{center}
%\end{figure}

%\clearpage

%\begin{figure}[htb!]
%\begin{center}
%<<gt_subjects, fig=TRUE, echo=F>>=
%sampleplot(gt.cc)
%@
%\vspace{-0.4cm}
%\caption{{\small Subjects Plot for the Golub data with \Rpackage{globaltest}. The bar height 
%indicates the influence of the respective sample on the test result. If an individual does 
%not fit into the pattern of its phenotype, negative values can occur. Bars are coloured 
%corresponding to groups. The reference line shows the expected influence of the samples 
%under the null hypothesis. Marks on the bars indicate the standard deviation of the influence 
%of the sample under the null hypothesis.}}
%\end{center}   
%\end{figure}

The function \Rfunction{Plot.subjects} can be invoked by the three alternative function calls (see section \ref{onepathway}) and hence 
also plots corresponding to more complex testing challenges can be produced as well. To give just one example we consider again the 
influence of the tumour grade, which can take three possible values, on gene expression (figure \ref{ga.subjects2}).

<<subjectsplot2>>=
Plot.subjects(xx=vantVeer, formula.full=~grade, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="grade", legendpos="topleft")
@

\begin{figure}[htb!]
\begin{center}
<<subjects2, fig=TRUE, echo=F, width=6, height=14>>=
Plot.subjects(xx=vantVeer, formula.full=~grade, formula.red=~1, model.dat=phenodata, test.genes=p53, Colorgroup="grade", legendpos="topleft")
@
\vspace{-0.4cm}
\caption{{\small \label{ga.subjects2} Subjects Plot for the van't Veer data with \Rpackage{GlobalAncova}. Tumour grade is the clinical variable of interest. The bar height 
indicates the influence of the respective sample on the test result. If an individual does 
not fit into the pattern of its phenotype, negative values can occur. Bars are colored 
corresponding to phenotype groups.}}
\end{center}
\end{figure}      


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\clearpage

\section{Acknowledgements}

We thank Sven Knüppel who took the initiative and contributed a C-code
implementation of the permutation test to our package.\\
This work was supported by the NGFN project 01 GR 0459, BMBF, Germany.


\bibliographystyle{plain}
\bibliography{references}

\end{document}

