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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{E:/biocbld/bbs-2.2-bioc/R/share/texmf/Sweave}
\begin{document}


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

\begin{Schunk}
\begin{Sinput}
> library(nem)
> data("BoutrosRNAi2002")
\end{Sinput}
\end{Schunk}

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}
\begin{Schunk}
\begin{Sinput}
> res.disc <- nem.discretize(BoutrosRNAiExpression, neg.control = 1:4, 
+     pos.control = 5:8, cutoff = 0.7)
\end{Sinput}
\begin{Soutput}
discretizing with respect to POS and NEG controls
\end{Soutput}
\end{Schunk}



\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:
\begin{Schunk}
\begin{Sinput}
> res.disc$para
\end{Sinput}
\begin{Soutput}
   a    b 
0.19 0.07 
\end{Soutput}
\end{Schunk}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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}). 

\begin{Schunk}
\begin{Sinput}
> result <- nem(res.disc$dat, type = "mLL", para = res.disc$para, 
+     inference = "search", verbose = FALSE)
> result
\end{Sinput}
\begin{Soutput}
Object of class 'score' generated by 'score()'

$graph:  phenotypic hierarchy on genes
$mLL:    a vector of length 355 
$pos:    a list   of length 355 
$mappos: a vector of length 355 
$lambda: 0 
\end{Soutput}
\end{Schunk}
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:
\begin{Schunk}
\begin{Sinput}
> plot(result, what = "graph")
> plot(result, what = "mLL")
> plot(result, what = "pos")
\end{Sinput}
\end{Schunk}

%--- plotting the graph


\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}{
\includegraphics{nem-scores1}
}
{scores1}{The best 30 scores}


%--- plotting the posterior positions
\myincfig{.4}{
\includegraphics{nem-pos1}
}
{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.

\begin{Schunk}
\begin{Sinput}
> result2 <- nem(res.disc$dat, type = "FULLmLL", hyperpara = c(1, 
+     9, 9, 1), inference = "search", verbose = FALSE)
> result2
\end{Sinput}
\begin{Soutput}
Object of class 'score' generated by 'score()'

$graph:  phenotypic hierarchy on genes
$mLL:    a vector of length 355 
$pos:    a list   of length 355 
$mappos: a vector of length 355 
$lambda: 0 
\end{Soutput}
\end{Schunk}



%--- plotting the graphs


\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}{
\includegraphics{nem-scores2}
}
{scores2}{The best 30 scores by full marginal likelihood}

%--- plotting the posterior positions
\myincfig{.4}{
\includegraphics{nem-pos2}
}
{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).

\begin{Schunk}
\begin{Sinput}
> resultPairs <- nem(res.disc$dat, para = res.disc$para, inference = "pairwise", 
+     verbose = FALSE)
> resultPairs
\end{Sinput}
\begin{Soutput}
Object of class 'pairwise' generated by 'pairwise.posterior()'

$graph:  phenotypic hierarchy on genes)
$scores: posterior distributions of local models

Summary of MAP estimates:
all  ..  -> <-> 
  6   2   3   1 
\end{Soutput}
\end{Schunk}

\myincfig{.3}{
\includegraphics{nem-graph3}
}{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
\begin{Schunk}
\begin{Sinput}
> resultTriples <- nem(res.disc$dat, para = res.disc$para, inference = "triples", 
+     verbose = FALSE)
\end{Sinput}
\begin{Soutput}
4 perturbed genes -> 4 triples to check (lambda =  0 )
\end{Soutput}
\begin{Sinput}
> resultTriples
\end{Sinput}
\begin{Soutput}
Object of class 'triples' generated by 'triples.posterior()'

$graph:  phenotypic hierarchy on genes
$avg:    matrix of edge frequencies in triple models
\end{Soutput}
\end{Schunk}

\myincfig{.3}{
\includegraphics{nem-graph4}
}{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.

\begin{Schunk}
\begin{Sinput}
> resultGreedy <- nem(res.disc$dat, para = res.disc$para, inference = "nem.greedy", 
+     verbose = FALSE)
\end{Sinput}
\begin{Soutput}
Greedy hillclimber for 4 S-genes (lambda = 0 )...
\end{Soutput}
\begin{Sinput}
> resultGreedy
\end{Sinput}
\begin{Soutput}
Object of class 'nem.greedy' generated by 'nem.greedy()'

$graph:  phenotypic hierarchy on genes)
\end{Soutput}
\end{Schunk}

\myincfig{.3}{
\includegraphics{nem-graph44}
}{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{Frohlich2007RNAiBMC, 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
\begin{Schunk}
\begin{Sinput}
> resultMN <- nem(res.disc$dat, para = res.disc$para, inference = "ModuleNetwork", 
+     verbose = FALSE)
\end{Sinput}
\begin{Soutput}
Estimating module network of 4 S-genes (lambda = 0 )...
\end{Soutput}
\begin{Sinput}
> resultMN
\end{Sinput}
\begin{Soutput}
Object of class 'ModuleNetwork' generated by 'moduleNetwork()'

$graph:  phenotypic hierarchy on genes)
\end{Soutput}
\end{Schunk}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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 (c.f. \cite{Frohlich2007RNAiBMC}).

\begin{Schunk}
\begin{Sinput}
> resultRegularization <- nem(res.disc$dat, para = res.disc$para, 
+     Pm = matrix(0, ncol = 4, nrow = 4), lambda = 10, inference = "search", 
+     verbose = FALSE)
> resultRegularization
\end{Sinput}
\begin{Soutput}
Object of class 'score' generated by 'score()'

$graph:  phenotypic hierarchy on genes
$mLL:    a vector of length 355 
$pos:    a list   of length 355 
$mappos: a vector of length 355 
$lambda: 10 
\end{Soutput}
\end{Schunk}

\myincfig{.3}{
\includegraphics{nem-graph6}
}{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 unknown edges) in the highest scoring 
network structure $\Phi_{opt}$.

\begin{Schunk}
\begin{Sinput}
> resultModsel <- nemModelSelection(c(0.01, 1, 100), res.disc$dat, 
+     para = res.disc$para, Pm = matrix(0, ncol = 4, nrow = 4), 
+     verbose = FALSE)
\end{Sinput}
\begin{Soutput}
Greedy hillclimber for 4 S-genes (lambda = 0.01 )...

Greedy hillclimber for 4 S-genes (lambda = 1 )...

Greedy hillclimber for 4 S-genes (lambda = 100 )...
\end{Soutput}
\end{Schunk}


\myincfig{.3}{
\includegraphics{nem-graph7}
}{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. Furthermore, there is evidence, that compared to the frequentistic method the Bayesian approach using the same amount of prior knowledge yields a higher increase of the reconstructed network's sensitivity 

\begin{Schunk}
\begin{Sinput}
> resultBayes <- nem(res.disc$dat, para = res.disc$para, Pm = matrix(0, 
+     ncol = 4, nrow = 4), verbose = FALSE)
\end{Sinput}
\begin{Soutput}
Greedy hillclimber for 4 S-genes (lambda = 0 )...
\end{Soutput}
\begin{Sinput}
> resultBayes
\end{Sinput}
\begin{Soutput}
Object of class 'nem.greedy' generated by 'nem.greedy()'

$graph:  phenotypic hierarchy on genes)
\end{Soutput}
\end{Schunk}


\myincfig{.3}{
\includegraphics{nem-graph77}
}{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, especially, if not both, positive and negative controls are available. 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 alternative to the $p$-value density, is to use the log-odds (B-values), which compare the likelihood of differential expression with the likelihood of non-differential expression of a gene. The \texttt{nem} package supports such a data format using the option \texttt{type = "CONTmLLRatio"} in the call of the function \texttt{nem}. 

A further 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. Note that for this purpose both, positive and negative controls should be available.

\begin{Schunk}
\begin{Sinput}
> logdensities = getDensityMatrix(myPValueMatrix, dirname = "DiagnosticPlots")
> nem(logdensities[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")
\end{Sinput}
\end{Schunk}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Selection of E-Genes}

There exists two mechanisms to select E-genes implemented in \texttt{nem}, which can be used complementary to each other: 
\begin{itemize}
\item A-priori filtering of E-genes, which show a pattern of differential expression, which can be expected to be non-random
\item Selection of most relevant E-genes for a given network hypothesis. 
\end{itemize}
In case of the log-odds inference scheme (\texttt{type = "CONTmLLRatio"}) selection of most relevant E-genes works as follows: A virtual "null" S-gene is used, to which E-genes are assigned, that are irrelevant. The prior probability of assigning an E-gene to the "null" S-genes is $\frac{\delta}{n}$, where $\delta>0$ and $n$ is the number of S-genes. The parameter $\delta$ is tuned via the AIC model selection criterion.

In all other cases selection of relevant E-genes is performed via an alternating optimization scheme: First, the network with the highest likelihood for the given set of E-genes is determined. Then, the E-genes with the highest marginal likelihood under that network hypothesis are calculated. In case of the $p$-value profile inference scheme (\texttt{type = "CONTmLLDens"}) these are simply all E-genes with a positive contribution to the log marginal likelihood. In all other cases by default the top $10n$ E-genes are selected. With the new set of E-genes a new network is inferred and the whole procedure is repeated until convergence.

\begin{Schunk}
\begin{Sinput}
> mydat = filterEGenes(Porig, logdensities)
> resAuto = nem(mydat, type = "CONTmLLDens", selEGenes = TRUE)
\end{Sinput}
\end{Schunk}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Statistical Stability and Significance}

An important step after network inference is some kind of validation. In abscence of the true network this step is rather difficult. We have implemented several possible tests, which may also be used in addition to each other:
\begin{itemize}
\item Statistical significance: Possibly the least demand we have for our inferred network is, that it should be significantly superior to a random one. We sample $N$ (default 1000) random networks from a null distribution and compare their marginal likelihood to that of the estimated network. The fraction of how often a random network is better than the inferred one yields an exact p-value.
\item Bootstrapping: We wish our network to be stable against small changes in our set of E-genes. Therefore, we use bootstrapping: We sample $m$ E-genes with replacement for $N$ times (default 1000) and run the network inference on each bootstrap sample. At the end we count for each edge how often it was inferred in all $N$ bootstrap runs. This yields a probability for each edge. We then may only consider edges with a probability above some threshold (e.g. 50%).
\item Jackknife: We also wish our network to be stable against removal of S-genes. Therefore, we use the jackknife: Each S-gene is left out once and the network inference run on the other S-genes. At the end we count for each edge, how often it was inferred. This yields a jackknife probability for each edge, which then may be used to filter edges.
\end{itemize}

\begin{Schunk}
\begin{Sinput}
> significance = nem.calcSignificance(res, disc$dat, as(resultGreedy$graph, 
+     "matrix"), resultGreedy$mLL, para = res.disc$para)
> bootres = nem.bootstrap(res.disc$dat, para = res.disc$para)
> jackres = nem.jackknife(res.disc$dat, para = res.disc$para)
\end{Sinput}
\end{Schunk}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\section{Visualization}
\begin{Schunk}
\begin{Sinput}
> plot.effects(res.disc$dat, result)
\end{Sinput}
\end{Schunk}

\myincfig{.4}{
\includegraphics{nem-plot_effects}
}{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.


\begin{Schunk}
\begin{Sinput}
> result.scc <- SCCgraph(result$graph, name = TRUE)
> plot(result.scc$graph)
\end{Sinput}
\end{Schunk}

\myincfig{.3}{
\includegraphics{nem-scc}
}{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:


\begin{itemize}
  \item R version 2.7.0 (2008-04-22), \verb|i386-pc-mingw32|
  \item Locale: \verb|LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252|
  \item Base packages: base, datasets, graphics, grDevices, grid,
    methods, stats, tools, utils
  \item Other packages: bootstrap~1.0-21, class~7.2-42,
    cluster~1.11.10, e1071~1.5-18, graph~1.18.1, nem~2.4.0,
    RBGL~1.16.0, Rgraphviz~1.18.0, time~1.0
\end{itemize}

\end{document}

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