% -*- mode: noweb; noweb-default-code-mode: R-mode; -*-
%\VignetteIndexEntry{SPIA}
%\VignetteKeywords{Pathway Analysis}
%\VignettePackage{SPIA}
\documentclass[11pt]{article}

%\usepackage{amsmath,epsfig,psfig,fullpage} 
\usepackage{amsmath,epsfig,fullpage} 
%\usepackage{graphicx,pstricks}
%\usepackage{ifpdf}
\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}
\usepackage{url}

\parindent 0in

\bibliographystyle{abbrvnat}
\begin{document}

\title{\bf Bioconductor's SPIA package}
\author{Adi L. Tarca$^{1,2,3}$, Purvesh Khatri$^{1}$ and Sorin Draghici$^{1}$}

\maketitle

$^1$Department of Computer Science, Wayne State University\\
$^2$Bioinformatics and Computational Biology Unit of the NIH Perinatology Research Branch\\
$^3$Center for Molecular Medicine and Genetics, Wayne State University \\


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

\section{Overview}

This package implements the Signaling Pathway Impact Analysis (SPIA) algorithm described in \cite{TarcaSPIA:2008}, \cite{Khatri:2007a} and \cite{DraghiciPE:2007}.
SPIA uses the information from a set of differentially expressed genes and  
their fold changes, as well as pathways topology in order to assess the significance of the pathways in the condition under the study.
The current version of SPIA algorithm uses KEGG signaling pathway data. SPIA ready KEGG pathway data for homo sapiens is included in the package and also 
available at \\
 \url{http://bioinformaticsprb.med.wayne.edu/SPIA/}.\\
The pathways included for each organism are those containing only directed relations between genes/proteins and no reactions.\\
The KEGG data that was preprocessed for SPIA analysis was downloaded from KEGG's ft repository on: 09/22/2010.


\section{Pathway analysis with SPIA package}

This document provides basic introduction on how to use the {\tt SPIA} package. For extended description of the methods used by 
this package please consult these references: \cite{ TarcaSPIA:2008,Khatri:2007a,DraghiciPE:2007}.\\ 

We demonstrate the functionality of this package using a colorectal cancer dataset obtained using Affymetrix GeneChip technology and available 
through GEO (GSE4107). The experiment contains 10 normal samples and 12 colorectal cancer samples and is described by \cite{Hong:2007}. 
RMA preprocessing of the raw data was performed using the {\tt affy} package, and a two group moderated t-test was applied using the {\tt limma} package. The data frame obtained as an end result 
from the function topTable in limma is used as starting point for preparing the input data for SPIA. This data frame called {\tt top} was made available in the 
{\tt colorectalcancer} dataset included in the SPIA package:
 
<<eval=TRUE, echo=TRUE>>=
library(SPIA)
data(colorectalcancer)
options(digits=3)
head(top)
@

For SPIA to work, we need a vector with log2 fold changes between the two groups for all the genes considered to be differentially expressed.
The names of this vector must be Entrez gene IDs.
The following lines will add one additional column in the {\tt top} data frame annotating each affymetrix probeset to an Entrez ID. 
Since there may be several probesets for the same Entrez ID, there are two easy ways to obtain one log fold change per gene. The first option is to use
the fold change of the most significant probeset for each gene, while the second option is to average the log fold-changes of all probestes of the same gene.
In the example below we used the former approach. The genes in this example are called differentially expressed provided that their FDR p-value is less than 0.05.
The following lines start with the {\tt top} data frame and produce two vectors that are required as input by {\tt spia} function:
 
<<eval=TRUE, echo=TRUE>>=
library(hgu133plus2.db)
x <- hgu133plus2ENTREZID 
top$ENTREZ<-unlist(as.list(x[top$ID]))
top<-top[!is.na(top$ENTREZ),]
top<-top[!duplicated(top$ENTREZ),]
tg1<-top[top$adj.P.Val<0.05,]
DE_Colorectal=tg1$logFC
names(DE_Colorectal)<-as.vector(tg1$ENTREZ)
ALL_Colorectal=top$ENTREZ
@

The {\tt DE\_Colorectal} is a vector containing the log2 fold changes of the genes found to be 
differentially expressed between cancer and normal samples, and {\tt ALL\_Colorectal} is a vector 
with the Entrez IDs of all genes profiled  on the microarray. The names of the {\tt DE\_Colorectal} are the Entrez gene IDs corresponding 
to the computed log fold-changes.  

<<eval=TRUE,echo=TRUE,results=verbatim>>=
DE_Colorectal[1:10]
ALL_Colorectal[1:10]
@


The SPIA algorithm takes as input the two vectors above and produces a table of pathways ranked from the most to the least significant. 
This can be achieved by calling the {\tt spia} function as follows:

<<eval=TRUE, echo=TRUE>>=
# pathway analysis based on combined evidence; # use nB=2000 or more for more accurate results
res=spia(de=DE_Colorectal,all=ALL_Colorectal,organism="hsa",nB=2000,plots=FALSE,beta=NULL)
#make the output fit this screen
res$Name=substr(res$Name,1,10)
#show first 15 pathways, ommit KEGG links
res[1:15,-12]
@


\begin{figure}
\begin{center}
\includegraphics{PFs}
\end{center}
\caption{Perturbations plot for colorectal cancer pathway (KEGG ID hsa:05210) using the {\tt colorectalcancer} dataset. The perturbation of all genes 
in the pathway are shown as a function of their initial log2 fold changes (left panel).
Non DE genes are assigned 0 log2 fold-change. The null distribution of the net accumulated perturbations is also given (right panel). The observed net accumulation tA with the real data is shown as a red vertical line.}
\label{fig:PFs}
\end{figure}

If the {\tt plots} argument is set to {\tt TRUE} in the function call above, a plot like the one shown in Figure ~\ref{fig:PFs} 
is produced for each pathway  on which there are differentially expressed genes. These plots are saved in a pdf file in the current directory.

An overall picture of the pathways significance according to both the over-representation evidence and perturbations based evidence 
can be obtained with the function {\tt plotP} and shown in Figure ~\ref{fig:pPlot}. The Colorectal cancer pathway is shown in green.

\begin{figure}[htbp]
\label{fig:pPlot}
  \begin{center}
<<fig=TRUE>>=
plotP(res,threshold=0.05)
points(I(-log(pPERT))~I(-log(pNDE)),data=res[res$ID=="05210",],col="green",pch=19,cex=1.5)

@
\caption{SPIA evidence plot for the colorectal cancer dataset. Each pathway is represented by one dot. The pathways at the right of the red oblique line are significant after 
Bonferroni correction of the global p-values, pG. The pathways at the right of the blue oblique line are significant after 
a FDR correction of the global p-values, pG.}
\end{center}
\end{figure}


In this plot, the horizontal axis represents the p-value (minus log of) corresponding to the probability of obtaining at least the observed
number of genes (NDE) on the given pathway just by chance. The vertical axis represents the p-value (minus log of) corresponding to
the probability of obtaining the observed total accumulation (tA) or more extreme on the given pathway just by chance. The computation of pPERT
is described in \cite{TarcaSPIA:2008}. 
In Figure ~\ref{fig:pPlot} each pathway is shown as a bullet point, and those significant at 5\% (set by the {\tt threshold} argument in {\tt plotP}) 
after Bonferroni correction are shown in red.\\

SPIA algorithm is illustrated also using the Vessels dataset:

<<eval=TRUE, echo=TRUE>>=
data(Vessels)
# pathway analysis based on combined evidence; # use nB=2000 or more for more accurate results
res<-spia(de=DE_Vessels,all=ALL_Vessels,organism="hsa",nB=500,plots=FALSE,beta=NULL,verbose=FALSE)
#make the output fit this screen
res$Name=substr(res$Name,1,10)
#show first 15 pathways, ommit KEGG links
res[1:15,-12]
@

The pathway image as provided by KEGG having the differentially expressed genes highlighted in red can be obtained by pasting
in a web browser the links available in the KEGGLINK column of the data frame produced by the function spia. For example, 

<<eval=TRUE, echo=TRUE>>=
res[,"KEGGLINK"][20]
@

is the link that would display the image of the 20th pathway in the res dataframe above.   

Note that the results for these datasets my differ from the ones described in \cite{TarcaSPIA:2008} since a) the pathways database used herein 
was updated and b) the default beta values were changed. 

The directed adjacency matrices of the graphs describing the different types of relations between genes/proteins (such as activation or repression)
used by SPIA are available in the {\tt extdata/hsaSPIA.RData} file for the homo sapiens organism.
The types of relations considered by SPIA and the default weight (beta coefficient) given to them are:
 
<<eval=TRUE, echo=TRUE>>=
rel<-c("activation","compound","binding/association","expression","inhibition","activation_phosphorylation","phosphorylation",
"indirect","inhibition_phosphorylation","dephosphorylation_inhibition","dissociation","dephosphorylation","activation_dephosphorylation",
"state","activation_indirect","inhibition_ubiquination","ubiquination","expression_indirect","indirect_inhibition","repression",
"binding/association_phosphorylation","dissociation_phosphorylation","indirect_phosphorylation")
beta=c(1,0,0,1,-1,1,0,0,-1,-1,0,0,1,0,1,-1,0,1,-1,-1,0,0,0)
names(beta)<-rel
cbind(beta)
@

A 0 value for a given relation type results in discarding those type of relations from the analysis for all pathways. The default values of {\tt beta} can changed by the user at any
time by setting the {\tt beta} argument of the {\tt spia} function call.    

Other organisms' KEGG pathway data can be downloaded from \url{http://bioinformaticsprb.med.wayne.edu/SPIA} as a "[org]SPIA.RData" file 
and copied into the {\tt extdata} directory of the SPIA package, and therefore make it available to the function {\tt spia}.

The user has the ability to generate his own gene/protein relation data and put it in a list format as the 
one shown in the hsaSPIA.RData file. 
In this file, each pathway data is included in a list:
<<eval=TRUE,echo=TRUE,results=verbatim>>=
load(file=paste(system.file("extdata/hsaSPIA.RData",package="SPIA")))
names(path.info[["05210"]])
path.info[["05210"]][["activation"]][25:35,30:40]
@
In the matrix above, only 0 and 1 values are allowed. 1 means the gene/protein given by the column has a relation of type "activation" with the gene/protein given by the row of the matrix.

Using other R packages such as {\tt graph} and {\tt Rgraphviz} one can visualize the richness of gene/protein relations of each type in each 
pathway. 
Firstly we load the required packages and create a function that can be used to plot as a graph each type of relation of 
any pathway, as used by SPIA. 

<<eval=TRUE,echo=TRUE,results=verbatim>>=
library(graph)
library(Rgraphviz)

plotG<-function(B){
 nnms<-NULL;colls<-NULL
 mynodes<-colnames(B)
 L<-list();
 n<-dim(B)[1]
 for (i in 1:n){
 L[i]<-list(edges=rownames(B)[abs(B[,i])>0])
 if(sum(B[,i]!=0)>0){
 nnms<-c(nnms,paste(colnames(B)[i],rownames(B)[B[,i]!=0],sep="~"))
 }
 }
 names(L)<-rownames(B)
 g<-new("graphNEL",nodes=mynodes,edgeL=L,edgemode="directed")
 plot(g)
}

@


We plot then the "activation" relations in the ErbB signaling pathway, based on the {\tt hsaSPIA} data.

\begin{figure}[htbp]
\label{fig:Graph}
  \begin{center}
<<fig=TRUE>>=
plotG(path.info[["04012"]][["activation"]])
@
\caption{Display of the "activation" relations in the ErbB signaling pathway, based on the hsaSPIA data.}
\end{center}
\end{figure}


For more details on how to use the main function in this package use "?spia".

\bibliography{SPIA} 
\end{document}