%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LaTeX Template for a Method chapter of the Bacterial Molecular
% Networks volume to be published in the series Methods in Molecular
% Biology.
%
% Please rename this file as well as the .bib file according to the
% name of the first author.

%\documentclass[11pt,onecolumn]{article}
%\fontfamily{cmss}
\documentclass{bmc_article}
\usepackage{common/Sweave}
\usepackage{common/BMnetR}
%\documentclass{jobim}
\usepackage{url}
\usepackage{graphicx}
\usepackage{subfigure} 
\usepackage{cite}
\usepackage{longtable,lscape}
%\usepackage[latin1]{inputenc}

\begin{document}

%% Please systematically use the following styles for indicating file
%% names, program names, options,

\input{MiMB_BMN_commands}
\renewenvironment{Schunk}{\begin{shaded}\topsep=0ex\relax}{\end{shaded}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Title, authors and affiliations

\title{Analyzing biological data using R: methods for graphs and networks.} 
\author{Nolwenn Le Meur$^{1,2,*}$, Robert Gentleman$^2$
}

%% Autor's affiliations
\address{
\iid(*) corresponding author\\
\iid(1) IRSET EA SeRAIC 4427 - IRISA - Equipe Symbiose, Universit\'e de Rennes I
Campus de Beaulieu, 35042 Rennes cedex, France. email:nlemeur@irisa.fr\\
\iid(2) Fred Hutchinson Cancer Center Research, Program in
    Computational Biology,Fred Hutchinson Cancer Center Research, Program in
    Computational Biology,M2-B876, P.O. Box 19024, Seattle, WA 98109, USA. email:rgentlem@fhcrc.org\\
}

\maketitle


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ABSTRACT SECTION
%% Author's affiliations
\section*{Abstract}

R is a powerful language and widely used software tool for the
analysis and visualization of data.  Its core capabilities can be
extended through many different add-on packages.  Among the many
packages are some which offer a broad range of facilities for
analyzing statistical properties of graphs. This chapter will provide
a practical tutorial covering the use of R methods for graphs and
networks to examine biological data and analyze their topological and
statistical properties.

\subsection*{keywords}
%% 5 to 10 keywords
graph, random graph, network, statistic,  systems biology


\section{Introduction}
\label{section_intro}
\subsection{Context}

Recent advances in molecular biology and associated
technologies have expanded our opportunities to explore and understand the
complexity of biological systems at different levels; from
DNA to proteins to signaling networks and beyond. To integrate
these data, to build and explore biological networks, we now need
computational tools, efficient flexible data structures and
statistical methods designed with the particular challenges of systems
biology and the analysis of high-throughput data in mind.

For the past 10 years, the Bioconductor community has worked on
developing such infrastructure.  The tools come in the form of
different packages for R \cite{Rref}, a free software environment for
statistical computing and graphics.  Each package provides software,
or data, to carry out a specific set of tasks.  Developing workflows
that rely on multiple packages is an effective strategy for the rapid
development of tools for analyzing biological data.  Packages are
being developed to address the specific needs of systems biology
approaches, such as network modeling, simulation, formal verification
and graph visualization \cite{Huber:2007}.  For instance, Castelo and
Roverato developed a R package providing reverse engineering
procedures for molecular network discovery from microarray expression
data \cite{Castelo:2009}.  Similarly, using graph theoretical
concepts, Le Meur and Gentleman \cite{LeMeur:2008} identified
multi-protein complexes and pairs of multi-protein complexes that
share an unusually high number of synthetic genetic interactions.
Network data computed and analyzed in R can directly be visualized
within R using for instance the \Rpackage{graph} and
\Rpackage{Rgraphviz} Bioconductor libraries that will be discussed in
this chapter, or \Rpackage{igraph}, an interesting and valuable R
library \cite{Gabor2006}. Additionally, the network data generated by
R can easily be exported to other visualization tools. For example,
Cytoscape \cite{shannon2003} can be seen as an analyst-friendly and
interactive alternative. Among many features, Cytoscape notably allows
zooming in and out or using \textit{bird's eye view} for browsing
large graphs. It also offers classical network analysis such as
selecting subnetworks based on graph properties (\textit{e.g.,} node
or edge annotation).

In this chapter, we review methods available in R for working with
graphs.  We discuss how to retrieve and integrate biological data to
build graphs of gene networks.  We also consider some of the
techniques available to explore and analyze the topological and
statistical properties of these graphs.  This chapter was written
using the Sweave system \cite{Sweave} , and all the code used to
produce the outputs, graphics, etc. can be retrieved from our on-line
supplements, located at
\url{http://www.bioconductor.org/pub/MiMB}.

\subsection{Case Studies}
Case Study 1: \textit{Bacillus subtilis protein-protein interaction
  data.}
\\
\textit{Bacillus subtilis} known as the hay bacillus or grass
bacillus, is a Gram-positive bacterium commonly found in soil. It may
contaminate food but rarely causes food poisoning. It is not
considered a human pathogen. \textit{B.  subtilis} and the enzymes it
produces have many applications such as immuno-stimulatory agents in
alternative medicine or as additives in laundry detergents. In this
chapter, we will explore some of the protein-protein interaction (PPI)
data available for \textit{B. subtilis}.
\\

\noindent Case Study 2: \textit{Response of Escherichia coli during an
  oxygen shift.}
\\
\textit{Escherichia coli} is a Gram-negative bacterium that is
commonly found in the lower intestine of warm-blooded organisms. Most
\textit{E. coli} strains are harmless, but some, such as serotype
O157:H7, can cause serious food poisoning in humans. \textit{E. coli}
can grow in aerobic as well as in anaerobic conditions (\textit{e.g.,}
the intestine).  In this case study, we will explore gene expression
results of various knock-out (KO) mutants under aerobic conditions to
characterize elements of the transcriptional network.

\section{Material}
\label{section_material}
We assume that the reader is acquainted with R and only give a quick
overview on the installation procedures and help functionality. For
more details about R and its programming language, we refer the reader
to the many manuals and courses freely available on the R and
Bioconductor websites
(\url{http://www.r-project.org};\url{ http://bioconductor.org/}).  In
addition there are many good books on R and Bioconductor such as
\cite{MASS, RBioinf, JMCBook, BiocCaseStudies}.

The Graphviz package (\url{http://www.graphviz.org/}) and the Boost
graph library (\url{http://www.boost.org/}) are additional tools
required to run graph examples of this chapter. Their installation
might not be trivial; it depends on the operating system. One should
be especially careful to the location of those libraries and the
definition of their path in the environmental variables. For details,
we refer the reader to the project websites or the Bioconductor
packages that depend on them, namely \Rpackage{Rgraphviz} and
\Rpackage{RBGL}.

\subsection{R and Bioconductor}
R is freely available and runs on a wide variety of commonly used
operating systems.  To install R, download the most recent version of
R from \url{http://www.r-project.org/}. Detailed instructions of
installation are available in the R FAQ and the R Installation and
Administration Manual on the R project website.  The base installation
of R consists of the programming language itself together with a
number of packages that are mainly oriented to statistical analyses.
R can be extended and its capabilities expanded through the use of
add-on packages that are notably available through the R and
Bioconductor websites.  Bioconductor is comprised of several hundred R
packages that provide solutions to a large number of bioinformatic and
computational biology problems.  The scope of Bioconductor packages is
broad covering many aspects of the analysis of microarrays but also
including tools for genomic annotation, high-throughput genotyping,
cell-based assays (\textit{e.g.}, RNAi, cytometry), handling graphs,
and more recently high-throughput sequencing.

\subsection{Using R and Bioconductor packages}
Starting R on Windows or OS X is usually accomplished by
double-clicking on the R application icon. On UNIX-like systems users
typically type \texttt{R} at a shell prompt. Because there are many
complex dependencies between packages and the version of R being used
the recommended way to install Bioconductor packages is to use the
\texttt{biocLite.R} installation script in a R session:
%
\begin{Schunk}
\begin{Sinput}
> source("http://bioconductor.org/biocLite.R") 
> biocLite() 
\end{Sinput}
\end{Schunk}
%
To find, download, and install add-on packages, we recommend that you
use built-in capabilities, either in the GUI or through the command
line.  On Windows and OS X packages can be installed \textit{via} a
menu.  You can also use the R function \Rfunction{install.packages} or
the \Rfunction{biocLite} script to download a specified list of
packages and their dependencies.  There is a difference between
installing and using a package. Typically you install a package once,
using the tools described above and use it multiple times.  For
example we could install the \Rpackage{graph} and \Rpackage{Rgraphviz}
packages using the code below.
%
\begin{Schunk}
\begin{Sinput}
> biocLite(c("graph", "Rgraphviz"))
\end{Sinput}
\end{Schunk}
%
To use those packages you must load them using the \Rfunction{library}
function:
\begin{Schunk}
\begin{Sinput}
> library("graph")
> library("Rgraphviz")
\end{Sinput}
\end{Schunk}
%

\subsection{R Documentation}
%
To get started it is important to know where you can find
documentation such as help for functions, classes, and concepts.  The
\Rfunction{?} operator followed by the name of a function will display
the manual page of this function. \Rfunction{class?foo} will show the
manual page for the \texttt{foo} class and most often the available
methods for this class.  The function \Rfunction{apropos} can be used
to find objects in the search path with names that partially match the
character string given as an argument to \Rfunction{apropos}. The
function \Rfunction{help.search} will search the help system for
documentation matching a given character string, and again partial
matching will be used.  You should compare the output for the three
commands below on your system.
%
\begin{Schunk}
\begin{Sinput}
> class?graphNEL
> apropos("graphNEL")
> help.search("graphNEL")
\end{Sinput}
\end{Schunk}
%
In addition, each Bioconductor package contains at least one
\textit{vignette}, which is a document that provides a task-oriented
description of the functionality provided by the package.  Vignettes
should contain executable examples and are intended to be used
interactively.

\subsection{R data structures}
The basic R data structure is the atomic vector. It contains an
indexed set of values that are all of the same type: logical, numeric,
complex or character. The numeric type can be further broken down into
integer, single or double types. In the code below we create a vector
named \Robject{x}, and while it may seem that we have mixed integer
and real types, the integers have been coerced to their real
representation, which can be seen in the output of the call to
\Rfunction{typeof}.
%
\begin{Schunk}
\begin{Sinput}
> x <- c(10, 5, 6, 1.5, 9)
> typeof(x)
\end{Sinput}
\begin{Soutput}
[1] "double"
\end{Soutput}
\end{Schunk}
%
Additional data structures are lists, arrays and data frames.  An
\Robject{array} is a vector with a dimension attribute and matrices
are a special variant where there are exactly two dimensions.  Arrays
are like atomic vectors in that all elements must be of the same basic
type.  A \Robject{list} is an object consisting of an ordered
collection of objects. It is different from atomic vectors because its
elements can be of any type, and for example a list could contain a
numeric vector, a logical value, and a matrix.  A \Robject{data.frame}
is much like a matrix in that it is two dimensional. Data frames are
widely used for holding sample level data since that is typically
rectangular (some variables for each sample or individual) but the
data types (columns) can be of mixed types (\textit{e.g.,} age, gender,
and so on).  

In the code below we create a matrix, \Robject{y}, of integers with 2
rows and 5 columns.  We also create a data.frame with two columns, one
contains numeric values and the other contains character values.
%
\begin{Schunk}
\begin{Sinput}
> y <-c(1:10) 
> dim(y) <- c(2, 5) 
> age <- c(33, 44, 55)
> gender <- c("M", "F", "F")
> myDF <- data.frame(age, gender)
\end{Sinput}
\end{Schunk}
%
While creating the data.frame with the default parameters, we allow
the conversion of the character vector to factors. Factors are used to
encode vectors as a category (or enumerated type). Levels are then an
optional vectors of the values that in our case might be taken by
\Robject{gender}.  The default is the unique set of values taken by
as.character(gender), sorted into increasing order of \Robject{gender}.
\begin{Schunk}
\begin{Sinput}
> is.factor(myDF$gender)
\end{Sinput}
\begin{Soutput}
[1] TRUE
\end{Soutput}
\begin{Sinput}
> levels(myDF$gender)
\end{Sinput}
\begin{Soutput}
[1] "F" "M"
\end{Soutput}
\end{Schunk}
%
Subsets of the elements of a vector may be selected using the square
bracket operator.  If the object being subscripted has multiple
dimensions, then each dimension can be subscripted independently with
arguments separated by commas.  If a subscript is not supplied (see
the examples for \Robject{y} and \Robject{myDF} in the code below)
then it is presumed that all elements are selected.  In the code below
we select the first 2 elements of the vector \Robject{x}, then the
first two columns of the matrix \Robject{y}, and finally the first row
(patient) of the data.frame \Robject{myDF}.
%
\begin{Schunk}
\begin{Sinput}
> x[1:2]
> y[,1:2]
> myDF[1,]
\end{Sinput}
\end{Schunk}
%
R also supports a number of paradigms for object oriented programming.
Object-oriented programming provides a coherent way to describe and
manipulate rich data structures. For instance in the
\Rpackage{Biobase} package, the \Rclass{ExpressionSet} class, designed
to standardize data structures to represent DNA microarray data, has
become the basic data structure for most gene expression data analysis
packages.  Similarly, recent developments propose classes and basic
methods to support many other types of high-throughput data such as
cell-based assays \cite{Boutros2006,Hahne:2009}, genotyping, or DNA
sequencing \cite{Morgan2009}.

More generally, we distinguish packages that are mainly infrastructure
(sets of tools that can be used to create other pieces of software)
and packages that are designed to provide an end-user application.
For instance, the \Rpackage{graph} and \Rpackage{Rgraphviz} packages
are infrastructure, while the \Rpackage{GOstats} and
\Rpackage{apComplex} are specific applications built on that
infrastructure.  In addition there is a large set of packages that
provide annotation for whole organisms and a number of real data sets
are also available as self-contained packages suitable for method
development and testing.

In the following sections, we will make use of the packages listed in
Table~\ref{tab1}.  We will focus on data structures suited to graphs
and demonstrate their usefulness to manipulate and represent
biological networks.

\section{Methods}
\label{section_methods}
\subsection{Working with graphs in R}
\subsubsection{Graph types}
There are many equivalent ways to represent a graph. One is as a pair
of sets, $G = (V , E)$ where $V$ is the set of nodes (or vertices) and
$E$ is the set of edges.  An edge is used to represent a particular
type of relationship (interaction) between two nodes.  There are many
different types of graphs and their definitions primarily depend on
their edge and node properties.  For instance, a graph can be
undirected or directed.  For an undirected graph the relationship
determined by the edges is symmetric, and hence does not have direction
assigned to its edges like \Robject{g1} in
Figure~\ref{fig:rendergraph}.  For a directed graph (or
digraph) the relationship is asymmetric, and node A can have an edge
to node B, but B need not have an edge to A.  These relationships are
thought of as being directed and are often represented by edges with
arrowheads indicating the direction, see panel a) of
Figure~\ref{fig:Cliques} for an example.  A directed graph with no
directed cycles is called a \textit{directed acyclic graph}, or DAG.

The \textit{degree} of a node is the number of edges that are incident
on that node.  If the graph is directed then we can consider in-degree
and out-degree separately.

In data analysis, directed graphs are useful for representing
interactions, or experiments, that are asymmetric such as protein
interaction experiments where one protein is the \textit{bait} and
other proteins are \textit{prey}.  In many cases the proteins need to
be modified to perform the assay and it is not uncommon to observe
that protein A interacts with B, when A is the bait, but B fails to
interact with A when B is the bait.  Another situation where directed
graphs may be helpful is in modeling the relationship between
transcription factors (TFs) and their targets, we explore this use in
our second case study below.

Special types of graphs can be used and developed to
answer specific biological questions. For instance DAGs are used by
the Gene Ontology consortium \cite{GO:2000} to store and explore gene
annotation data.  One useful type of graph is one that represents the
output of a clustering algorithm.  We call such graphs \textit{cluster
  graphs} and they have the property that nodes are separated into
groups or clusters and within a cluster all nodes are connected (a
complete graph) but between clusters there are no edges.  A network is
a type of graph that can be used to represent flow.  In particular a
network is a directed graph with a vertex, called the source, whose
in-degree is 0, a vertex, called the sink, whose out-degree is 0, and
there can be a numerical capacity associated with each edge in a network.

\subsubsection{Building and exploring graphs}
There are many ways to computationally represent graphs.  Some data
structure are designed for efficiency of different types of access or
manipulation while others are efficient for handling very large, but
sparse, graphs.  The \Rpackage{graph} package provides different
representations such as matrix representations and a node and an edge
list representation.  A number of methods are available to coerce
graphs between the different representations.  There are other
packages, such as \Rpackage{hypergraph} that handle even more
specialized operations or types of graphs.

The \Rpackage{graph} package offers 3 basic strategies for building
random graphs.  The \Rfunction{randomEGraph} function generates a
random edge graph. In this graph edges are either randomly generated
according to a specified probability, or the user can specify a number
of edges that are then randomly assigned. \Rfunction{randomNodeGraph}
builds a random graph with a pre-specified node degree distribution.
And for the \Rfunction{randomGraph} function the user provides both
the node labels and a factor with some fixed number of levels. Each
node is randomly assigned levels of the factor and edges are created
between nodes that share the same levels of the factor.  As an
example, we create a random graph \Robject{g1} using the
\Rfunction{randomEGraph} function from the \Rpackage{graph}
package. \Robject{g1} is an undirected graph composed of 15 nodes
(labeled by the first 15 letters of the alphabet) and 20 edges.
%
\begin{Schunk}
\begin{Sinput}
> set.seed(123) 
> g1 <-  randomEGraph(LETTERS[1:15], edges= 20) 
> class(g1)
\end{Sinput}
\begin{Soutput}
[1] "graphNEL"
attr(,"package")
[1] "graph"
\end{Soutput}
\begin{Sinput}
> g1 
\end{Sinput}
\begin{Soutput}
A graphNEL graph with undirected edges
Number of Nodes = 15 
Number of Edges = 20 
\end{Soutput}
\end{Schunk}
%
The \Rpackage{graph} package also offers a collection of functions and
methods to convert various forms of matrices into graph objects
(\textit{e.g.,} \Rfunction{ftM2graphNEL}, \Rfunction{aM2bpG}). 


Other graph-related packages provide methods for computing graph
properties.  Some packages are aimed at specific areas of application
such as the \Rpackage{ppiStats} package which is aimed at PPI data and
the \Rpackage{rsbml} package which provides support for SBML modeling
(Systems Biology Markup Language \url{http://www.sbml.org}).  The
\Rpackage{RBGL} package is an interface from \Rpackage{R} to the Boost
Graph Library \cite{Siek}. It is the primary source of algorithms to
explore and analyze graphs.  For instance, it contains methods to
search for shortest paths, connected components, and cliques.  As an
example, if a graph represents a metabolic network, we might be
interested in computing the shortest-path to generate a product
downstream of a cascade of interlacing reactions. Alternatively, for a
PPI network or for a gene expression network we might want to retrieve
the connected components and cliques.

In the code below, we determine the connected components in the random
graph \Robject{g1} using the \Rfunction{connectedComp} function. Then
once we know which nodes are connected we ask about paths in the
graph, from one node to another.
%
\begin{Schunk}
\begin{Sinput}
> library("RBGL")
> s1 <-  connectedComp(g1)
> s1
\end{Sinput}
\begin{Soutput}
$`1`
 [1] "A" "B" "C" "D" "E" "G" "H" "I" "J" "K" "L" "M" "N"
[14] "O"

$`2`
[1] "F"
\end{Soutput}
\begin{Sinput}
> sp.between(g1, "A", "O")
\end{Sinput}
\begin{Soutput}
$`A:O`
$`A:O`$length
[1] 2

$`A:O`$path_detail
[1] "A" "D" "O"

$`A:O`$length_detail
$`A:O`$length_detail[[1]]
A--D D--O 
   1    1 
\end{Soutput}
\end{Schunk}
%
We can extract the large connected component using the
\Rfunction{subGraph} function and then perform other computations
using just that subgraph.  In the code below we look for a solution to
the minimum cut problem. In graph theory, a cut is a partition of the
vertices of a graph into two non-empty but disjoint subsets by cutting
the minimum number of edges (for more details see the \Rfunction{minCut}
man page).
%
\begin{Schunk}
\begin{Sinput}
> gsub <- subGraph(s1[[1]], g1)
> minCut(gsub)
\end{Sinput}
\begin{Soutput}
$mincut
[1] 1

$S
[1] "B"

$`V-S`
 [1] "A" "C" "D" "E" "G" "H" "I" "J" "K" "L" "M" "N" "O"
\end{Soutput}
\end{Schunk}
%

\subsection{Graph layout}
In Bioconductor,  we have separated the process of drawing graphs into
two distinct phases. The first is \textit{layout} which involves
determining the locations of the nodes and edges, typically in a two
dimensional space.  This can then be followed by the process of
rendering the graph, which is the actual plotting of the graph.  There
are two reasons for separating these tasks.  First, we would like to
use R's graphics capabilities for rendering, so that graphs can easily
be intermingled with other plots and so that the full range of R's
advanced graphics can be used.  Second, we often want to plot graphs
based on the same set of nodes (and sometimes edges) multiple
times. Being able to control the location of the nodes and edges helps
us achieve that goal.  The main graph layout algorithms are provided by the
\Rpackage{Rgraphviz} package, which is based on the \textit{Graphviz
  project} \cite{graphviz} and others are available in the \Rpackage{RBGL} 
package.

A typical workflow of a graph layout is to pass a graph
object to the layout function (\Rfunction{layoutGraph}), which
returns another graph object containing all the necessary information
for subsequent rendering.  Figure~\ref{fig:rendergraph}(a) shows the
results of the rendering function (with default parameters) used
below to represent our random graph \Robject{g1}.
%
\begin{Schunk}
\begin{Sinput}
> g1 <- layoutGraph(g1) 
> renderGraph(g1)
\end{Sinput}
\end{Schunk}
%
There is a hierarchy to set rendering parameters for a graph. The
levels of this hierarchy are the session, the rendering operations,
and the individual nodes and edges. At the session level, defaults
parameters are used if not set somewhere further down the
hierarchy. For single rendering operation, \textit{i.e.,} a single
call to \Rfunction{renderGraph}, defaults can be set using its
\texttt{graph.pars} argument.
%
\begin{Schunk}
\begin{Sinput}
> renderGraph(g1, 
+             graph.pars=list(edges= list(col= "lightblue", 
+                               lty= 1))) 
\end{Sinput}
\end{Schunk}
%
All parameters set using \Rfunction{renderGraph}'s \texttt{graph.pars}
argument are transient, whereas setting session-wide parameters, using
the function \Rfunction{graph.par} will affect all subsequent
rendering operations.  Figure~\ref{fig:rendergraph}(b) is the result
of modifying the \texttt{graph.par} as follow:
%
\begin{Schunk}
\begin{Sinput}
> graph.par(list(nodes = list(fill= "lightgray", 
+                  textCol= "blue", cex= 1), 
+                edges=list(col= "red", lty= 2))) 
> renderGraph(g1)
\end{Sinput}
\end{Schunk}
%
Parameters for individual nodes or edges in a graph object can also be
set using the \Rfunction{nodeRenderInfo} and
\Rfunction{edgeRenderInfo} functions. These changes will be retained
in all subsequent layout or rendering operations of that particular
graph.
\begin{Schunk}
\begin{Sinput}
> nodeRenderInfo(g1) <- list(fill= c(A= "lightyellow", 
+                              B= "lightblue")) 
> renderGraph(g1) 
\end{Sinput}
\end{Schunk}
%
There are many options for layout type, node color, node shape, edge
representation etc., and we refer the reader to the vignette of the
\Rpackage{Rgraphviz} package for further details.

\section{Case studies}
\subsection{Case study 1: \textit{B. subtilis} protein-protein
  interaction data.}
In this case study we are interested in \textit{B. subtilis} PPI
data. We will focus on physical binary interactions between protein
pairs obtained via a yeast 2-hybrid (Y2H) assay \cite{Bruckner2009}
and complex co-membership between sets of proteins detected via
Affinity Purification - Mass Spectrometry (AP-MS) \cite{Wingren2009}).

It is worth briefly commenting on the data that arise from these two
assays.  For Y2H both bait and prey proteins need to be modified. One
of the two is modified to contain a DNA binding domain, while the
other is modified to contain an activation domain.  A typical
experiment involves a generally smaller set of bait proteins crossed
with a generally larger set of prey.  The resulting data are usually
processed by the experimenter to be binary and usually only positive
results are reported.  Thus the data generally consist of a set of
bait-prey pairs.  For these experiments both the set of baits that
were used and the set of prey they were targeted to are known, but are
seldom reported.  For AP-MS like experiments only the bait is modified
and hence the potential set of prey is very large. Each
\textit{pull-down} provides data on the set of interactors with the
bait protein (the protein that was pulled-down).  For these
experiments the set of baits is known, but the set of prey is
typically the set of expressed proteins in the tissue being assayed.


Curated, PSI-MI 2.5 XML, \textit{B. subtilis} PPI data are available
from IntAct (\url{http://www.ebi.ac.uk/intact/}).  The PSI-MI 2.5 XML
schema is a recent initiative by the \textit{Protein Standardization
  Initiative} to formally standardize the way biologists should report
molecular interaction data.  In Bioconductor, the \Rpackage{RpsiXML}
package serves as a programmatic interface between the R statistical
environment and the PSI-MI 2.5 XML files. It implements the PSI-MI 2.5
standard and currently supports 8 databases (unfortunately each
database has its own implementation of the PSI-MI 2.5 standard).  The
code chunk below demonstrates how to download and parse the
\textit{B. subtilis} PPI data using \Rpackage{RpsiXML}.  First, the
URL is specified, then that is used to open the connection and
download the text and finally that text is parsed to extract the
different components.
%
\begin{Schunk}
\begin{Sinput}
> library("RpsiXML") 
> site <- "ftp://ftp.ebi.ac.uk/pub/databases/intact/current/psi25"
> burl <- paste(site, "species/bacsu_small.xml", sep="/")
> bfile <- paste(readLines(url(burl)), collapse= "")
> bacsu <- parsePsimi25Interaction(bfile, 
+                                  INTACT.PSIMI25, 
+                                  verbose= FALSE) 
\end{Sinput}
\end{Schunk}
%
The R object \Robject{bacsu} is of class
\Robject{psimi25InteractionEntry}. You can access its attributes using
the associated methods like \Rfunction{interactions}, that lists each
interaction, or \Rfunction{organismName}, that returns the name of the
organism (see man page for the complete list of methods).

To explore the data, we might first want to understand what types of
interactions are present and in what quantities, as well as how many
different experiments are represented.  The interaction types refer to
the methods used to detect the interactions as stored in the IntAct
database like 'Y2H', representing binary physical interactions, or
'pull-down', representing complex co-membership (See the IntAct
annotation rules for more details
\url{http://www.ebi.ac.uk/~intact/site/doc/IntActAnnotationRules.pdf}). The
code below shows how to tabulate that information.
%
\begin{Schunk}
\begin{Sinput}
> ii = interactions(bacsu)
> table(sapply(ii, function(x) x@interactionType))
\end{Sinput}
\begin{Soutput}
            2 hybrid    adenylate cyclase 
                  41                    1 
      affinity chrom       anti bait coip 
                   7                    1 
       anti tag coip  confocal microscopy 
                   3                    1 
           crosslink  electron microscopy 
                   1                    4 
                emsa fluorescence imaging 
                   4                    3 
                fret     light scattering 
                   1                    2 
protein kinase assay            pull down 
                   5                   12 
                 tap  transcription compl 
                   4                    3 
   x-ray diffraction 
                   4 
\end{Soutput}
\begin{Sinput}
> table(sapply(ii, function(x) x@expPubMed))
\end{Sinput}
\begin{Soutput}
10997907 11722727 12065423 12100548 12359875 12644235 
       1        4        5        3        1        1 
12738765 14711369 14993308 15045029 15104138 15362849 
       1        3        9        4       13        1 
15491161 16796675 17662947 17803906 17994626 18854156 
       4       19        7        8        3        4 
19450516 19800110 
       1        5 
\end{Soutput}
\end{Schunk}
%
We can see that 41 of the interactions are from Y2H screens
and 12 are from AP-MS experiments.  From the second
call to \Rfunction{table} we see that 18 different PubMed citations
are involved and that the two largest sources of data are 13
interactions from PMID:15104138 and 19 interactions from
PMID:16796675.

We next can verify the quality of the downloaded interactions. For
instance, we can create a function that test whether Y2H dataset
have all their bait and prey information and apply this function to
the entire list of interactions.
\begin{Schunk}
\begin{Sinput}
> filterfun <- function(x) { 
+   x@interactionType=="2 hybrid" && is.na(x@bait) && is.na(x@prey)
+ }
> table(sapply(ii, filterfun))
\end{Sinput}
\begin{Soutput}
FALSE  TRUE 
   91     6 
\end{Soutput}
\end{Schunk}
%
%
The result indicates that there may be some issues: 6
out of 41 have missing information.

More generally quality assessment (QA) is an crucial processing step
in (any) data analysis and exploration. Experimental data (especially
high-throughput) are often biased due to systematic or stochastic
errors \cite{Schuchhardt2000, Gentleman2007}. QA allows identifying
biases (and potentially their causes) and minimizing the effects by
correcting or removing erroneous data. Several Bioconductor packages
offers QA methods specific to a variety technology. For example, the
\Rpackage{ppiStats} provides tools for the analysis of Y2H and AP-MS
data.  The package is specifically designed to summarize graph data,
identify systematically biased baits, and estimate stochastic error
probabilities \cite{Chiang2007, Chiang2009}.  In particular, its
authors \cite{Chiang2007} introduced the concepts of viable baits
(VBs) and viable prey (VP) to help take into account the problem of
untested interactions versus tested but not observed interactions.
Experimental results are typically only reported for proteins that are
observed in at least one interaction, leaving uncertainty as to which
proteins were actually tested.  VB is thus the set of bait proteins
that successfully detect at least one interaction partner, VP is the
set of all proteins detected by at least one bait, and VBP is the set
of proteins that are both viable baits and viable preys. In the
following paragraph we will get an overview of the usefulness this
package.

Upon reading the PPI data into R, one intuitive exploratory approach
is to represent it as a graph.  To extract a graph from XML data,
\Rpackage{RpsiXML} provides 2 methods: \Rfunction{psimi25XML2Graph}
that generates one large graph composed of all the bait-prey data
across all experiments while \Rfunction{separateXMLDataByExpt} yields
a list of graphs indexed by the PubMed ID.  Here we choose the latter
as the data originate from many various experimental techniques
(mostly low throughput) and it might be misleading to mix them
together.  For similar reasons, the \Rfunction{separateXMLDataByExpt}
function does not mix interaction data based on direct, e.g. binary
physical interactions, with that obtained from indirect methods, such
as AP-MS.  Here we ask for \texttt{indirect}, e.g. AP-MS,
interactions.
%
\begin{Schunk}
\begin{Sinput}
> graphs <- separateXMLDataByExpt(bfile, INTACT.PSIMI25, 
+                                 type="indirect", directed= TRUE,
+                                 abstract= TRUE, verbose= FALSE) 
> names(graphs) 
\end{Sinput}
\end{Schunk}
%
This subset of the PPI data originated from 6 different works. One way
of getting a descriptive summary of these experiments is to use the
\Rfunction{createSummaryTables} function of the \Rpackage{ppiStats}
package. This function takes a list of directed graph objects and
creates a summary table based on the directed connectivity, as shown
in Table~\ref{tab:summaryTable}.
%
\begin{Schunk}
\begin{Sinput}
> library("ppiStats")
> directedGraph <- sapply(graphs, isDirected)
> aSummary <- createSummaryTables(graphs[directedGraph])
\end{Sinput}
\end{Schunk}
%
%
We next consider the paper 16796675 \cite{Ishikawa2006} in more
detail. First, we can read the abstract information.
%
\begin{Schunk}
\begin{Sinput}
> abstract(graphs[["16796675"]])
\end{Sinput}
\begin{Soutput}
An object of class 'pubMedAbst':
Title: A new FtsZ-interacting protein, YlmF,
     complements the activity of FtsA during
     progression of cell division in Bacillus
     subtilis.
PMID: 16796675
Authors: S Ishikawa, Y Kawai, K Hiramatsu, M
     Kuwano, N Ogasawara
Journal: Mol Microbiol
Date: Jun 2006
\end{Soutput}
\end{Schunk}
Then we verify the number and type of interaction reported.
\begin{Schunk}
\begin{Sinput}
> Ishikawa <-  sapply(ii, function(x) x@expPubMed == 16796675)
> sum(Ishikawa)
\end{Sinput}
\begin{Soutput}
[1] 19
\end{Soutput}
\begin{Sinput}
> IshikawaInts <- ii[Ishikawa]
> table(sapply(IshikawaInts, function(x) x@interactionType))
\end{Sinput}
\begin{Soutput}
 2 hybrid pull down 
        8        11 
\end{Soutput}
\end{Schunk}
%
They report 19 interactions, 11 detected by pull down (most likely
AP-MS technology) and 8 detected via Y2H.  We next extract the graph
for their data from the \Robject{graphs} object.  Recall that in
creating this we only accessed \texttt{indirect} interactions, so only
the AP-MS-like interactions are represented.
%
\begin{Schunk}
\begin{Sinput}
> IshikawaGraph <- graphs$`16796675`
> IshikawaGraph
\end{Sinput}
\begin{Soutput}
A graphNEL graph with directed edges
Number of Nodes = 15 
Number of Edges = 32 
[1] "psimi25Graph"
attr(,"package")
[1] "RpsiXML"
\end{Soutput}
\end{Schunk}
%
The resulting graph, \Robject{IshikawaGraph}, is a directed graph composed
of 15 nodes (representing proteins) and 32 edges (representing
observed indirect interactions).  The IntAct data reported 11
interactions, but recall that each such interaction is between one
bait and one or more prey.  The nodes in the graph correspond to the
unique set of proteins detected as either bait or prey, while the
edges indicate which baits detected which prey.  Thus, there were 32
interactions detected among 15 different proteins.

In a biological network it is often useful to know the number of
interactors associated with a node.  This information can be used for
a variety of purposes, such as monitoring flux between components of
the network, characterizing protein complexes, or in assessing the
quality of proteomics bait-prey technologies by modeling experimental
errors on node degree \cite{Scholtens:2008}.  This corresponds to
computing the degree of the nodes in the graph.  Since the graph is
directed the in-degree and out-degree are reported separately.
%
\begin{Schunk}
\begin{Sinput}
> degree(IshikawaGraph)
\end{Sinput}
\begin{Soutput}
$inDegree
O31728 O34894 P17865 P28264 P37469 P94542 P54420 P33166 
     4      4      4      4      0      0      2      2 
P45694 P80868 P21465 P39751 P06574 P24327 P06567 
     2      2      2      1      2      2      1 

$outDegree
O31728 O34894 P17865 P28264 P37469 P94542 P54420 P33166 
     3      3      3     11      8      4      0      0 
P45694 P80868 P21465 P39751 P06574 P24327 P06567 
     0      0      0      0      0      0      0 
\end{Soutput}
\end{Schunk}
%
Two proteins, P28264 and P37469 have large out-degrees, 11 and 8
respectively.  However these relationships are not reciprocated (the
in-degrees are 4 and 0, respectively).  In fact, in PPI experiments a
major issue is that not all proteins are employed as both bait and
prey and we cannot observe interactions that were not tested.
In the Ishikawa experiment only four proteins were observed as both
bait and prey (see \cite{Chiang2007} for terminology).
%
\begin{Schunk}
\begin{Sinput}
> idViableProteins(IshikawaGraph)$VBP
\end{Sinput}
\begin{Soutput}
[1] "O31728" "O34894" "P17865" "P28264"
\end{Soutput}
\end{Schunk}
%
The \Rfunction{assessSymmetry} function from the \Rpackage{ppiStats}
package can be used to test the symmetry of your data by calculating
the reciprocated degree (nr), the nonreciprocating out (no) and in
(ni) degrees.  Using this information and a Binomial error model you
can assess the $p$-value for the in and out degree of each protein
\cite{Scholtens:2008}.

In protein interaction graphs, an other question of special interest
is the identification of functionally related modules (\textit{e.g.,}
protein complexes). In graph analysis this might be achieved by
searching for cliques. A clique is a complete subgraph, \textit{i.e.},
there is an edge between every pair of vertexes.  Here we search for
the maximum clique, \textit{i.e.,} the largest subgraph(s), using the
\Rfunction{maxClique} function implemented in \Rpackage{RBGL}. Though
this algorithm needs to be applied on an undirected graph. We thus
extract the underlying structure of \Robject{IshikawaGraph} using the
\Rfunction{ugraph} function from \Rpackage{graph}.
\begin{Schunk}
\begin{Sinput}
> library("RBGL")
> xu <- ugraph(IshikawaGraph) 
> cls <- maxClique(xu)$maxCliques 
> cs <- sapply(cls, length) 
> maxC <- cls[cs== max(cs)]
> maxC
\end{Sinput}
\begin{Soutput}
[[1]]
[1] "P28264" "O34894" "P17865" "O31728" "P94542"
\end{Soutput}
\end{Schunk}
%
We identify a subgraph that contains the four VBPs plus one additional
protein.  Figure~\ref{fig:Cliques} illustrates the result of the
maximum clique algorithm.  We do note that a clique is a very extreme
notion of a cohesive subgroup and that there are many others that
might be more suitable when the data are measured with error, as is
the case for most PPI experiments.  In particular we note that
\Rpackage{RBGL} provides implementations of \textit{k-cliques},
\textit{k-cores} and many other concepts described in the social
network literature \cite{SNA}.  Also, \cite{Scholtens2005} consider
this problem and provide both an algorithm and an implementation via
the \Rpackage{apComplex} package.
Figure~\ref{fig:Cliques} Panel (a) displays all the
interactions observed by \cite{Ishikawa2006} using the AP-MS
technology. Panel (b) renders only the maximum clique.  Note, that
P9542 only has out edges indicating it was a bait.
%
%
%
Perhaps the next step in an analysis would be to identify and
understand the roles and function of the proteins in the clique.
The protein interaction data are recorded using UNIPROT IDs but most
papers will use common names, or symbols, so we need some way to
navigate between these.   For many important organisms the
Bioconductor project offers database-oriented packages that provide
maps between the different identifiers used for that organism.  In
addition one can rely on external resources such as BioMart
\cite{biomart} which can be accessed from within R via the 
\Rpackage{biomaRt} package.

Every analysis with \Rpackage{biomaRt} starts with selecting a BioMart
database and one of the associated datasets (see the
\Rpackage{biomaRt} vignette for details). Here, we query the
"uniprot-mart" database and its only available dataset, "UNIPROT",
with the 5 selected proteins.  We define the list of filters and
attributes we are interested in and use the \Rfunction{getBM}
to retrieve annotation information we would like.
%
\begin{Schunk}
\begin{Sinput}
> library("biomaRt")
> uniprot <- useMart("uniprot_mart")
> uniprot <- useDataset("UNIPROT",mart=uniprot)
> filters <- listFilters(uniprot)
> attributes <- listAttributes(uniprot)
> attributes[c(1:4,23), ]
> selectedProt <- maxC[[1]]
> getBM(attributes =  c(attributes[c(1:4, 23), 1]),
+   filters = "accession", values = selectedProt, mart = uniprot) 
\end{Sinput}
\end{Schunk}
%
The name and comments columns of the results show us that the selected
proteins are indeed tightly related physically or functionally. The
proteins O31728 (\textit{sepF}), P17865 (\textit{ftsZ}) and P28264
(\textit{ftsA}) are part of the divisome complex (or \textit{FtsZ
  complex}), while P94542 (\textit{zapA}) is a divisome activator and
O34894 (\textit{ezrA}) is a divisome inhibitor
(Figure~\ref{fig:Cliques}b). We also remark that the protein
\textit{sepF} corresponds to the gene name \textit{YlmF} used in
\cite{Ishikawa2006}.

\subsection{Case Study 2: \textit{Response of Escherichia coli during
    an oxygen shift}}

In this study case, the data consist of 43 microarrays for 6 mutants
in two conditions aerobic and anaerobic. The mutants are knockout
strains of 5 key transcriptional regulators of the oxygen response
network: $\Delta${\it arcA}, $\Delta${\it appY}, $\Delta${\it fnr},
$\Delta${\it oxyR}, $\Delta${\it soxS} and the double knockout
$\Delta${\it arcA}$\Delta${\it fnr}. These data are part a
genome-scale analysis \cite{Covert2004}.  We restrict our attention to
the effect of the single transcription factor knockout data in aerobic
conditions.  We realize that a comparison of aerobic to anaerobic
conditions would be of some importance, but such an analysis would be
well beyond the scope of this chapter, where our primary goal is to
inform the reader of the available tools for graph analysis and to
demonstrate their use in a setting that is not especially complex.
Even with that simplification, our analysis is incomplete and quite
cursory.  However, the reader could easily extend our example in many
interesting directions.

First, we consider a model that will help to identify the set of genes
that are regulated by our transcription factors.  The nodes in our
network will be the genes in \textit{E. coli} that were assayed. We
define an edge to exist between one of the five transcription factors
listed above and any other gene if the gene expression between the
wild type condition and the condition with the transcription factor
knocked out is significantly different.  Expression of the second gene
can either increase (in which case the TF was likely a repressor) or
it can decrease (in which case the TF was likely an enhancer).  We
note that the biological effect need not to be direct; that is, there is
no reason to assume that the TF which acts in \textit{cis}, could not
affect an intermediate regulator of a second gene.  One could
determine which TFs are acting in \textit{cis} by using either known
TF binding sites or by performing a ChIP-seq experiment. 

As a second type of analysis we use the genes identified by the
knock-out experiments to identify cellular modules regulated by those
transcription factors.  We used the biological process component of GO
\cite{GO:2000} and KEGG pathways (\url{www.genome.jp/kegg/}) to
identify the candidate pathways and sets of genes for this part of the
analysis.

To get started, we create a vector with the names of the KO genes.
\begin{Schunk}
\begin{Sinput}
> KOgene <- c("arcA", "appY", "oxyR", "soxS", "fnr")
\end{Sinput}
\end{Schunk}
%
Next we load the expression data that we formatted from the online
supplementary information of \textit{Covert et al.}  \cite{Covert2004}
we downloaded at
(\url{http://www.nature.com/nature/journal/v429/n6987/suppinfo/nature02456.html}).
\textit{Covert et al.}  provided expression estimates for each probe
based on dChip \cite{dChip}.  This is slightly problematic for our
purposes since dChip can give negative expression values, so we
renormalized the data using the VSN algorithm \cite{VSN} from the
\Rpackage{vsn} package.  The normalized data are provided on the web
page that accompanies this article
(\url{http://www.bioconductor.org/pub/MiMB}). The data were saved in a
*.rda file (specific binary format for external representation of R
objects). Once downloaded from the  and can be reloaded in R using the function \Rfunction{load}
and the location of the file (\textit{e.g.,}
\verb+load(path/to/bES.rda)+ ).%
\begin{Schunk}
\begin{Sinput}
> load("data/bES.rda")
> bES
\end{Sinput}
\begin{Soutput}
ExpressionSet (storageMode: lockedEnvironment)
assayData: 4345 features, 37 samples 
  element names: exprs 
phenoData
  sampleNames: aerobic_wild_rep1, aerobic_wild_rep2, ...
  , anaerobic_soxS_rep3  (37 total)
  varLabels and varMetadata description:
    strain: 
    growth: 
    genotype: 
    name: 
featureData
  featureNames: thrL_b0001_at, thrA_b0002_at, ..., lasT_
  b4403_at  (4345 total)
  fvarLabels and fvarMetadata description:
    genename: NA
experimentData: use 'experimentData(object)'
Annotation:  
\end{Soutput}
\begin{Sinput}
> pData(bES)[1:5,1:3]
\end{Sinput}
\begin{Soutput}
                     strain  growth genotype
aerobic_wild_rep1 wild type aerobic       wt
aerobic_wild_rep2 wild type aerobic       wt
aerobic_wild_rep3 wild type aerobic       wt
aerobic_appY_rep1    mutant aerobic     appY
aerobic_appY_rep2    mutant aerobic     appY
\end{Soutput}
\end{Schunk}
%
The object \Robject{bES} is an object of class \verb+ExpressionSet+, a
container that holds both gene expression data and the corresponding,
per sample, metadata.

Statistical testing is essential to evaluate contrasts between
experimental conditions.  Since there are relatively few arrays we use
the \Rpackage{limma} package to test for differentially expressed
genes between KO and wild-type strains. 

Using the \Robject{bES} data set, we first subset the data so that
only the arrays corresponding to the aerobic condition are retained
(in code not shown we have also removed the double knock-out arrays).
%
\begin{Schunk}
\begin{Sinput}
> library("limma")
> aerobia <- bES[, bES$growth=="aerobic"]
\end{Sinput}
\end{Schunk}
%
Every \Rpackage{limma} analysis starts with the definition of a design
matrix, which we define below. A design matrix is a matrix of
explanatory variables which here allows you to formalize your
experimental design (observed data) and biological question for
statistical analysis. In our example, it describes which samples have
been applied to each array. Complex experimental design implies
complex definition of design matrix, requiring some statistical
knowledge (for more details we refer the reader to the well-documented
vignette of the \Rpackage(limma) package).
%
\begin{Schunk}
\begin{Sinput}
> aGeno <- as.factor(pData(aerobia)$genotype)
> aDesign <- model.matrix(~ 0 + aGeno)
> colnames(aDesign) <- levels(aGeno)
\end{Sinput}
\end{Schunk}
%
In the following code we fit the model, define the contrast to test,
and apply the moderated empirical Bayes test with Benjamini Hocheberg
$p$-value adjustment. This last step is to adjust for multiple
testing concerns.
%
\begin{Schunk}
\begin{Sinput}
> afit <- lmFit(aerobia,  aDesign) 
> a.cont.matrix <- makeContrasts("appY"=appY-wt, 
+                                "arcA"=arcA-wt, 
+                                "fnr"=fnr-wt, 
+                                "oxyR"=oxyR-wt, 
+                                "soxS"=soxS-wt, levels=aDesign) 
> afit2 <- contrasts.fit(afit, a.cont.matrix) 
> afit2 <- eBayes(afit2) 
> aRes <- c()
> for(i in 1:ncol(a.cont.matrix)){
+   aDE <- topTable(afit2, coef= i, n=500, adjust="BH", 
+                   p.value=0.001)
+   intRes <- cbind(colnames(a.cont.matrix)[i], 
+                   aDE$genename , round(aDE$logFC, 2))
+   aRes <- rbind(aRes, intRes)
+ }
> colnames(aRes) <- c("KO","targets","logFC")
\end{Sinput}
\end{Schunk}
%
We select those genes (probes) with $p$-value $\leq 0.001$.  The
number of significantly differentially expressed genes for each TF KO
condition can be obtained using the code shown next.
\begin{Schunk}
\begin{Sinput}
> table(aRes[,1])
\end{Sinput}
\begin{Soutput}
appY arcA  fnr oxyR soxS 
  33    1   65  147    1 
\end{Soutput}
\end{Schunk}
%
We note that three of the TFs have a large number of targets
identified, while two (\textit{arcA} and \textit{soxS}) have only one
target each. This may be due to the rather stringent $p$-value threshold
we used or the biological functions of those genes that depend on
the oxygen condition.  Another interesting feature of the data is that
the transcription factor KOs $\Delta${\it appY}, $\Delta${\it fnr},
$\Delta${\it oxyR} seem to have similar effects on fairly large sets
of their targets suggesting that perhaps they affect similar sets of
intermediate (Figure~\ref{graphAero}).

We can display the graph of these interactions using the
\Rpackage{graph} and \Rpackage{Rgraphviz} packages.  The type of graph
that is most helpful is one that is essentially a bipartite graph.  A
bipartite graph is a graph where nodes can be separated into two sets
and all edges go from a member of one set to a member of the other
set.  The mappings between GO terms and genes form a bipartite graph,
as do the mappings between KEGG pathways and genes.  In this case the
mapping is from the five candidate TFs to their targets.  But it is
not really a bipartite graph since the TFs can target themselves
(self-loops).  However, as you will see, rendering the graph yields an
essentially bipartite graph view.  There are some initial efforts to
build specialized tools for bipartite graphs, or the equivalent
representation as \textit{hypergraphs} in the \Rpackage{hypergraph}
package, but much remains to be done.

In practical terms we will construct this graph using the information
from \Rpackage{limma} which tells us both the TFs and their targets.
This is what can be referred to as a \textit{from-to} representation
of a graph and we use the function \Rfunction{ftM2graphNEL} to
construct an instance of the \Rclass{graphNEL} class.   In the call we
specify that the log fold change should be included in the graph as
edge weights. 
\begin{Schunk}
\begin{Sinput}
> graphAERO <- ftM2graphNEL(aRes[, c("KO","targets")],  
+                         W=as.numeric(aRes[, "logFC"]), 
+                           edgemode="directed")
\end{Sinput}
\end{Schunk}
%
We remove self-loops and thus the \textit{arcA} node, which only
regulates itself. The fact that \textit{arcA} only regulates itself is
intriguing. However \textit{ArcA} is known to act primarily as a
negative transcriptional regulator under anaerobic conditions thus
might not be expressed (or at very low level) in aerobic condition,
explaining its auto-regulation and the mild effect of its deletion
\cite{iuchi1988}.
\begin{Schunk}
\begin{Sinput}
> aeroG <- removeSelfLoops(graphAERO)
> aeroG <- removeNode("arcA", aeroG)
\end{Sinput}
\end{Schunk}
% 
We use the fold change between KO and wild type to color the edges
(orange means higher expression in KO and blue means higher expression
in wild type).  In the code below we demonstrate how to layout this
graph using the tools in R.  First, the graph is laid out, then we
determine which of the edges correspond to positive log FC and which
to negative.  We then construct a vector named \Robject{col} that
contains the value \texttt{"orange"} for those edges that are higher
in KO, and \texttt{"blue"} otherwise. We apply that information to the
graph using the \Rfunction{edgeRenderInfo} function. Finally, for
readability we remove the labels from the TF targets using
the \Rfunction{nodeRenderInfo} function and we render the graph.
%
\begin{Schunk}
\begin{Sinput}
> aeroG <- layoutGraph(aeroG, layoutType="dot")
> col <- c(rep("orange", numEdges(aeroG)))
> names(col) <- edgeNames(aeroG)
> low <-  unlist(lapply(edgeWeights(aeroG), function(x) x < 0))
> col[low] <- "blue"
> edgeRenderInfo(aeroG) <-  list(col=col)  
> KOlabel <- nodes(aeroG)%in%KOgene
> nodeRenderInfo(aeroG)$labelY[!KOlabel] <- NA
> renderGraph(aeroG)
\end{Sinput}
\end{Schunk}
%
Somewhat surprisingly, there is a preponderance of orange edges and
indeed 158 genes have higher expression in the KO
strains while only 85 have lower expression
levels. Suggesting that these TFs were largely repressing gene
expression, or perhaps that their deletion caused other systems to
up-regulate.  Had these TFs been largely activating transcription then
we would have seen blue edges, showing that their deletion
corresponded to a decrease in expression.

To better understand this, we turned to the literature and found
evidence that these TFs have dual roles.  \textit{fnr} is the primary
transcriptional regulator mediating the transition from aerobic to
anaerobic growth through the regulation of hundreds of genes.
Generally speaking, this protein activates genes involved in anaerobic
metabolism and represses genes involved in aerobic metabolism which
may partly explained the up-regulation observed in aerobic conditions
\cite{salmon2003}. In \textit{E.coli}, \textit{oxyR} was until now
mostly described as an activator (except for the \textit{flu} gene
\cite{correnti2002}) but it has recently been characterized as a
repressor in several other bacteria (such as
\textit{Neisseria gonorrhoeae} and \textit{Deinococcus
  radiodurans}\cite{chen2008}).  \textit{appY} induces the expression
of energy metabolism genes under anaerobiosis, stationary phase, and
phosphate starvation \cite{brondsted1996}. Finally, \textit{soxS} is a
transcriptional activator of oxidative stress genes, which might
explain that in physiological aerobic condition its depletion has a
minor effect. 

Of course we can not rule out the possibility that there was somewhat 
more variation in some of the experiments than in others. This would make it
harder to detect differential expression under some conditions and hence
could reduce the target set.  A thorough analysis would engage this
idea and ensure that there were no obvious differences between the
conditions.

We can also use KEGG and GO categories to try to better understand
whether there are particular pathways or processes that are affected.
This analysis consists of testing for over-representation of
categories among genes in a selected gene set based upon the
Hypergeometric distribution \cite{falcon2007}. To this aim, we make
use of the methods in the \Rpackage{GOstats} and \Rpackage{Category}
packages. To perform this test we need the Entrez Gene identifiers of
each \textit{E. coli} gene on the array and in each set selected from
the \Rpackage{limma} analysis. The \Rpackage{org.EcK12.eg.db}
annotation package allows such mapping and more, providing various
annotations such as GO categories, REFSEQ or PMID. To map
\textit{E. coli} gene name symbols to their Entrez Gene Identifiers we
first select the gene name from the \Robject{bES} \verb+ExpressionSet+
and get the unique Entrez Gene Identifiers. We apply the
same procedure to the genes associated with each KO and perform
the GO and KEGG \Rfunction{hyperGTest} function to search for
over-representation with a $p$-value cut-off of $0.001$. 
We do note that this approach, while widely used is slightly problematic. 
Among the many issues that arise is the fact that we are combining genes that
are both up-regulated and down-regulated by the TF KO.
%
%
First we consider the results for the GO category analysis.  We
restricted our attention to the biological process (BP) ontology.
Very similar categories are identified for \textit{fnr},
\textit{oxyR}, and \textit{appY}.  The set of GO categories for
\textit{fnr} was a subset of those identified for \textit{appY}, with
\texttt{GO:0006066} being associated with \textit{appY} but not
\textit{fnr}.  Similarly the set of GO categories for \textit{appY}
was essentially a subset of those for \textit{oxyR}, with
\texttt{GO:0006066} and \texttt{GO:0006725} being associated with
\textit{appY} but not \textit{oxyR}.  These GO categories are labeled
\textit{alcohol metabolic process} and \textit{cellular aromatic
  compound metabolic process} respectively.  We list the complete set
of GO categories for \textit{oxyR} using the code below.
%
\begin{Schunk}
\begin{Sinput}
> sigGO <- lapply(GO, sigCategories)
> sapply(sigGO, length)
\end{Sinput}
\begin{Soutput}
arcA appY oxyR soxS  fnr 
   0   10   12    0    9 
\end{Soutput}
\begin{Sinput}
> lapply(GO["oxyR"], function(x) summary(x)[,c(2,3,7)])
\end{Sinput}
\begin{Soutput}
$oxyR
         Pvalue OddsRatio
1  2.158614e-17 72.053571
2  2.158614e-17 72.053571
3  2.810665e-13 23.839286
4  6.354012e-11 14.196429
5  4.449673e-09       Inf
6  4.449673e-09       Inf
7  4.449673e-09       Inf
8  3.379156e-07 49.753846
9  1.452331e-04  5.043394
10 3.069493e-04 19.283582
11 3.069493e-04 19.283582
12 9.962969e-04  3.042772
                                       Term
1                        iron ion transport
2                       metal ion transport
3                             ion transport
4             establishment of localization
5         enterobactin biosynthetic process
6                catechol metabolic process
7          siderophore biosynthetic process
8              peptide biosynthetic process
9                       aerobic respiration
10             thiamin biosynthetic process
11 thiamin and derivative metabolic process
12                    anaerobic respiration
\end{Soutput}
\end{Schunk}

Next we consider the results for the KEGG pathways. 
\begin{Schunk}
\begin{Sinput}
> sigKEGG <-lapply(KEGG, sigCategories)
> sapply(sigKEGG, length)
\end{Sinput}
\begin{Soutput}
arcA appY oxyR soxS  fnr 
   0    1    4    0    1 
\end{Soutput}
\end{Schunk}
%
We see that there four significant KEGG pathways for \textit{oxyR} and one
for each of \textit{appY} and \textit{fnr}.  All 3 TFs
share the one significant pathway, \textit{Biosynthesis of siderophore group 
nonribosomal peptides},  which is in accordance with the GO
categories identified.  We can then return to the original data to get
a sense of which of the TFs had genes involved in this pathway.
%
\begin{Schunk}
\begin{Sinput}
> gnodes <- nodes(graphAERO)
> geneInKEGG <- geneIdsByCategory(KEGG[[2]])$`01053`
> symbol <- unlist(mget(geneInKEGG, 
+                       org.EcK12.egSYMBOL, ifnotfound=NA))
> siderophore <- gnodes[gnodes%in%symbol]
> # level of differentially expressed gene
> as.data.frame(aRes[aRes[,2]%in%siderophore,])
\end{Sinput}
\begin{Soutput}
     KO targets logFC
1  appY    entC  2.16
2  appY    entE  2.06
3  appY    entB  2.09
4  appY    entF  2.36
5   fnr    entC  2.72
6   fnr    entE   2.6
7   fnr    entB  2.56
8   fnr    entF  2.93
9  oxyR    entC   3.1
10 oxyR    entE   2.9
11 oxyR    entB  2.91
12 oxyR    entF  3.43
\end{Soutput}
\end{Schunk}
%
We note that several genes involved in iron metabolism are up
regulated between WT and KO. And indeed it seems that KOs of
\textit{appY}, \textit{fnr}, and \textit{oxyR} seem to have similar
effects on a largely overlapping set of genes.  Microorganisms growing
under aerobic conditions need iron for a variety of functions
including reduction of oxygen for synthesis of ATP, reduction of
ribotide precursors of DNA, for formation of heme, and for other
essential purposes.  Accumulation of the intracellular siderophore has
also be shown to be up-regulated by oxidative stress in
\textit{Aspergillus nidulans} \cite{eisendle2006}, which underscores
the intertwining of iron metabolism and oxidative stress.

\section{Discussion }
\label{section_notes}
We have provided a fairly brief and concise overview of how some of
the many Bioconductor packages can be used to help analyze network
data.  The tools are quite varied, and certainly much remains to be
done, but most of the building blocks needed for a comprehensive tool
suite now exist and can fairly easily be used to accomplish complex
analyses.  

We believe that the current capabilities are extensive, but that
perhaps the most compelling aspect of these tools is the ease with
which they can be extended and adapted to different analysis
problems.  Interested readers can obtain our data and scripts from the
accompanying web-site. They can verify our results, perhaps point out
errors or bugs in our code, and then go on to carry out much more
detailed analyses of these.  We encourage them to do so.

\section*{Acknowledgments}
\label{section_acks}
Funding for this research was provided by grants from La Ligue Contre
Le Cancer (RAB08008NSA; project R08010NS) . We also gratefully
acknowledge the Intitut national de la sant\'e et de la recherche
m\'edicale (INSERM) for supporting this project.  Funding for RG was
provided by NIH grant P41 HG004059.

The authors would like to thank Wolfgang Huber, Tony Chiang, Shailesh
Date and Denise Scholtens for helpful comments and for providing a
stimulating research environment that has led to the many tools used
here and to the ideas that underlie this research.

\newpage
%\section*{References}

%\begin{thebibliography}{99}

%\bibliography{bilbiography}
%\bibliographystyle{unsrt}
\bibliography{NLeMeur_MiMB_MBN_2009_biblio.bib}
\bibliographystyle{bmc_article}

%\end{thebibliography}

\newpage
\section*{Figure captions}

\begin{figure}[htbp] 
  \begin{center}
    \subfigure{(a)\includegraphics[width=0.45\textwidth]{figure/Fig-rendergraph1}}
    \subfigure{(b)\includegraphics[width=0.45\textwidth]{figure/Fig-rendergraph2}}
    \caption{\label{fig:rendergraph} Flexible graph rendering with
      R. The graph was generated by the \Rpackage{graph} package.  In
      panel (a) the default parameters of the \Rfunction{layoutGraph}
      and \Rfunction{renderGraph} functions were used.  In panel (b)
      the node and edge parameters were modified using the
      \Rfunction{graph.par} function. This affected the node color,
      line color and type etc. }
    \end{center}
\end{figure}

\begin{figure}[htbp] 
  \begin{center}
    \subfigure{(a)\includegraphics[width=0.45\textwidth]{figure/Fig-Clique}}
    \subfigure{(b)\includegraphics[width=0.45\textwidth]{figure/Fig-CliqueAlone}}
    \caption{\label{fig:Cliques} Graphs of interacting
      proteins in \textit{B. subtilis}, during cell division, as
      reported by \cite{Ishikawa2006}.  Only those interactions
      detected by AP-MS are rendered here.   Panel a) shows all
      interactions and their directions (the arrows go from a bait
      protein to a prey protein).  In orange
      is the maximum clique found by the \Rfunction{maxClique}
      function of the \Rpackage{RBGL} package.  In panel (b) only the maximal
      clique is rendered.  The proteins O31728, P17865 and P28264 are part of the
      divisome complex; P94542 is a divisome activator; and O34894 is a
      divisome inhibitor.}
    \end{center}
\end{figure}

\begin{figure}[htbp] 
  \begin{center}
    \includegraphics[width=0.95\textwidth]{figure/Fig-graphAero}
    \caption{\label{graphAero} Unbalanced effect of the deletion of
      \textit{Escherichia coli} key transcriptional regulators of the
      oxygen response, in aerobic condition. In this figure the TF
      knock-outs (top) are linked to their targets
      (bottom). Self-loops have been removed (auto-regulation of the
      TFs) as well as the node for the TF \textit{arcA}, which only
      regulates itself and affects no other gene in the aerobic
      conditions. The edges are colored according to whether there is
      higher expression in WT than in KO (blue) or higher expression
      in KO than WT (orange).  There are a surprising number of orange
      edges. We can also see the high level of overlap in targets
      between TFs.}
    \end{center}
\end{figure}




%% 
%% Please read carefully the instructions about figure formatting and
%% numbering in the section ``Illustrations'' of the instructions to
%% authors (\file{Author\_Instructions\_Humana\_MiMB.doc}).
%% 
%% In LaTex, Figures are cross-referenced in the text by using the
%% command $\\$ref (Figure~\ref{fig_my_first_figure}).
%% 
%% \begin{figure}[h!]
%%   \caption{Legend of the first figure.}
%%   \label{fig_my_first_figure}
%% \end{figure}
%% 
%% \begin{figure}[h]
%%   \caption{Legend of the second figure.}
%%   \label{fig_my_second_figure}
%% \end{figure}
%% 
%% 
\newpage
\section*{Tables}
%% 
%% 
\begin{table}[!htbp]
  \centering
  \caption{\textbf{Bioconductor's graph-related packages}. Packages used in this chapter.}
  \label{tab1}
  \begin{center}
 \begin{tabular}{p{3cm}p{5cm}p{3cm}}
   \hline
   Package & Description & Type  \\
   \hline
   \Rpackage{graph} &  basic class definitions and graph handling capabilities & infrastructure \\ \hline
   \Rpackage{Rgraphviz} & interfaces R with the AT and T Graphviz library for plotting R graph objects from the graph package & infrastructure \\ \hline
   \Rpackage{RBGL} & interface to the graph algorithms contained in the \textit{BOOST} library & infrastructure\\ \hline
   \Rpackage{RpsiXML} & R interface to PSI-MI 2.5 files & infrastructure \\ \hline
   \Rpackage{ppiStats} & tools for the analysis of protein interaction data & application \\ \hline
   \Rpackage{limma} & Data analysis, linear models and differential expression for microarray data & application \\ \hline
   \Rpackage{biomaRt} & interface to the BioMart software suite & annotation \\ \hline
   \Rpackage{GOstats} & basic manipulation tools for graphs, hypothesis testing and other simple calculations on Gene Ontology. & application \\ \hline
   \Rpackage{Category} & statistical methods for performing category analysis.& application \\ \hline
   \Rpackage{org.EcK12.eg.db} & genome wide annotation for \textit{E. coli} strain K12, primarily based on mapping using ENTREZ GENE identifiers from NCBI & annotation\\
   \hline
 \end{tabular}
\end{center}
\end{table}


\input{table/table2}

\end{document}



