% \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 \and Sven Knueppel}


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

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


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


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

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.

\vspace{0.2cm}
\noindent {\bf 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, test.genes=NULL)
@ 

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.  
The number of permutations is 10,000 for default but in the current version the test for many genes can take quite a long time with 10,000 permutations. 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.

\vspace{0.2cm}
\noindent {\bf 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. 

\vspace{0.2cm}
\noindent {\bf 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 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} \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. 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
@


\vspace{0.2cm}
\noindent {\bf Analysis of Arbitrary Clinical Variables}\\
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
@ 

\vspace{0.2cm}
\noindent {\bf 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" and again 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.

\vspace{0.2cm}
\noindent {\bf 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.

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


\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
@

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 <- gt.pw@res[ ,"P-value"]
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.

<<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")
@

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

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

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}

This work was supported by the NGFN project 01 GR 0459, BMBF, Germany.


\bibliographystyle{plain}
\bibliography{references}

\end{document}

