%\VignetteIndexEntry{Nested Effects Models - An example in Drosophila immune response}
%\VignetteDepends{}
%\VignetteKeywords{Pathways}
%\VignettePackage{nem}


\documentclass[11pt,a4paper]{article}

%\usepackage[round]{natbib}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage[latin1]{inputenc}
\newcommand{\gene}[1]{\emph{#1}}

\setlength{\parskip}{1.5ex}
\setlength{\parindent}{0cm}

% NEW COMMANDS
% ------------------------------------------------
\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\Rpackage}[1]{\textit{#1}}
\newcommand{\Rclass}[1]{\textit{#1}}
\newcommand{\Rfunction}[1]{{\small\texttt{#1}}}

\newcommand{\myincfig}[4]{
  \setkeys{Gin}{width=#1\textwidth}
  \begin{figure}[htbp]
    \begin{center}
      #2
      \caption{\label{#3}#4}
    \end{center}
  \end{figure}
  \setkeys{Gin}{width=.8\textwidth}
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

<<no.nonsense,echo=FALSE,results=hide>>=
rm(list=ls())
@

\title{--- Nested Effects Models --- \\ An example in \emph{Drosophila} immune response}
\author{Florian Markowetz\footnote{Lewis-Sigler Institute for Integrative Genomics, Pinceton, NJ 08544, USA. eMail: florian@genomics.princeton.edu; URL: http://genomics.princeton.edu/$\sim$florian} \and Holger Fr\"ohlich\footnote{German Cancer Research Center, Im Neuenheimer Feld 580, 69120 Heidelberg, Germany. eMail: h.froehlich@dkfz-heidelberg.de}}
\date{\today}
\maketitle

\begin{abstract}
Cellular signaling pathways, which are not modulated on a transcriptional level,
cannot be directly deduced from expression profiling experiments. The situation
changes, when external interventions like RNA interference or gene knock-outs
come into play. 

In \cite{Markowetz2005Inference}, \cite{Markowetz2006Thesis}, \cite{Markowetz2007} and \cite{Frohlich2007GCBRNAi} we introduced an algorithm to infer non-transcriptional pathway features based on differential
gene expression in silencing assays. 
The method is implemented in the Bioconductor package \texttt{nem}.
Here we demonstrate its practical use in
the context of an RNAi data set investigating the response to microbial
challenge in Drosophila melanogaster. 

We show in detail how the data is pre-processed and discretized,
how the pathway can be reconstructed by different approaches, and how the final result can be post-processed 
to increase interpretability.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Drosophila RNAi data}

We applied our method to data from a study on innate
immune response in {\em Drosophila} \cite{Boutros2002Data}. 
Selectively removing
signaling components blocked induction of all, or
only parts, of the transcriptional response to LPS.

%-----------------------------------------------------------
\paragraph{Dataset summary}
The dataset consists of 16 Affymetrix-microarrays: 4 replicates of control
experiments without LPS and without RNAi (negative controls), 4 replicates of
expression profiling after stimulation with LPS but without RNAi (positive
controls), and 2 replicates each of expression profiling after applying LPS and
silencing one of the four candidate genes \gene{tak}, \gene{key}, \gene{rel},
and \gene{mkk4/hep}.


%-----------------------------------------------------------
\paragraph{\label{preprocess}Preprocessing and E-gene selection}
For preprocessing, we perform normalization on probe level using a variance
stabilizing transformation (Huber \emph{et al.}, 2002), and probe set summarization
using a median polish fit of an additive model (Irizarry \emph{et al.}, 2003).
The result is included as a dataset in the package \texttt{nem}.

<<results=hide>>=
library(nem)
data("BoutrosRNAi2002")
@

The function \texttt{nem.discretize} implements the following two preprocessing steps:
First, we select the genes as effect reporters (E-genes), which are more then two-fold upregulated by LPS
treatment.
Next, we transform the continuous expression data to binary values. We set an
E-genes state in an RNAi experiment to \texttt{1} if its expression value is
sufficiently far from the mean of the positive controls, \emph{i.e.} if the
intervention interrupted the information flow. If the E-genes expression is
close to the mean of positive controls, we set its state to \texttt{0}.

Let $C_{ik}$ be the continuous expression level of $E_i$ in experiment $k$. Let
$\mu^+_i$ be the mean of positive controls for $E_i$, and $\mu^-_i$ the mean of
negative controls. To derive binary data $E_{ik}$, we defined individual
cutoffs for every gene $E_i$ by:

\begin{equation}
E_{ik} =
\begin{cases}
\ 1 & \text{if }C_{ik} < \kappa\cdot\mu^+_i + (1-\kappa)\cdot\mu^-_i ,\\
\ 0 & \text{else}.
\end{cases}
\label{discretization}
\end{equation}
<<>>=
res.disc <- nem.discretize(BoutrosRNAiExpression,neg.control=1:4,pos.control=5:8,cutoff=.7)
@

<<echo=FALSE, results=hide>>=
disc <- cbind(res.disc$neg,res.disc$pos,res.disc$dat)
e.2fold <- BoutrosRNAiExpression[res.disc$sel,]

#--- hierarchisch clustern
hc    <- hclust(as.dist(hamming.distance(disc[,9:16])))
e.2fold <- e.2fold[hc$order, ]
disc  <- disc [hc$order, ]
@

<<label=data_cont,fig=TRUE,echo=FALSE,results=hide,include=FALSE,width=5,height=13>>=
#--- CONTINUOUS DATA
#pdf("data_cont.pdf",width=5,height=13)
par(las=2,mgp=c(5.5,1,0),mar=c(6.7,7,4,1),cex.lab=1.7,cex.main=2)
image(x   = 1:ncol(e.2fold),
y   = 1:nrow(e.2fold),
z   = scale(t(e.2fold)),
main= "Original data",
xlab= "Experiments",
xaxt= "n",
ylab= "E-genes",
yaxt= "n",
col = gray(seq(0,.95,length=10))
)
abline(v=c(4,8,10,12,14)+.5)
axis(1,1:ncol(e.2fold),colnames(e.2fold))
axis(2,1:nrow(e.2fold),rownames(e.2fold))
#dev.off()
@

<<label=data_disc,fig=TRUE,echo=FALSE,results=hide,include=FALSE,width=5,height=13>>=
#--- DISCRETE DATA
#pdf("data_disc.pdf",width=5,height=13)
par(las=2,mgp=c(5.5,1,0),mar=c(6.7,7,4,1),cex.lab=1.7,cex.main=2)
image(x   = 1:ncol(disc),
      z   = t(disc),
      main= "Discretized data",
      xlab= "Experiments",
      xaxt= "n",
      ylab= "",
      yaxt= "n",
      col = gray(seq(.95,0,length=10))
      )
abline(v=c(4,8,10,12,14)+.5)
axis(1,1:ncol(e.2fold),colnames(e.2fold))
#dev.off()
@
\begin{figure}[t]
\begin{center}
\includegraphics[width=.4\textwidth]{nem-data_cont}
%\hspace{-1cm}
\includegraphics[width=.4\textwidth]{nem-data_disc}
\end{center}
\caption{\label{data}Continuous and discretized data}
\end{figure}

%-----------------------------------------------------------
\paragraph{Estimating error probabilities}
From the positive and negative controls, we can  estimate the error
probabilities $\alpha$ and $\beta$.
The type I error $\alpha$ is the number of positive controls discretized to state \texttt{1},
and the type II error $\beta$ is the number of negative controls in state \texttt{0}.
To guard against unrealistically low estimates we add pseudo counts.
The error estimates are included into the discretization results:
<<>>=
res.disc$para
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage\section{Applying Nested Effects Models}

Which model explains the data best?  With only four S-genes, we can exhaustively enumerate all \
pathway models and search the whole space for the best-fitting model.
To score these models use either the marginal likelihood depending on $\alpha$
and $\beta$ (details found in Markowetz \emph{et al}. (2005)) or the full
marginal likelihood depending on hyperparameters (details in
Markowetz, 2006). 

In cases, where
exhaustive search over model space is infeasible (i.e. when we have
more than 4 or 5 perturbed genes) several heuristics have been developed and integrated into the \emph{nem} package:
\begin{itemize}
\item edge-wise and triplets learning \cite{Markowetz2007}
\item greedy hillclimbing
\item module networks \cite{Frohlich2007GCBRNAi}
\end{itemize}
An interface to all inference techniques is provided by the function \texttt{nem()}.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Exhaustive search by marginal likelihood}

Scoring models by marginal log-likelihood is implemented in function
\texttt{score}.  Input consists of models and data, the type of the score
(\texttt{"mLL"} or \texttt{"FULLmLL"}), the corresponding paramters
(\texttt{para}) or hyperparameters (\texttt{hyperpara}) and a prior for E-gene
positions (\texttt{P}). 

<<>>=
result <- nem(res.disc$dat,type="mLL",para=res.disc$para,inference="search",verbose=FALSE)
result
@
The output is the highest scoring model (\texttt{result\$graph}), a vector of scores (\texttt{result\$mLL}) and a list of
E-gene position posteriors (\texttt{result\$pos}), and a MAP estimate of E-gene
positions (\texttt{result\$mappos}). 
We can plot the results using the commands:
<<results=hide>>=
plot(result,what="graph")
plot(result,what="mLL")
plot(result,what="pos")
@

%--- plotting the graph
<<echo=FALSE,results=hide>>=
Sgenes <- unique(colnames(res.disc$dat))
models <- enumerate.models(Sgenes)
best5 <- -sort(-unique(result$mLL))[1:5]
col<-c("yellow","yellow","green","blue")
names(col) = Sgenes
for (i in 1:5) {
   graph <- as(models[[which(result$mLL == best5[i])[1]]]-diag(4),"graphNEL")
   pdf(file=paste("topo",i,".pdf",sep=""))
   par(cex.main=5)
   plot(graph,
        nodeAttrs=list(fillcolor=col),
        main=paste("-",i, "-"))        
   dev.off()
   }
@


\begin{figure}[ht]
\begin{center}
\framebox{\includegraphics[width=.18\textwidth]{topo1}}
\hfill
\includegraphics[width=.18\textwidth]{topo2}
\includegraphics[width=.18\textwidth]{topo3}
\includegraphics[width=.18\textwidth]{topo4}
\includegraphics[width=.18\textwidth]{topo5}
\end{center}
\caption{\label{topologies}The five silencing schemes getting high scores in Fig.~\ref{scores1}. It takes a second to see it, but Nr.2 to 5 are not that different from Nr.1. The main feature, ie. the branching downstream of \gene{tak} is conserved in all of them.}
\end{figure}


%--- plotting the scores
\myincfig{.4}{
<<label=scores1,fig=TRUE,echo=FALSE,results=hide>>=
plot(result,what="mLL")
@
}
{scores1}{The best 30 scores}


%--- plotting the posterior positions
\myincfig{.4}{
<<label=pos1,fig=TRUE,echo=FALSE,results=hide,width=5,height=13>>=
plot(result,what="pos")
@
}
{Egene_positions}{Posterior distributions of E-gene positions given
the highest scoring silencing scheme (Nr. 1 in Fig.~\ref{topologies}). The MAP estimate
corresponds to the row-wise maximum.}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Exhaustive search Full marginal likelihood}
Additionally to what we did in the paper \cite{Markowetz2005Inference}
the PhD thesis \cite{Markowetz2006Thesis} contains equations for a ``full marginal
likelihood'' in which error probabilities $\alpha$ and $\beta$ are integrated
out. This section shows that using this score we learn the same pathways as
before.

<<>>=
result2 <- nem(res.disc$dat,type="FULLmLL",hyperpara=c(1,9,9,1),inference="search",verbose=FALSE)
result2
@



%--- plotting the graphs
<<echo=FALSE,results=hide>>=
best5 <- -sort(-unique(result2$mLL))[1:5]
for (i in 1:5) {
   graph <- as(models[[which(result2$mLL == best5[i])[1]]]-diag(4),"graphNEL")
   pdf(file=paste("topo2",i,".pdf",sep=""))
   par(cex.main=5)
   plot(graph,
        nodeAttrs=list(fillcolor=col),
        main=paste("-",i, "-"))
   dev.off()
   }
@


\begin{figure}[ht]
\begin{center}
\framebox{\includegraphics[width=.18\textwidth]{topo21}}
\hfill
\includegraphics[width=.18\textwidth]{topo22}
\includegraphics[width=.18\textwidth]{topo23}
\includegraphics[width=.18\textwidth]{topo24}
\includegraphics[width=.18\textwidth]{topo25}
\end{center}
\caption{\label{topologies2}Same topologies as before.}
\end{figure}


%--- plotting the scores
\myincfig{.4}{
<<label=scores2,fig=TRUE,echo=FALSE,results=hide>>=
plot(result2,what="mLL")
@
}
{scores2}{The best 30 scores by full marginal likelihood}

%--- plotting the posterior positions
\myincfig{.4}{
<<label=pos2,fig=TRUE,echo=FALSE,results=hide,width=5,height=13>>=
plot(result2,what="pos")
@
}
{Egene_positions2}{Posterior distributions of E-gene positions given
the highest scoring silencing scheme (Nr. 1 in Fig.~\ref{topologies2}). The MAP estimate
corresponds to the row-wise maximum.}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Edge-wise learning}

Instead of scoring whole pathways, we can learn the model edge by edge \cite{Markowetz2007}.
For each pair of genes $A$ and $B$ we infer the best of four possible models:
$A \cdot \cdot B$ (unconnected),
$A \rightarrow B$ (effects of A are superset of effects of B),
$A \leftarrow B$ (subset), and
$A \leftrightarrow B$ (undistinguishable effects).

<<>>=
resultPairs <- nem(res.disc$dat,para=res.disc$para,inference="pairwise",verbose=FALSE)
resultPairs
@

\myincfig{.3}{
<<label=graph3,fig=TRUE,echo=FALSE>>=
col<-c("yellow","yellow","green","blue")
names(col) = nodes(resultPairs$graph)
plot(resultPairs, nodeAttrs=list(fillcolor=col))
@
}{graph3}{Result of edge-wise learning. Compare this to the result from global search. It looks exactely the same.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Inference from triples}

Edge-wise learning assumes independence of edges. But this is not true in transitively closed graphs, where a direct edge must exist whenever there is a longer path between two nodes.
Natural extension of edge-wise learning is inference from triples of nodes \cite{Markowetz2007}.
In the package \texttt{nem} we do it by
<<>>=
resultTriples <- nem(res.disc$dat,para=res.disc$para,inference="triples",verbose=FALSE)
resultTriples
@

\myincfig{.3}{
<<label=graph4,fig=TRUE,echo=FALSE>>=
col<-c("yellow","yellow","green","blue")
names(col) = nodes(resultTriples$graph)
plot(resultTriples$graph, nodeAttrs=list(fillcolor=col))
@
}{graph4}{Result of triple learning. Compare this to the result from global search and pairwise learning}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Inference with greedy hillclimbing}

Greedy hillclimbing is a general search and optimization strategy known from the literature \cite{RusNor95}. Given an initial network hypothesis (usually an empty graph) we want to arrive at a local maximum of the likelihood function by successively adding that edge, which adds the most benefit to the current network's likelihood. This procedure is continued until no improving edge can be found any more.

<<>>=
resultGreedy <- nem(res.disc$dat,para=res.disc$para,inference="nem.greedy",verbose=FALSE)
resultGreedy
@

\myincfig{.3}{
<<label=graph44,fig=TRUE,echo=FALSE,results=hide>>=
col<-c("yellow","yellow","green","blue")
names(col) = nodes(resultGreedy$graph)
plot(resultGreedy$graph, nodeAttrs=list(fillcolor=col))
@
}{graph44}{Result of greedy hillclimbing. It is exactly the same as for the exhaustive search.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Inference with module networks}

Rather than looking for a complete network hypothesis at once the
idea of the module network is to build up a graph from smaller subgraphs,
called \emph{modules} in the following \cite{Frohlich2007GCBRNAi}. The module network is thus
a divide and conquer approach: We first perform a hierarchical clustering
of the preprocessed expression profiles of all S-genes, e.g. via average linkage. The
idea is that S-genes with a similar E-gene response profile (here:
with regard to the Pearson correlation similarity) should be close in the signaling
path. We now successively move down the cluster tree hierarchy until we find a cluster with 
only 4 S-genes at most. Figure \ref{fig:MNIdea} illustrates the idea with an assumed network of 10 S-genes.
At the beginning we find $S_8$ as a cluster singleton. Then by successively moving down the hierarchy we identify 
clusters $S_6, S_7, S_1, S_10$, $S_3, S_2, S_5$ and $S_4, S_9$. All these clusters (modules) contain 4 S-genes at most and can 
thus be estimated by taking the highest scoring of all possible network hypotheses. 

\begin{figure}
\includegraphics[scale=0.5]{ModuleNetworks1}

\caption{\label{fig:MNIdea} Basic idea of module networks: By successively moving down the cluster hierarchy we identify the clusters (modules) of S-genes, which are marked in red. They contain 4 S-genes at most and can be estimated via exhaustive search.}
\end{figure}

Once all modules have been estimated their connections are constructed. This is done in a constraint greedy hillcimbing fashion: We successively add that edge between any pair of S-genes being contained in different modules, which increases the
likelihood of the complete network most. This procedure is continued until no improvement can be gained any more, i.e. we have reached a local maximum of the likelihood function. 

In the package \texttt{nem} we call the module network by
<<>>=
resultMN <- nem(res.disc$dat,para=res.disc$para,inference="ModuleNetwork",verbose=FALSE)
resultMN # will do exactly the same as exhaustive search in this case
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Incorporating prior Assumptions}

\subsubsection{Regularization}

The \texttt{nem} package allows to specify a prior on the network structure itself. 
This can be thought of biasing the score
of possible network hypotheses towards prior knowledge. It is crucial
to understand that in principle in any inference scheme there exist
two competing goals: Belief in prior assumptions / prior knowledge
versus belief in data. Only trusting the data itself may lead to overfitting,
whereas only trusting in prior assumptions does not give any new information
and prevents learning. Therefore, we need a trade-off between both
goals via a regularization constant $\lambda > 0$, which specifies the belief in our prior assumptions. In the simplest case our 
assumption could be that the true network structure is sparse, i.e. there are only very few edges. More complex priors could involve separate prior edge probabilities.

<<>>=
resultRegularization <- nem(res.disc$dat,para=res.disc$para,Pm=matrix(0,ncol=4,nrow=4), lambda=10, inference="search", verbose=FALSE)
resultRegularization
@

\myincfig{.3}{
<<label=graph6,fig=TRUE,echo=FALSE,results=hide>>=
col<-c("yellow","yellow","green","blue")
names(col) = nodes(resultRegularization$graph)
plot(resultRegularization, nodeAttrs=list(fillcolor=col))
@
}{graph6}{Result of module network learning with regularization towards sparse graph structures ($\lambda=10$).}


In practice we would like to choose a $\lambda$ in an automated fashion.
This leads to an instance of the classical \emph{model selection}
problem (e.g. \cite{HastieTibshiraniBook2001}) in statistical learning.
One way of dealing with it is to tune $\lambda$ such that the \emph{Akaike
information criterion} (AIC)\begin{equation}
AIC(\lambda,\Phi_{opt})=-2\log P(D|\Phi_{opt})+2d(\lambda,\Phi_{opt})\end{equation}
becomes minimal \cite{HastieTibshiraniBook2001}. Here $d(\lambda,\Phi_{opt})$
denotes the number of parameters (i.e. the number of edges) in the highest scoring 
network structure $\Phi_{opt}$.

<<>>=
resultModsel <- nemModelSelection(c(0.1,1,10),res.disc$dat,para=res.disc$para,Pm=matrix(0,ncol=4,nrow=4),inference="search", verbose=FALSE)
@


\myincfig{.3}{
<<label=graph7,fig=TRUE,echo=FALSE,results=hide>>=
col<-c("yellow","yellow","green","blue")
names(col) = nodes(resultModsel$graph)
plot(resultModsel, nodeAttrs=list(fillcolor=col))
@
}{graph7}{Result of module network learning with regularization towards sparse graph structures and automatic model selection.}


\subsubsection{Bayesian Model Selection}

Searching for an optimal regularization constant relates to a frequentistic point of view to incorporate prior knowledge. Instead, from a Bayesian perspective one should define a prior on the regularization parameter and integrate it out. Here, this is done by assuming an inverse gamma distribution prior on $\nu=\frac{1}{2\lambda}$ with hyperparameters $1,0.5$, which leads to a simple closed form of the full prior \cite{Frohlich2007GCBRNAi}. An advantage of the Bayesian perspective is that no explicit model selection step is needed.

<<>>=
resultBayes <- nem(res.disc$dat,para=res.disc$para,Pm=matrix(0,ncol=4,nrow=4),inference="search", verbose=FALSE)
resultBayes
@


\myincfig{.3}{
<<label=graph77,fig=TRUE,echo=FALSE,results=hide>>=
col<-c("yellow","yellow","green","blue")
names(col) = nodes(resultBayes$graph)
plot(resultBayes, nodeAttrs=list(fillcolor=col))
@
}{graph77}{Result of module network learning with a Bayesian network prior favoring sparse graphs.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Omitting the Data Discretization Step}

In general performing a data discretization on the expression profiles as described in Sec. \ref{preprocess} can be critical. An alternative is given by taking the 
raw $p$-value profiles obtained from testing for differential gene expression. In this situation we assume the individual $p$-values in the data matrix to be drawn from a mixture of a uniform, a Beta$(1,\beta)$ and a Beta$(\alpha,1)$ distribution. The parameters of the distribution are fitted via an EM algorithm \cite{Frohlich2007GCBRNAi}. The \texttt{nem} package supports such a data format using the option \texttt{type = "CONTmLLDens"} in the call of the function \texttt{nem}. Moreover there is a function \texttt{getDensityMatrix}, which conveniently does all the fitting of the $p$-value profiles and produces diagnostic plots into a user specified directory.

A second possibility to omit the data discretization step is to calculate the effect probability for each gene based on given the empirical distributions of the controls. 

<<eval=FALSE>>=
densities = getDensityMatrix(myPValueMatrix,dirname="DiagnosticPlots")
nem(densities[res.disc$sel,],type="CONTmLLDens",inference="search")

preprocessed <- nem.cont.preprocess(BoutrosRNAiExpression,neg.control=1:4,pos.control=5:8)
nem(preprocessed$prob.influenced,type="CONTmLL",inference="search")
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Visualization}
<<results=hide,fig=FALSE>>=
plot.effects(res.disc$dat,result)
@

\myincfig{.4}{
<<label=plot_effects,fig=TRUE,echo=FALSE,results=hide,width=5,height=13>>=
plot.effects(res.disc$dat,result)
@
}{plot_effects}{plotting data according to inferred hierarchy}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Post-processing of results}


\paragraph{Combining strongly connected components}
First, we identify all nodes/genes which are not distinguishable given the data.
This amounts to finding the strongly connected components in the graph. Relish
and Key are now combined into one node.


<<fig=FALSE>>=
result.scc <- SCCgraph(result$graph,name=TRUE)
plot(result.scc$graph)
@

\myincfig{.3}{
<<label=scc,fig=TRUE,echo=FALSE,results=hide>>=
# col2<-c("yellow","blue","green")
# names(col2) = nodes(result.scc$graph)
plot(result.scc$graph)#,nodeAttrs=list(fillcolor=col2))
@
}{scc}{The undistinguishable profiles of \gene{key} and \gene{rel} are summarized into a single node.}


\paragraph{Transitive reduction}
Additionally, in bigger graphs \texttt{transitive.reduction()} helps to see the
structure of the network better. In this simple example there are no shortcuts to remove.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\bibliographystyle{abbrv}
\bibliography{references}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Session Information}

The version number of R and packages loaded for generating the vignette were:


<<echo=FALSE,results=tex>>=
toLatex(sessionInfo())
@


\end{document}

%
%   end of file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
