%\VignetteIndexEntry{Dynamic Deterministic Effects Propagation Networks - exemplary workflow}
%\VignetteDepends{}
%\VignetteKeywords{Pathways}
%\VignettePackage{ddepn}


\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{Dynamic Deterministic Effects Propagation Networks (DDEPN) - exemplary workflow}
\author{ Christian Bender \footnote{German Cancer Research Center, Im Neuenheimer Feld 580, 69120 Heidelberg, Germany. eMail: c.bender@dkfz-heidelberg.de} }
\date{\today}
\maketitle


%\begin{abstract}
%Network modelling in systems biology has become an important tool to
%study molecular interactions, especially in the medical field like
%cancer research. The understanding of the interplay of proteins in
%cellular signalling is the basis for the development of novel drugs
%and therapies. Here, we set up a new method for the reconstruction of
%signalling networks from time course protein data after external
%perturbation. We show how to use protein expression and
%phosphorylation data measured on Reverse Phase Protein Arrays to infer
%a signalling network among proteins of the ERBB signalling cascade in
%a human breast cancer cell line. 
%
%Our method (\cite{Bender2010}) models the signalling dynamics by a boolean signal
%propagation mechanism that defines a sequence of state transitions for
%a given network structure. A likelihood score is proposed that
%describes the probability of our measurements given a particular state
%transition matrix. We identify the optimal sequence of state
%transitions via a Hidden Markov Model. Network structure search is
%performed by a genetic algorithm that optimises the overall likelihood
%of a population of candidate networks. We test our method on simulated
%networks and data and show its increased performance in comparison to
%another Dynamical Bayesian Network approach. The reconstruction of a
%network in our real data results in several known signalling chains
%from the ERBB network, showing the validity and usefulness of our
%approach.
%\end{abstract}
This document describes the package 'ddepn' implementing our proposed 
network reconstruction algorithm (\cite{Bender2010,Bender2011,Bender2011a}). Section \ref{sec:quickstart}
gives a quick introduction on how to use the package. Section \ref{sec:formatdata}
shows how to format the input arguments for the ddepn function. Section 
\ref{sec:prior} gives an overview on how to include prior knowledge on 
network structures and in section \ref{sec:usecases} code snippets are given 
for the possible types of calling ddepn. Finally, in section \ref{sec:application}
an example for inference on real data is given.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{QuickStart: using DDEPN for network inference on simulated data sets}
\label{sec:quickstart}
This section shows an exemplary workflow to reconstruct a signalling network from 
simulated data. Details on formatting the input data matrix as well as arguments
for the function calls can be found in subsequent sections. An analysis on real
data can be performed analogously and an example is given at the end of this
document, showing a simple analysis of longitudinal data measuring protein 
phosphorylation in HCC1954 breast cancer cells, generated on Reverse Phase
Protein Arrays.

%% chunk 1
<<label=Load_package,echo=FALSE, eval=FALSE>>=
library(ddepn)
@

\subsection{Simulating data}
In this section we show how to generate artificial networks and data. A reference
signalling network is simulated and used to sample measurements that incorporate
the network structure.

First, simulate a network with 6 nodes and 2 distinct input stimuli. 

<<SimulateNetwork, eval=FALSE>>=
set.seed(12345)
n <- 6
signet <- signalnetwork(n=n, nstim=2, cstim=0, prop.inh=0.2)
net <- signet$phi
stimuli <- signet$stimuli
weights <- signet$weights
@

Second get intensities for each protein that are based on the network structure generated above.

<<SimulateData, eval=FALSE>>=
#plotdetailed(net,weights=weights,stimuli=stimuli,fontsize=15)
## sample data
dataset <- makedata(net, stimuli, mu.bg=1200, sd.bg=400, mu.signal.a=2000, sd.signal.a=1000)
data <- dataset$datx
@


\subsection{Running the Genetic Algorithm (GA)}
Now run the genetic algorithm to reconstruct the network from the data generated above and 
compare it to the originally sampled network \emph{net}.

<<RunGA, results=hide, eval=FALSE>>=
ret <- ddepn(data, phiorig=net, inference="netga",
              maxiterations=150, p=50, q=0.3, m=0.8,
              usebics=TRUE)
@
After the reconstruction, the generated network can be viewed as follows:
<<PlotTheResults, eval=FALSE>>=
plotrepresult(ret)
@
\begin{figure}[htp]
\centering
\includegraphics[width=12cm,height=4cm]{netgaresult}
\caption{Result plot of the genetic algorithm example. Left: the original graph; Middle: the inferred graph; 
	     Right: sensitivity/specificity plot for comparing the original and inferred graphs}
\label{fig:resultga}
\end{figure}

The destination file of the output of the \emph{netga} function can be specified by argument
\emph{outfile}. Figure \ref{fig:scorenetgaoutput} at the end of this document shows the output
for the GA inference of figure \ref{fig:resultga}.


\subsection{Running Markov Chain Monte Carlo Sampling (inhibMCMC)}
An example run for MCMC sampling follows. Here, the package \emph{multicore} is needed, since two 
parallel and independent MCMC runs are performed. If \emph{multicore} is not available on the 
machine, just set \emph{multicores=FALSE} to perform sampling for a single chain.

<<MCMCsampling, echo=TRUE, results=hide, eval=FALSE>>=
# use the original network as prior probability matrix
B <- net
B[B==2] <- -1
# construct a prior matrix with uniform probabilities for each edge
if(require(multicore)) {
	ret <- ddepn(data, phiorig=net, inference="mcmc",
		maxiterations=15010, burnin=1000,
		usebics=FALSE, lambda=0.01, B=B,
		multicores=TRUE, cores=2, priortype="laplaceinhib")
}
@
After the sampling one can examine the sampling run:

<<Examine MCMC sampling, echo=TRUE, results=hide, eval=FALSE>>=
plotrepresult(ret$samplings[[1]])
@
\begin{figure}[htp]
\centering
\includegraphics[width=12cm,height=4cm]{netmcmcresult}
\caption{Result plot of MCMC sampling, analogous to figure \ref{fig:resultga}.}
\label{fig:resultmcmc}
\end{figure}

The returned list \emph{ret} contains two elements,
another list with name samplings (\emph{ret\$samplings}), which holds the different runs
for the MCMC. In case of \emph{multicores=FALSE}, only one run is performed and \emph{ret\$samplings}
holds only one element. Otherwise \emph{cores} runs are performed independently in parallel, 
and \emph{ret\$samplings} holds \emph{cores} elements. The second element in \emph{ret} with 
name \emph{ltraces} is a matrix and holds the score traces of all runs, where each column
corresponds to one trace. Output diagnostics can be produced using the R-package \emph{coda}. See
figure \ref{fig:mcmcoutput} for some example plots.

<<MCMC output, echo=TRUE, results=hide, eval=FALSE>>=
mcmc1 <- mcmc(data=ret$ltraces[-c(1:burnin,mi),1])
mcmc2 <- mcmc(data=ret$ltraces[-c(1:burnin,mi),2])
mcmcl <- mcmc.list(mcmc1,mcmc2)
plot(mcmcl)
gelman.plot(mcmcl)
@
\begin{figure}[h]
\centering
\includegraphics[width=10cm,height=8cm]{mcmcoutput}
\caption{Using package coda for some MCMC output analysis. A: Traces for 2 MCMC runs; B: 
Distributionplot of the sampling; C: gelman.plot of two MCMC samplings}
\label{fig:mcmcoutput}
\end{figure}

To check for convergence, Gelman's potential scale reduction can be computed with the function 
\emph{gelman\_diag} (a wrapper for the gelman.diag function in the coda package). Note that 
multiple chains have to be run for using this convergence statistic.

<<Gelmans potential scale reduction as convergence diagnostic, echo=TRUE, results=hide, eval=FALSE>>=
gelman_diag(ret)
@


Further, the output of the MCMC sampler can be analysed using the \emph{plot\_edgeconfidences} 
function. For each possible interaction, the number of activating and inhibiting edges 
(normalised to the total number of sampled edges) is shown as boxplot. If the counts 
for activations are significantly higher than the counts for inhibitions, an inclusion of 
this activation is suggested for generating a final network (and vice versa). If there is no 
difference in the number of activating and inhibiting edges, it is suggested to not include 
any interaction between the two nodes. Figure \ref{fig:mcmcedgeconf} shows the output of the 
\emph{plot\_edgeconfidences} function.  

<<Plotting edge confidences, echo=TRUE, results=hide, eval=FALSE>>=
plot_edgeconfidences(ret)
@

\begin{figure}[h]
\centering
\includegraphics[width=12cm,height=8cm]{mcmcedgeconf}
\caption{In each of the inhibMCMC chains, activation and inhibition edges are sampled. The boxes 
show the frequencies of activations (marked in red) and inhibitions (marked in blue) to the
total number of sampled edges (denoted as confidence, counted after the burn-in phase). Each 
panel represents the destination node of an edge. For example, in the upper left panel, all 
edges pointing to node X4 are shown. The box in the first column (denoted by X1-) represents 
the inhibition from X1 to X4, the box in the fourth column (denoted by X2+) the activation
from X2 to X4. The higher the difference between activation and inhibition, the more support
is found in the data for the respective edge.}
\label{fig:mcmcedgeconf}
\end{figure}

This information is used in the function \emph{create\_signetwork}, which creates a network 
from all significant interactions using statistical tests (wilcoxon rank-sum or t-test) 
for each interaction. Note that multiple parallel chains have to be run for valid testing 
(say at least 4 chains):
 
<<Plotting edge confidences, echo=TRUE, results=hide, eval=FALSE>>=
net <- create_signetwork_cv(ret, alpha=0.05, adj.method="none", plot=FALSE, type="wilcox",
		alternative="one.sided", paired=FALSE, ord=NULL, sel_policy="strict")
@


Also, the distributions of the data and the inferred active and passive state parameters can 
be visualised by the function \emph{plot\_profiles(ret)}. The plot is shown in 
figure \ref{fig:mcmcprofiles}.

<<Plotting the profiles and parameters, echo=TRUE, results=hide, eval=FALSE>>=
plot_profiles(ret,mfrow=c(3,4))
@

\begin{figure}[h]
\centering
\includegraphics[width=10cm,height=10cm]{mcmcprofiles}
\caption{Visualisation of the inferred model parameters for active and passive protein state. 
Shown are the data (as boxplots) and a spline fit for each protein over time (upper two rows), 
as well as the data and parameters as histogram plot (lower row).}
\label{fig:mcmcprofiles}
\end{figure}


As in \emph{netga}, the destination file of the output of the \emph{mcmc\_ddepn} function can
be specified by argument \emph{outfile}. Figure \ref{fig:resultmcmc} at the end of 
this document shows the output of one sampling run for the MCMC inference.


\subsection{Resuming the inference}

After an inference run is finished, it might be necessary or useful to resume the inference and 
continue, for example an MCMC run that did not converge. To do so, the function 
\emph{resume\_ddepn} can be used. For example, resume an inhibMCMC run and add another 100
iterations:

<<Resuming inference of inhibMCMC, echo=TRUE, results=hide, eval=FALSE>>=
## resuming the inference from an inhibMCMC run and add another 100 iterations
ret4 <- ddepn(dataset$datx,phiorig=phit, inference="mcmc", maxiterations=100,
		burnin=30, lambda=0.01, B=B, priortype="laplaceinhib", usebics=FALSE)
ret4 <- resume_ddepn(ret4,maxiterations=100)
@

Or resume a GA run and add another 30 iterations:
<<Resuming inference of netga, echo=TRUE, results=hide, eval=FALSE>>=
## resuming the inference from an netga run and add another 30 iterations
ret5 <- ddepn(dataset$datx,phiorig=phit, inference="netga", maxiterations=20, p=10, q=0.3, m=0.8, lambda=0.01, B=B, priortype="laplaceinhib", usebics=FALSE)
ret5 <- resume_ddepn(ret5,maxiterations=30)
@


\section{Notes on formatting constraints for the arguments of \emph{ddepn}}
\label{sec:formatdata}
There is only one neccessary argument to the function call of \emph{ddepn}: The data 
matrix \emph{dat}. Optionally,  a reference 
network \emph{phiorig} and seed networks \emph{phi} can be passed to \emph{ddepn}.
Each of these arguments is described briefly below.

\subsection{Input data matrix \emph{dat}}
The data matrix contains all measurements for the nodes in the rows (e.g. proteins or genes), 
and the experiments and time points in the columns. There are some special needs on how to name
the columns. We allow several treatments to be included in the data matrix. Examples for these
treatments are stimulation by growth factors or inhibition by application of a drug. We
refer generally to each of these as 'treatment'.

Each of the treatments will be included in the final network as a node, e.g. stimulation
by the growth hormone \emph{EGF} is added to the data matrix as row with name \emph{EGF} (and thus
appears as node \emph{EGF} in the final network). The expression values for the stimuli nodes are 
set to $0$ in each column of the data matrix, but are never used in the algorithm and regarded as
dummy values. Effects originating in these nodes are estimated in the
reconstruction process.

To distinguish the different experimental conditions in the matrix, the
columns of the data matrix have to be named in the format \emph{treatment\_time}, where treatment
also can be a combination of several treatments, e.g. stimulation by \emph{EGF} and simultaneous
inhibition by a drug \emph{X}. In this case, each stimulus has to be separated by an ampersand
(\emph{$\&$}). The time point is separated from the stimuli via an underscore character
(\emph{$\_$}), and can be on any scale (minutes, hours etc.). An example data matrix is
shown below. In this table, the dummy rows for the treatments are already included (rows EGF and X).
However, they are not mandatory as input to \emph{ddepn} and, if missing, will be generated 
automatically, only requiring the correct labeling of the columns. An example for a data matrix is 
given in table \ref{tab:exampledatamatrix}.

%\begin{scriptsize}
\begin{table}[h]
\begin{tiny}
\label{tab:exampledatamatrix}
\begin{tabular}{l|cccccccc}
	& EGF\_1& EGF\_1 & EGF\_2 & EGF\_2 & EGF\&X\_1 & EGF\&X\_1 & EGF\&X\_2 & EGF\&X\_2\\
	\hline
	EGF & 0 &  0 &  0  & 0 &  0 &  0  & 0 &  0 \\
	X   & 0 &  0 &  0  & 0  & 0 & 0 &  0  & 0 \\
	AKT & 1.45 & 1.8 &  0.99 & 1.6 & 1.78 & 1.8  & 1.56 & 1.58 \\
	ERK & 1.33 & 1.7 &  1.57 & 1.3  & 0.68 & 0.34 & 0.62 & 0.47 \\
	MEK & 0.45 & 0.8 &  0.99 & 0.6  & 0.78 & 0.8  & 0.56 & 0.58
\end{tabular}
\caption{Example data matrix for 3 nodes and 2 stimuli.}
\end{tiny}
\end{table}

%\end{scriptsize}

%\subsection*{\emph{stimuli}} 
% The user might pass an argument stimuli directly, describing the structure of the treatments,
% i.e. which treatments were done simultaneously. This is usually generated automatically using the 
% properly labeled column names of the data matrix and will be removed soon. However a short example
% exemplifies the use of \emph{stimuli}:\\
% Consider two experimental conditions, first stimulation by \emph{EGF}, second 
% simultaneous stimulation by \emph{EGF} and inhibition by \emph{X}. The \emph{stimuli} list holds 
% two elements:\\
% 
%\begin{verbatim}
% > stimuli
%[[1]]
%EGF 
% 1 
%
%[[2]]
%EGF  X 
% 1   2
%\end{verbatim}
%  
%Each element is a named numeric vector, holding the row numbers of the treatments in the data table 
%labeled by the names of the treatments.


\subsection{Reference network \emph{phiorig}} 

If desired, a reference network \emph{phiorig} can be given, used to compare the edges of the 
inferred network to it. The user \textit{must} ensure that all treatments are included as nodes,
since the inference will estimate effects from these.\\
The format of the network is an adjacency matrix, where each entry corresponds to an edge
from the node specified by the rowname to the node specified by the column name. Two types of
edges are allowed: $1$ for activation, $2$ for inhibition. $0$ means no edge between the pair of
nodes. An example network and corresponding adjacency matrix are shown in figure \ref{fig:examplenet}.

<<ExampleNet, echo=FALSE, eval=FALSE>>=
phi <- matrix(0,nrow=5,ncol=5,dimnames=list(c("EGF","X","AKT","MEK","ERK"),c("EGF","X","AKT","MEK","ERK")))
phi[1,c(3,4,5)] <- 1
phi[2,4] <- 2
phi[4,5] <- 1
layout(t(matrix(c(1,2))),widths=c(1,2))
plotdetailed(phi,stimuli=list(c(EGF=1,X=2)))
plotmatrix(phi,"phi")
@

\begin{figure}
\centering
%\includegraphics[width=8cm,height=5cm]{examplenet}
\includegraphics{examplenet}
\caption{Example network. Left: graph representation, Right: adjacency matrix}
\label{fig:examplenet}
\end{figure}

%\begin{verbatim}
%> phi
%    EGF X  AKT MEK ERK
%EGF  0  0   1   0   0
%X    2  0   1   0   0
%AKT  0  0   0   0   0
%MEK  0  0   1   0   0
%ERK  1  1   0   0   0
%\end{verbatim}

\subsection{Seed networks \emph{phi}}

Both the GA and MCMC sampler require single or multiple seed networks. If not given, an 
unconnected network is used as seed for each individual in the GA population or the start 
networks for each MCMC run, respectively. However, the user can provide own seed networks using 
the argument \emph{phi}. This can either be an adjacency matrix or a list of adjacency matrices.
Again, the treatments \textit{must} be included as nodes. If given a single adjacency matrix, it
is used as seed network for each of the individuals in the population of networks during the
genetic algorithm, or as independent seeds for parallel MCMC samplings. If given as list,
its length must equal the number of individuals in the population in the GA (specified by the 
function argument \emph{p}), or the number of independent runs in the MCMC sampler (specified by
the argument \emph{cores}). The format of the adjacency matrices is the same as for \emph{phiorig}.


\section{Prior knowledge inclusion}
\label{sec:prior}

Currently, two methods for using biological prior knowledge are implemented. We refer to the 
first as \emph{Laplace} prior (\cite{Froehlich2007,Werhli2007}), and to the second as
\emph{ScaleFree} prior (\cite{Kamimura}).

\subsection{Laplace prior}
\label{sec:laplace}
The laplace prior penalises deviations of edges in the inferred network from prior edge
confidences, which can be aquired from external network databases
(e.g. KEGG \cite{Kanehisa2008KEGG}). The package includes a snapshot of the KEGG
database (containing signalling and disease pathways), downloaded in October 2010.

<<kegggraphs, echo=TRUE, results=hide, eval=FALSE>>=
data(kegggraphs)
length(kegggraphs)
@
\begin{verbatim}
[1] 78
\end{verbatim}

The list \emph{kegggraphs} includes 78 elements, each of which has 3 members, a string \emph{id} 
specifying the KEGG pathway id, a graphNEL object \emph{g} and an adjacency matrix \emph{phi}.

<<kegggraphs name, echo=TRUE, results=hide, eval=FALSE>>=
names(kegggraphs)[1]
@
\begin{verbatim}
[1] "MAPK signaling pathway"
\end{verbatim}

<<kegggraphs element, echo=TRUE, results=hide, eval=FALSE>>=
kegggraphs[[1]]
@
\begin{verbatim}
$id
[1] "04010"

$g
A graphNEL graph with directed edges
Number of Nodes = 267 
Number of Edges = 882 

$phi
              hsa:5923 hsa:5924 hsa:11072 hsa:11221 hsa:1843...
hsa:5923             0        0         0         0        0...
...
\end{verbatim}

Each graphNEL object was converted to a detailed adjacency list phi (including inhibitions as entries
with value 2) using the 'kegggraph.to.detailed.adjacency' function:

<<Convert graphNEL to detailed adjacency matrix, echo=TRUE, results=hide, eval=FALSE>>=
kegggraph.to.detailed.adjacency(gR=kegggraphs[[1]]$g)
@

To obtain prior probabilites for each edge between all pairs of nodes present in \emph{kegggraphs},
one can follow the approach of \cite{Werhli2007} and count all pairs of nodes occuring in the 
reference networks. Further the number, how often each pair is connected by an edge is counted. The 
support for an edge is then the ratio of the number of edges divided by the number of pairs. If 
the edge is an inhibitory edge, the ratio is multiplied by $-1$, leaving a negative confidence score. 

This kind of prior matrix can the be used together with the \emph{laplaceinhib} prior:
In the function call to \emph{ddepn}, just pass the arguments \emph{B} and \emph{lambda} and
set \emph{usebics=FALSE} and \emph{priortype="laplaceinhib"} to use the laplace prior for the inference.

<<Use Laplaceinhib prior, echo=TRUE, results=hide, eval=FALSE>>=
ddepn(data, lambda=0.01, B=B, usebics=FALSE,
		priortype="laplaceinhib")
@

If no information is available on the type of the edges in the reference networks, one should use the 
\emph{laplace} prior type. Here, all entries in $\mathcal{B}$ are positive (describing the belief 
in existence of an edge), and for the prior calculation, the edge type in the inferred network is 
ignored, too.

<<Use Laplace prior, echo=TRUE, results=hide, eval=FALSE>>=
ddepn(data, lambda=0.01, B=B, usebics=FALSE,
		priortype="laplace")
@


%We compute a normalising factor $Z$ as described in \cite{Werhli2007} for a given hyperparameter 
%$\lambda$. The prior probability density for each edge between node $i$ and $j$ is calculated
%as laplace type distribution:
%The prior density for each edge between node $i$ and $j$ is calculated as laplace type distribution:

%\[
%P(\Phi_{ij}|\lambda) = \frac{1}{2\lambda}exp(\frac{-|\Phi_{ij} - \mathcal{B}_{ij}|^\gamma}{\lambda})
%\]

%In the function call to \emph{ddepn}, just pass the arguments \emph{B} and \emph{lambda} and
%set \emph{usebics=FALSE} and \emph{priortype="laplaceinhib"} to use the laplace prior for the inference.
%

%<<Use Laplace prior, echo=TRUE, results=hide, eval=FALSE>>=
%ddepn(data, lambda=2, B=B, usebics=FALSE, priortype="laplaceinhib")
%@



\subsection{ScaleFree prior}
\label{sec:scalefree}

According to \cite{Kamimura} we set up a prior distribution that penalises high node degrees
in the inferred network. The assumption is that for biological networks the degree of a 
node follows a power law distribution, i.e. the probability of seeing $k$ nodes follows 

\[P(k) \propto k^{-\gamma}.\]

We set up the prior distribution as described in \cite{Kamimura}. To use the ScaleFree prior,
just pass the arguments \emph{gam} (the exponent $\gamma$), \emph{it} (the number of 
permutations) and factor \emph{K} to the function call of \emph{ddepn}, and again set 
argument \emph{usebics=FALSE}.

<<Use ScaleFree prior, echo=TRUE, results=hide, eval=FALSE>>=
ddepn(data, gam=2.2, it=500, K=0.8, usebics=FALSE)
@
 

\section{Use cases for GA and MCMC inference}
\label{sec:usecases}

This section shows the various types of calls to \emph{ddepn} with all of the different settings
(inference type, prior type).\\

\subsection{Data generation:}
<<Data generation,echo=TRUE, eval=FALSE>>=
library(ddepn)
set.seed(12345)
n <- 6
signet <- signalnetwork(n=n, nstim=2, cstim=0, prop.inh=0.2)
net <- signet$phi
stimuli <- signet$stimuli
weights <- signet$weights
dataset <- makedata(net, stimuli, mu.bg=1200, sd.bg=400, mu.signal.a=2000, sd.signal.a=1000)
data <- dataset$datx

# netga arguments
minetga <- 15
p <- 30
q <- 0.3
m <- 0.8

# mcmc arguments
mimcmc <- 1000
burnin <- 100
lambda <- 5

#laplace prior arguments
lambda=0.01
# for priortype="laplace": use original net as prior and reset inhibition
# edges, such that their type is ignored
Bpos <- net
Bpos[Bpos==2] <- 1
# for priortype="laplaceinhib": use original net as prior and reset
# inhibition edge entries in B from 2 to -1
B <- net
B[B==2] <- -1

# scale free prior arguments
gam <- 2.2
it <- 500
K <- 0.8

@


\subsection{GA, use BICs optimisation and no prior}
<<GA BIC,echo=TRUE,eval=FALSE>>=
ret <- ddepn(data, phiorig=net, inference="netga",
             maxiterations=minetga, p=p, q=q, m=m,
             usebics=TRUE)
@
\subsection{GA, use laplaceinhib prior}
<<GA Laplace,echo=TRUE,eval=FALSE>>=
ret <- ddepn(data, phiorig=net, inference="netga",
		maxiterations=minetga, p=p, q=q, m=m,
		usebics=FALSE, lambda=lambda, B=B, priortype="laplaceinhib")
@
\subsection{GA, use laplace prior}
<<GA Laplace,echo=TRUE,eval=FALSE>>=
ret <- ddepn(data, phiorig=net, inference="netga",
		maxiterations=minetga, p=p, q=q, m=m,
		usebics=FALSE, lambda=lambda, B=Bpos, priortype="laplace")
@
\subsection{GA, use scalefree prior}
<<GA ScaleFree,echo=TRUE,eval=FALSE>>=
ret <- ddepn(data, phiorig=net, inference="netga",
		maxiterations=minetga, p=p, q=q, m=m,
		usebics=FALSE, gam=gam, it=it, K=K, priortype="scalefree")
@
\subsection{MCMC, use laplaceinhib prior}
<<MCMC Laplace,echo=TRUE,eval=FALSE>>=
ret <- ddepn(data, phiorig=net, inference="mcmc",
		maxiterations=mimcmc, burnin=burnin,
		usebics=FALSE, lambda=lambda, B=B, priortype="laplaceinhib")
@
\subsection{MCMC, use laplace prior}
<<MCMC Laplace,echo=TRUE,eval=FALSE>>=
ret <- ddepn(data, phiorig=net, inference="mcmc",
		maxiterations=mimcmc, burnin=burnin,
		usebics=FALSE, lambda=lambda, B=Bpos, priortype="laplace")
@
\subsection{MCMC, use scalefree prior}
<<MCMC ScaleFree,echo=TRUE,eval=FALSE>>=
ret <- ddepn(data,phiorig=net, inference="mcmc",
		maxiterations=mimcmc, burnin=burnin,
		usebics=FALSE, gam=gam, it=it, K=K, priortype="scalefree")
@



%If desired, a reference network \emph{phiorig} can be given. This is used to compare edges to the inferred
%network in the algorithm and calculate sensitivity and specificity statistics. A seed network \emph{phi} can
%also be provided, which is used as starting point for inference by the genetic algorithm  


\section{Application to example data set from HCC1954 cell line}
\label{sec:application}
We show how we apply DDEPN to a data set from HCC1954 breast cancer cell line. Phosphorylation
of 16 proteins was measured on Reverse Phase Protein Arrays. We performed time courses over
10 time points (0,4,8,12,16,20,30,40,50,60 minutes). Each measurement was replicated five times
independently and for each of these five biological replicates, spotting was done in 
triplicate for each time point, leaving 15 replicates per time point and protein. Cells 
were stimulated with the EGFR ligand Epidermal Growth Factor (EGF) and the ERBB3 ligand 
Heregulin (HRG), both as single treatments as well as combined treatment (EGF\&HRG). The data
is included in the dataset:

<<Get data,echo=TRUE,eval=FALSE>>=
data(hcc1954)
@

First, change the column names. These contain, which of the replicates belong to the same
biological replicate, which is not used in the inference in DDEPN. 

<<Get data,echo=TRUE,eval=FALSE>>=
dat <- format_ddepn(hcc1954)
@

Now start a genetic algorithm to infer a signalling network between the 16 phosphoproteins
(takes a while):

<<Infer signalling network,echo=TRUE,eval=FALSE>>=
mi=1000
p=500
q=0.3
m=0.8
ret <- ddepn(dat, phiorig=NULL, inference="netga",
		maxiterations=mi, p=p, q=q, m=m,
		usebics=TRUE)
@



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

%\bibliographystyle{abbrv}
%\bibliography{references}

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

\begin{figure}[p]
\centering
\includegraphics{scorenetgatest}
\caption{Output plot of \emph{netga}. The left column shows the posterior, likelihood and 
prior traces (each point in the trace corresponds to the quantile of the respective score
that was specified by argument \emph{quantBIC}). The right column shows the traces of the 
posterior, likelihood and prior differences between iteration $i$ and $i+1$.}
\label{fig:scorenetgaoutput}
\end{figure}

\begin{figure}[p]
\centering
\includegraphics{mcmctest3_2_stats}
\caption{Output plot of \emph{mcmc\_ddepn}. The left column shows the posterior, likelihood and 
prior traces. The middle column shows the traces of the posterior, likelihood and prior
ratios between iteration $i$ and $i+1$. The right column shows the original and 
inferred networks, as well as a ROC curve for the comparison of both networks (only 
if a reference network is provided).}
\label{fig:resultmcmc}
\end{figure}

\clearpage

\section*{Session Information}

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


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

\bibliographystyle{plain}  % Style BST file
\bibliography{references}     % Bibliography file (usually '*.bib' ) 
\end{document}

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