%\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}
\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}}
\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 Markowetz \emph{et al.} (2005) and Markowetz (2006) 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} (Boutros \emph{et al.}, 2002). 
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{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:
<<echo=FALSE,results=hide>>=
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 et al (2005)) or the full
marginal likelihood depending on hyperparameters (details in
Markowetz, 2006). Additionally we show how to employ an edge-wise
inference heuristic, which can also be applied to cases where
exhaustive search over model space is infeasible (i.e. when we have
more than 4 or 5 perturbed genes).
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")
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{mLL}. 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 (Markowetz \emph{et al.}, 2005)
the PhD thesis (Markowetz, 2006) 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")
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_positions}{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.
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).

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

\myincfig{.3}{
<<label=graph3,fig=TRUE,echo=FALSE,results=hide>>=
col<-c("yellow","yellow","green","blue")
names(col) = nodes(result3$graph)
plot(result3, 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.
In the package \texttt{nem} we do it by
<<eval=FALSE>>=
result4 <- nem(res.disc$dat,para=res.disc$para,inference="triples")
result4
@

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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>>=
result3.scc <- SCCgraph(result3$graph,name=TRUE)
plot(result3.scc$graph)
@

\myincfig{.3}{
<<label=scc,fig=TRUE,echo=FALSE,results=hide>>=
col2<-c("yellow","blue","green")
names(col2) = nodes(result3.scc$graph)
plot(result3.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\begin{thebibliography}{14}

\bibitem[Boutros {\it et~al.} (2002)]{boutros02sequential}
Boutros M, Agaisse H, Perrimon N.
Sequential activation of signaling pathways during innate immune responses in \emph{Drosophila}.
{\em Developmental Cell}, 3(5):711--722, 2002.

\bibitem[Huber {\it et~al.} (2002)]{huber02variance}
Huber W, Heydebreck A, S{\"u}ltmann H, Poustka A, Vingron M.
Variance Stabilization Applied to Microarray Data Calibration and to the Quantification of Differential Expression.
{\em Bioinformatics}, 18:S96--S104, 2002.

\bibitem[Irizarry {\it et~al.} (2003)]{irizarry03summaries}
Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP.
Summaries of Affymetrix GeneChip probe level data.
{\em Nucleic Acids Res.} 31(4):e15, 2003.

\bibitem[Markowetz {\it et~al.} (2005)]{markowetz05non-transcriptional}
Markowetz F, Bloch J, Spang R.
Non-transcriptional pathway features reconstructed from secondary effects of RNA interference
{\em Bioinformatics}, 2005

\bibitem[Markowetz (2006)]{markowetz05thesis}
Markowetz F.
Probabilistic Models for Gene Silencing Data.
{\em PhD thesis}, Free University Berlin, 2005

\end{thebibliography}

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