BACKGROUND
in the modern field of systems biology scientists aim to get insights into the architecture and behavior of complex cellular and genomic processes.
an important task in this context is the detection of novel interdependencies between gene products.
this insight into the genomic networks is an important step towards a better understanding of the functional aspects of a biological system and of great value for drug target identification at a later stage.
within this context modern dna microarray technology plays an important role.
in addition, the advent of rna silencing techniques has further increased its power by allowing the selective knock-down of certain genes of interest.
this may enable us to detect interdependencies between gene products on a non-transcriptional level.
the genes of interest are knocked down individually, and the respective downstream effects on gene expression are measured by using genome-wide microarrays.
by observing the nested structure of significant up or down regulations of affected genes, this may allow one to reverse engineer features of the upstream signaling pathway  <cit> .
in a recent work markowetz et al.
<cit>  introduced a method to reverse engineer the signaling pathway between perturbed genes using the nested structure of secondary downstream effects.
they developed a bayesian statistical framework, in which for a given network hypothesis one can calculate a score and thus can reduce the set of all possible networks to the most likely ones.
a severe limitation of this method lies, however, in the restriction to small networks of up to  <dig> genes, because the method completely enumerates all possible network hypotheses.
furthermore, a difficulty in the practical use is the required binary discretization of the data .
in our work we therefore aim to extend the framework by markowetz et al.
in order to make it practically applicable for a broader range of real life problems.
we are thereby motivated by biological experiments conducted in our department:  <dig> genes in the er-α pathway in human mcf- <dig> breast cancer cells were silenced via small interfering rnas and the effects on gene expression were subsequently measured on cdna microarrays.
our extensions of the original approach go in three directions: first, we introduce a way to omit the data discretization step needed via a calculation based on p-values instead, which is more suitable for our data and makes the whole framework more flexible .
second, we show how prior assumptions on the network structure can be incorporated into the network scoring scheme via techniques from regularization theory  <cit> .
third and most important, we develop and investigate methods to scale up the network inference to large scale networks.
for this purpose two approaches are considered: simulated annealing on a restricted set of possible networks and our so-called module networks, which build the complete network recursively from smaller pieces that are connected subsequently.
in order to validate these approaches we conduct studies on artificially created networks and show that module networks offer the highest sensitivity and specificity in the reconstruction of edges in the networks.
finally, we demonstrate the applicability of our approach to real data by inferring the complete  <dig> genes er-α signaling pathway network.
using a bootstrapping approach this reconstruction can be found with good statistical stability and hence seems to be reliable.
RESULTS
statistical inference framework for signaling pathways from rnai data
we start with a brief review of the statistical inference framework for signaling pathways by markowetz et al.
: in general within this framework one distinguishes between silenced genes  and genes showing a downstream effect .
it is assumed that each e-gene is attached to a single s-gene only .
knocking down a specific s-gene sk interrupts signal flow in the downstream pathway, and hence an effect on the e-genes attached to sk and all s-genes depending on sk is expected.
let us assume n knock-downs are performed and there exist m e-genes in total.
the outcomes of these experiments are summarized in an m × n data matrix d. according to bayes' formula a specific network hypothesis Φ ∈ { <dig> }n × { <dig> }n can be scored as:
p ∝ pp
the position of the e-genes is introduced as a model parameter Θ = {θi|θi ∈ { <dig> ..., n}, i =  <dig> ..., m}, i.e θi = j, if e-gene i is attached to s-gene j. assuming independence of the observations  di in the data matrix d  one can write down the conditional likelihood p|Φ, Θ) as:
p=∏i=1mp
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebarcqgg8bafcqqhmogrcqggsaalcqqhyoqucqggpaqkcqgh9aqpdaqewbqaaiabdcfaqjabcicaoiabdseaenaabaaaleaacqwgpbqaaeqaaogaeiifawnaeuopdykaeiilawcccigae8hude3aasbaasqaaiabdmgapbqabagccqggpaqkasqaaiabdmgapjabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@4a79@
it is furthermore assumed that all parameters θi are statistically independent, i.e.
p=∏i=1mp
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqqhyoqucqgg8bafcqqhmogrcqggpaqkcqgh9aqpdaqewbqaaiabdcfaqjabcicaoggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiifawnaeuopdykaeiykakcaleaacqwgpbqacqgh9aqpcqaixaqmaeaacqwgtbqba0gaey4diunaaaa@4506@
the likelihood p can now be written as:
p = ∫p pdΘ
=∏i=1m∑j=1npp
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqgh9aqpdaqewbqaamaaqahabagaemiuaalaeiikagiaemiraq0aasbaasqaaiabdmgapbqabagccqgg8bafcqqhmogrcqggsaaliigacqwf4oqcdawgaawcbagaemyaakgabeaakiabg2da9iabdqgaqjabcmcapiabdcfaqjabcicaoiab=h7axnaabaaaleaacqwgpbqaaeqaaogaeyypa0jaemoaaomaeiifawnaeuopdykaeiykakcaleaacqwgqbgacqgh9aqpcqaixaqmaeaacqwgubgba0gaeyyeiuoaasqaaiabdmgapjabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@5613@
please note that the edges in network Φ can either represent transcriptional regulation events or phosphorilation or post-translational effects, as we reconstruct the signal flow in the network based on the nested structure of the measured effects.
the effects on the e-genes that are measured are transcriptional effects, which are ultimately regulated by transcription factors.
some e-genes may be regulated by kinases, as due to the inherent nature of microarray measurements, it is impossible to distinguish between direct and indirect effects.
our approach
generalized inference framework
in their original work markowetz et al.
suppose the data matrix d to consist of counts, how often a specific e-gene shows an effect in ℓ experiment repetitions.
this requires a data discretization step, for which user specified type-i and type-ii error rates are assumed.
the choice of these parameters is certainly critical for the inference procedure, because it directly influences  and appears to be difficult to estimate.
markowetz et al.
suppose to have both, positive and negative controls  for this procedure, which for our data is not available .
in contrast, in our approach we make the assumption that d is an m × n matrix of  p-values, which specify the likelihood of e-gene i being differentially expressed after knock-down of s-gene k. the p-values are calculated using a method for detecting differential gene expression, e.g.
limma  <cit> .
this way various experimental designs, including dye swaps, on arbitrary chip platforms can be used in a simple manner.
we now suppose a decomposition of p as follows:
p=∏k=1np
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakgabeaakiabcyha8jabfa6agjabcycasggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiykakiaeyypa0zaaebcaeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcyha8jabfa6agjabcycasiab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiykakcaleaacqwgrbwacqgh9aqpcqaixaqmaeaacqwgubgba0gaey4diunaaaa@4f3a@
in accordance to  <cit>  this makes the assumption that knock-down experiments are statistically independent from each other.
hence, eq.
can be written down as
p=∏i=1m∑j=1n∏k=1npp
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebarcqgg8bafcqqhmogrcqggpaqkcqgh9aqpdaqewbqaamaaqahabawaaebcaeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcyha8jabfa6agjabcycasggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeyypa0jaemoaaomaeiykakiaemiuaalaeiikagiae8hude3aasbaasqaaiabdmgapbqabagccqgh9aqpcqwgqbgacqgg8bafcqqhmogrcqggpaqkasqaaiabdugarjabg2da9iabigdaxaqaaiabd6gaubqdcqghpis1aawcbagaemoaaomaeyypa0jaegymaedabagaemoba4ganiabgghildaaleaacqwgpbqacqgh9aqpcqaixaqmaeaacqwgtbqba0gaey4diunaaaa@6543@
the only thing missing is the definition of p. for this purpose we suppose the dik to be drawn from a mixture of a uniform  <cit>  distribution reflecting the null hypothesis and another distribution f <dig> reflecting the alternative hypothesis  <cit> :
p = γk + ·f <dig>  γk ∈
under the alternative hypothesis there is a high density for small p-values and a strong decrease for increasing p-values.
both distributions overlap with mixing coefficient γk.
p can therefore be decomposed as:
p={f1if Φ predicts an effect1otherwise
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcyha8jabfa6agjabcycasggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiykakiaeyypa0zaaiqabeaafaqaaegacaaabagaemozay2aasbaasqaaiabigdaxaqabagccqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcmcapaqaaiabbmgapjabbagamjabbccagiabfa6agjabbccagiabbchawjabbkhayjabbwgaljabbsgakjabbmgapjabbogajjabbsha0jabbohazjabbccagiabbggahjabb6gaujabbccagiabbwgaljabbagamjabbagamjabbwgaljabbogajjabbsha0bqaaiabigdaxaqaaiabb+gavjabbsha0jabbigaojabbwgaljabbkhayjabbeha3jabbmgapjabbohazjabbwgalbaaaiaawuhaaaaa@6f57@
the density function f <dig> reflects the strength of the knock-down effect on e-gene i under the alternative hypothesis.
if it is greater than  <dig> the alternative hypothesis would be accepted, and if it is smaller than  <dig> rejected.
in this work we assume f <dig> to be a mixture of a beta distribution  and a small uniform component:
f <dig> = πk + beta
in practice we set πk =  <dig>  and tuned the parameter βk on the full distribution of raw p-values for knock-down experiment k  such that f <dig> >  <dig>  if the benjamini-hochberg false discovery rate  <cit>  for dik was ≤ 10% and f <dig> ≤  <dig> otherwise.
an alternative treatment using a fitting procedure with expectation maximization  <cit>  is described in our recent publication  <cit> .
regularization
eq.
allows one to specify a prior p on the network structure itself.
this can be thought of as biasing the score of possible network hypotheses towards prior knowledge.
it is crucial to understand that in principle in any inference scheme there exist two competing goals: belief in prior assumptions/prior knowledge versus belief in data.
only trusting the data itself may lead to overfitting, whereas only trusting in prior assumptions does not give any new information and prevents learning.
therefore, we need a trade-off between both goals.
this technique is known as regularization in the machine learning literature  <cit> .
we have to take into account at this point that our assumptions may only be true up to a certain degree.
hence, for each edge we should suppose a prior probability reflecting the degree of belief in its existence.
in principle, this degree of belief can be very different for each edge.
we summarize all prior edge probabilities in an n × n matrix Φ^
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuqhmogrgaqcaaaa@2e36@.
making the assumption that all edge priors p are independent, i.e.
p=∏i,jp
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqqhmogrcqggpaqkcqgh9aqpdaqeqbqaaiabdcfaqjabcicaoiabfa6agnaabaaaleaacqwgpbqacqwgqbgaaeqaaogaeiykakcaleaacqwgpbqacqggsaalcqwgqbgaaeqaniabg+givdaaaa@3ef4@
allows us to define the connection between Φij and Φ^ij
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuqhmogrgaqcamaabaaaleaacqwgpbqacqwgqbgaaeqaaaaa@311a@ for each edge separately.
note that Φij ∈ { <dig> } depending on whether we set the edge i → j or not.
hence, for each edge we have a certain difference |Φij−Φ^ij|
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaabdaqaaiabfa6agnaabaaaleaacqwgpbqacqwgqbgaaeqaaogaeyoei0iafuopdykbakaadawgaawcbagaemyaakmaemoaaogabeaaaogaay5bslaawia7aaaa@399b@ to our prior assumptions.
the smaller this difference, the higher p should be.
we can therefore model p as a laplacian distribution with parameter λ :
p=λ2exp⁡
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqqhmogrdawgaawcbagaemyaakmaemoaaogabeaakiabcyha8hgaciab=t7asjabcmcapiabg2da9maalaaabagae83udwgabagaegomaidaaigbcwgaljabciha4jabcchawnaabmaabagaeyoei0iae83udw2aaqwaaeaacqqhmogrdawgaawcbagaemyaakmaemoaaogabeaakiabgkhitiqbfa6agzaajawaasbaasqaaiabdmgapjabdqgaqbqabaaakiaawea7cagliwoaaiaawicacaglpaaaaaa@5029@
if we now write down the log-posterior of eq.
log p ∝ log p + log p
∝log⁡p−λ∑i,j|Φij−Φ^ij|
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqghdistcyggsbabcqggvbwbcqggnbwzcqwgqbaucqggoaakcqwgebarcqgg8bafcqqhmogrcqggpaqkcqghsisliigacqwf7oabdaaeqbqaamaaemaabagaeuopdy0aasbaasqaaiabdmgapjabdqgaqbqabagccqghsislcuqhmogrgaqcamaabaaaleaacqwgpbqacqwgqbgaaeqaaagccaglhwuaayjcsdaaleaacqwgpbqacqggsaalcqwgqbgaaeqaniabgghildaaaa@4e82@
we see that λ specifies the trade-off between the model's fit to our data and our prior assumptions.
an important special case of the latter would be Φ^
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuqhmogrgaqcaaaa@2e36@ =  <dig>  i.e.
the matrix consisting only of zeros.
the meaning of this prior would be to generally prefer sparse networks structures over dense ones.
setting λ → ∞ corresponds to completely trusting the prior, whereas λ =  <dig> leads to a maximum likelihood estimate, i.e.
complete trust in data.
in practice we would like to choose a λ balancing both goals.
this leads to an instance of the classical model selection problem  in statistical learning.
one way of dealing with it is to tune λ such that the akaike information criterion
aic = - <dig> log p + 2d
becomes minimal  <cit> .
here d denotes the number of free parameters  in the network structure Φopt optimizing .
large scale network inference
the inference framework does not answer the question how to come up with a candidate network topology, which we would like to score.
markowetz et al.
<cit>  completely enumerate all possible topologies.
this is, however, only suitable for small networks of up to  <dig> s-genes.
for  <dig> s-genes there already exist more than  <dig> , <dig> and for  <dig> genes more than  <dig> possible network topologies.
in this context it should be noted that the scoring scheme cannot distinguish between two network hypotheses, if they only differ in transitive edges.
this issue is known as prediction equivalence.
hence, it only makes sense to consider the set of all transitively closed network hypotheses.
however, restricting ourselves to this limited class of network structures does not generally solve the problem, since even then the number of networks to consider scales in a similar way with the number of s-genes .
hence, we have to resort to heuristics.
stochastic sampling
a quite obvious idea to prevent the computational effort to enumerate all possible network hypothesis is to sample from the set of all transitively closed network graphs randomly.
we decided to use simulated annealing  here  <cit> .
sa is rather similar to markov chain monte carlo  sampling  <cit> , but additionally makes use of a so-called cooling scheme, which gradually decreases the neighborhood size of a given state in search space.
sa has been successfully applied to many difficult optimization problems from various disciplines, including bioinformatics  <cit> .
in order to use sa, we have to define a state transition function t : s → s, which defines how to come from one graph to a modified one in search space.
a special challenge in this context is that we need to guarantee that in principle all possible transitively closed network topologies can be reached by our function t.
supposed we have functions add and del, which add and remove edges from a given transitively closed graph and produce a new one from this.
we will restrict ourselves to the set of all transitively closed directed graphs  here for reasons that will become clear soon.
we now define add and del in a formal way as follows:
definition 1
let tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ be the set of all transitively closed graphs dags with n nodes.
suppose the nodes for g ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ are indexed in some way by natural numbers.
we define add : tn×ℰi→tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqwgpbqaaeqaaogaeyokh4qae83ext1aasbaasqaaiabd6gaubqabaaaaa@43ff@ as a function, which takes a graph g and a pair of node indices  ∈ ℰ1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegymaedabeaaaaa@38a0@ ⊂ ℕ × ℕ, inserts the edge between i, j and transitively closes g then.
ℰ1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegymaedabeaaaaa@38a0@ is defined as the set of all pairs of node indices, for which there exists no edge in g and which after edge insertion do not violate the dag property.
likewise, del : tn×ℰ2→tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqaiyagmaeqaaogaeyokh4qae83ext1aasbaasqaaiabd6gaubqabaaaaa@4396@ is a function, which takes a graph g and a pair of node indices  ∈ ℰ2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegomaidabeaaaaa@38a2@ ⊂ ℕ × ℕ. g is transitively reduced, the edge between i, j deleted and then transitively closed again.
with a transitive reduction g* of a graph g ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ we mean a graph with a minimal number of edges such that the transitive closure of g* is g. ℰ2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegomaidabeaaaaa@38a2@ is defined as the set of all pairs of node indices, for which there exists an edge in g.
in contrast to general graphs the transitive reduction of a dag is unique  <cit> , which is the reason for our restriction.
this way we can guarantee that the del function is well defined and injective.
this gives rise to the following lemma:
lemma 2
the operations add : tn×ℰ1→t′n
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqaixaqmaeqaaogaeyokh4qaf83extlbauaadawgaawcbagaemoba4gabeaaaaa@43a0@and del : tn×ℰ2→t′n
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqaiyagmaeqaaogaeyokh4qaf83extlbauaadawgaawcbagaemoba4gabeaaaaa@43a2@have the following inverse property to each other: for any g ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@del, i, j) = g, provided that edge i, j does not exist in g and does not violate the dag property.
likewise, add, i, j) = g, provided that edge i, j exists in g.
proof.
the transitive reduction of a transitively closed dag is unique  <cit> .
hence, the del operation is a well defined injective function.
additionally note that in the add operation we can never insert an edge, which lies in the transitive hull of g ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ , since otherwise it would have been there already .
therefore, for any g ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@del, i, j) = g, provided that edge i, j does not exist in g and does not violate the dag property.
likewise, add, i, j) = g, provided that edge i, j exists in g.   □
theorem 3
from any graph a ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@in the set of all transitively closed dags we can reach any other graph b ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@using add and del operations.
proof.
let  <dig> ∈ tn
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ be the dag without any edges .
let us denote by x → y that from dag x we can get dag y using add and del operations on appropriate edges.
note that we have a →  <dig> and b →  <dig>  because we can use del to remove edges successively until no edges are left.
because add and del are inverse operations to each other according to the lemma, we have  <dig> → a,  <dig> → b. hence, we get a →  <dig> → b, which proves the theorem.
□
still, the sa approach suffers from a potential problem: both, the add and the del operation, at the bottom line perform a whole cascade of changes on the original graph.
thus there may be harsh changes in the scoring function when applying such an operation to a given candidate network.
this may make it difficult to come close to the optimal network hypothesis.
module networks
rather than looking for a complete network hypothesis at once the idea of the module network is to build up a graph from smaller subgraphs, called modules in the following.
the module network is thus a divide and conquer approach: we first split the complete node set into smaller subgroups.
this can be done by pam clustering  <cit>  on the p-value density profiles of the s-genes.
the idea is that s-genes with a similar e-gene response profile  should be close in the signaling path.
the number of clusters for the pam clustering is chosen between  <dig> and half of the number of s-genes such that the average silhouette index becomes maximal.
the silhouette value for each point in a cluster is a measure of how similar that point is to points in its own cluster vs. points in other clusters, and ranges from - <dig> to + <dig>  <cit> .
it is defined as:
s=min⁡j)−d¯wmax⁡,min⁡j)
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgtbwucqggoaakcqwgpbqacqggpaqkcqgh9aqpdawcaaqaaigbc2gatjabcmgapjabc6gaunaabaaaleaacqwgqbgaaeqaaogaeiikagiafmizaqmbaebadawgaawcbagaemoqaieabeaakiabcicaoiabdmgapjabcycasiabdqgaqjabcmcapiabcmcapiabgkhitiqbdsgakzaarawaasbaasqaaiabdefaxbqabagccqggoaakcqwgpbqacqggpaqkaeaacyggtbqbcqgghbqycqgg4baecqggoaakcuwgkbazgaqeamaabaaaleaacqwgxbwvaeqaaogaeiikagiaemyaakmaeiykakiaeiilawiagiyba0maeiyaakmaeioba42aasbaasqaaiabdqgaqbqabagccqggoaakcuwgkbazgaqeamaabaaaleaacqwgcbgqaeqaaogaeiikagiaemyaakmaeiilawiaemoaaomaeiykakiaeiykakcaaaaa@6349@
where d¯w
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgkbazgaqeamaabaaaleaacqwgxbwvaeqaaogaeiikagiaemyaakmaeiykakcaaa@328f@ is the average distance from the i-th point to the other points in its own cluster, and d¯b
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgkbazgaqeamaabaaaleaacqwgcbgqaeqaaogaeiikagiaemyaakmaeiilawiaemoaaomaeiykakcaaa@34a2@ is the average distance from the i-th point to points in another cluster j.
each cluster of s-genes now forms one module.
these modules are eventually further subdivided into smaller submodules until each submodule contains only  <dig> s-genes at most.
this way we obtain a tree structure of modules, where each node  has children .
we begin with estimating the leaves in the module tree.
as each leaf module can contain  <dig> s-genes at maximum this can be done by enumerating all possible transitively closed network hypotheses and taking the highest scoring one.
after the leaves in the module tree have been built, their connection is estimated.
for this purpose we score all pairwise connections between any pair of s-genes from leaves l <dig> and l <dig>  denoting by |l1| and |l2| the number of s-genes in both leaves, these are  <dig>  |l1|·|l2| tests altogether, because between any pair of s-genes  we can either have no edge, an edge from n <dig> to n <dig>  an edge from n <dig> to n <dig> or an edge in both directions.
after the best connection between l <dig> and l <dig> has been estimated, the corresponding subgraph is transitively closed.
after all connections between leaves belonging to the same submodule in the module tree have been established, we recursively continue with connecting submodules in the same fashion as we did for leaf modules until the topology for the total network is completed.
generalized inference framework: proof of principle
to show the correctness of our generalized inference framework, we conducted experiments on the drosophila dataset by boutros et al.
.
this dataset was also employed by markowetz et al.
<cit>  as a proof of principle with discretized data.
the dataset consists of expression profiles from  <dig> affymetrix microarrays:  <dig> genes  were stimulated by lipopolysaccharide  for  <dig> minutes and then knocked-down by rnai with  <dig> replicates for each expression profile.
additionally there were  <dig> replicates of control experiments without lps and rnai and  <dig> replicates of expression profiling with lps but without rnai.
the dataset is available in a preprocessed form as a supplement of  <cit> .
we took the same  <dig> genes showing a secondary effect  as used in this publication and calculated p-values for differential gene expression between lps stimulated and knock-down conditions by fitting an empirical bayes model using the limma package in the r statistical computing environment  <cit> .
we enumerated all possible  <dig> transitively closed network topologies and calculated their scores using .
the scores of the top  <dig> models and the best model are depicted in figure  <dig>  the score distribution of the  <dig> top models is slightly different, because of our modified inference scheme.
we had a closer look at the best  <dig> models and found them to be identical to those shown in  <cit>  .
the second best model differs from the best model only in the missing edge key → rel.
the next two models are either missing the edge tak → rel or tak → key.
the key feature is preserved in all of them: the signal runs through tak before splitting into two pathway branches, one containing mkk4/hep, the other both key and rel.
this fits exactly to the findings of boutros et al.
<cit> .
large scale inference: evaluation on artificial networks
to test our methods and to get better insights into the performance of our large scale inference methods, we generated data from artificial random networks.
the sampling procedure for artificial networks is described in section "methods".
we sampled networks with n =  <dig>   <dig>  s-genes.
for each number of s-genes we varied the number m =  <dig>   <dig> ..., 4n of e-genes and the parameter β =  <dig>   <dig>  describing the beta component of the f <dig> distribution .
we compared the sa approach with the module network.
we evaluated both methods in terms of average sensitivity  and specificity  over  <dig> generated networks for each parameter combination .
the initial temperature for the sa was set to  <dig> and the maximum number of iterations to 100n.
the initial network structure was always the graph with no edges.
a logarithmic temperature cooling scheme according to  <cit>  was used.
the results are shown in figure  <dig> – figure  <dig>  in general all methods achieve a higher specificity than sensitivity, which is due to our "p-value" sampling strategy, and they show a high robustness against a varying number of e-genes.
all in all the module network approach shows a superiority to the sa approach, especially in terms of sensitivity.
using module networks the sensitivity and specificity for n =  <dig> goes up to almost 100%.
for n =  <dig>  the sensitivity lies around 80%, while the specificity reaches more than 90%.
moreover, for all tested values of β the curves are relatively close together.
we also compared the computation times for both approaches and found the module network to be substantially faster for n =  <dig>  .
the average running time for network inference with n =  <dig> nodes was only 4s with the module network on our amd dual core opteron  <dig> ghz machine.
in conclusion we think that the module network offers the most reliable and fast mechanism for large scale network inference among our tested approaches and is therefore taken as our inference method in the following section.
application to rnai data from human er-α pathway
we applied the module network to infer the complete topology for a network of  <dig> silenced genes  in the er-α pathway.
the  <dig> genes were selected from previous microarray studies in our department to be influenced by er status in breast cancer patients.
each of the  <dig> genes was silenced individually using two different sirnas, and the effect on gene expression was studied on whole genome cdna microarrays.
the data were generated in our department, details on the data generation and preprocessing steps are described in section "methods".
we found several known interdependencies between e- and s-genes as well as among s-genes by an intensive literature screen.
the corresponding information was obtained from the ingenuity™ software and can be visualized in form of a interdependency graph .
it represents some prior knowledge, which can be used for the network inference with our module networks algorithm via the regularization technique .
we considered  <dig> situations for the network inference:  <dig>  no prior knowledge ,  <dig>  inclusion of known interdependencies between e- and s-genes : for known interdependencies we set p =  <dig> and p =  <dig>  while otherwise we have a uniform prior p = 1n
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabigdaxaqaaiabd6gaubaaaaa@2f11@.
in case of several interdependencies for one e-gene p is rescaled appropriately.
<dig>  additional inclusion of literature knowledge for interdependencies between s-genes with parameter λ chosen from the set { <dig> ,..., 10-2} via the aic criterion .
the prior for the network structure  is chosen such that  =  <dig>  if an interaction between s-gene i and s-gene j is known and  =  <dig>  otherwise.
that means missing a known interdependency is punished more than introducing an edge where nothing is known.
to ensure the statistical stability of the inferred network we employed bootstrapping: we sampled m e-genes from the total set of e-genes  <dig> times with replacement and each time ran the module network for topology induction.
at the end we only considered edges, which were found in more than 50% of all bootstrap trials.
CONCLUSIONS
we proposed a method for reconstructing signaling pathways from secondary effects, which were observed on microarray after silencing genes of interest via rnai.
our approach systematically extends and generalizes previous work by markowetz et al: an inference scheme was developed, which can deal with p-values for differential gene expression and does not rely on discretized data only.
regularization was employed to incorporate prior assumptions on the network architecture into our framework.
finally, new algorithms for large scale inference of signaling pathways were developed and evaluated in a systematic fashion on artificially created data.
thereby, our module network, which recursively build up the complete topology from smaller pieces, revealed the best performance in terms of sensitivity and specificity.
we used the module network to infer the signaling pathway for  <dig> genes in the er-α pathway in human mcf- <dig> breast cancer cells and used a bootstrapping approach to ensure the statistical stability of the result.
the induced edges in our inferred network were found with good consistency and could in many cases be also confirmed by the literature.
future biological experiments are planned to validate our reconstructed network in a systematic way.
in conclusion of our results we think that our approach offers a scale-able, reliable and fairly general way for large scale inference of signaling pathways from secondary effects and therefore provides researchers with a valuable tool to gain insight into complex cellular processes.
the code for the module network inference method is available in the latest version of the -r-package nem, which can be obtained from the bioconductor homepage .

 in the modern field of systems biology scientists aim to get insights into the architecture and behavior of complex cellular and genomic processes. an important task in this context is the detection of novel interdependencies between gene products. this insight into the genomic networks is an important step towards a better understanding of the functional aspects of a biological system and of great value for drug target identification at a later stage. within this context modern dna microarray technology plays an important role. in addition, the advent of rna silencing techniques has further increased its power by allowing the selective knock-down of certain genes of interest. this may enable us to detect interdependencies between gene products on a non-transcriptional level. the genes of interest are knocked down individually, and the respective downstream effects on gene expression are measured by using genome-wide microarrays. by observing the nested structure of significant up or down regulations of affected genes, this may allow one to reverse engineer features of the upstream signaling pathway. in a recent work markowetz et al. introduced a method to reverse engineer the signaling pathway between perturbed genes using the nested structure of secondary downstream effects. they developed a bayesian statistical framework, in which for a given network hypothesis one can calculate a score and thus can reduce the set of all possible networks to the most likely ones. a severe limitation of this method lies, however, in the restriction to small networks of up to genes, because the method completely enumerates all possible network hypotheses. furthermore, a difficulty in the practical use is the required binary discretization of the data. in our work we therefore aim to extend the framework by markowetz et al. in order to make it practically applicable for a broader range of real life problems. we are thereby motivated by biological experiments conducted in our department: genes in the er-α pathway in human mcf- breast cancer cells were silenced via small interfering rnas and the effects on gene expression were subsequently measured on cdna microarrays. our extensions of the original approach go in three directions: first, we introduce a way to omit the data discretization step needed via a calculation based on p-values instead, which is more suitable for our data and makes the whole framework more flexible. second, we show how prior assumptions on the network structure can be incorporated into the network scoring scheme via techniques from regularization theory. third and most important, we develop and investigate methods to scale up the network inference to large scale networks. for this purpose two approaches are considered: simulated annealing on a restricted set of possible networks and our so-called module networks, which build the complete network recursively from smaller pieces that are connected subsequently. in order to validate these approaches we conduct studies on artificially created networks and show that module networks offer the highest sensitivity and specificity in the reconstruction of edges in the networks. finally, we demonstrate the applicability of our approach to real data by inferring the complete genes er-α signaling pathway network. using a bootstrapping approach this reconstruction can be found with good statistical stability and hence seems to be reliable.  statistical inference framework for signaling pathways from rnai data we start with a brief review of the statistical inference framework for signaling pathways by markowetz et al.: in general within this framework one distinguishes between silenced genes and genes showing a downstream effect. it is assumed that each e-gene is attached to a single s-gene only. knocking down a specific s-gene sk interrupts signal flow in the downstream pathway, and hence an effect on the e-genes attached to sk and all s-genes depending on sk is expected. let us assume n knock-downs are performed and there exist m e-genes in total. the outcomes of these experiments are summarized in an m × n data matrix d. according to bayes' formula a specific network hypothesis Φ ∈ { }n × { }n can be scored as: p ∝ pp the position of the e-genes is introduced as a model parameter Θ = {θi|θi ∈ {..., n}, i =..., m}, i.e θi = j, if e-gene i is attached to s-gene j. assuming independence of the observations di in the data matrix d one can write down the conditional likelihood p|Φ, Θ) as: p=∏i=1mp mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebarcqgg8bafcqqhmogrcqggsaalcqqhyoqucqggpaqkcqgh9aqpdaqewbqaaiabdcfaqjabcicaoiabdseaenaabaaaleaacqwgpbqaaeqaaogaeiifawnaeuopdykaeiilawcccigae8hude3aasbaasqaaiabdmgapbqabagccqggpaqkasqaaiabdmgapjabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@4a79@ it is furthermore assumed that all parameters θi are statistically independent, i.e. p=∏i=1mp mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqqhyoqucqgg8bafcqqhmogrcqggpaqkcqgh9aqpdaqewbqaaiabdcfaqjabcicaoggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiifawnaeuopdykaeiykakcaleaacqwgpbqacqgh9aqpcqaixaqmaeaacqwgtbqba0gaey4diunaaaa@4506@ the likelihood p can now be written as: p = ∫p pdΘ =∏i=1m∑j=1npp mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqgh9aqpdaqewbqaamaaqahabagaemiuaalaeiikagiaemiraq0aasbaasqaaiabdmgapbqabagccqgg8bafcqqhmogrcqggsaaliigacqwf4oqcdawgaawcbagaemyaakgabeaakiabg2da9iabdqgaqjabcmcapiabdcfaqjabcicaoiab=h7axnaabaaaleaacqwgpbqaaeqaaogaeyypa0jaemoaaomaeiifawnaeuopdykaeiykakcaleaacqwgqbgacqgh9aqpcqaixaqmaeaacqwgubgba0gaeyyeiuoaasqaaiabdmgapjabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@5613@ please note that the edges in network Φ can either represent transcriptional regulation events or phosphorilation or post-translational effects, as we reconstruct the signal flow in the network based on the nested structure of the measured effects. the effects on the e-genes that are measured are transcriptional effects, which are ultimately regulated by transcription factors. some e-genes may be regulated by kinases, as due to the inherent nature of microarray measurements, it is impossible to distinguish between direct and indirect effects. our approach generalized inference framework in their original work markowetz et al. suppose the data matrix d to consist of counts, how often a specific e-gene shows an effect in ℓ experiment repetitions. this requires a data discretization step, for which user specified type-i and type-ii error rates are assumed. the choice of these parameters is certainly critical for the inference procedure, because it directly influences and appears to be difficult to estimate. markowetz et al. suppose to have both, positive and negative controls for this procedure, which for our data is not available. in contrast, in our approach we make the assumption that d is an m × n matrix of p-values, which specify the likelihood of e-gene i being differentially expressed after knock-down of s-gene k. the p-values are calculated using a method for detecting differential gene expression, e.g. limma. this way various experimental designs, including dye swaps, on arbitrary chip platforms can be used in a simple manner. we now suppose a decomposition of p as follows: p=∏k=1np mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakgabeaakiabcyha8jabfa6agjabcycasggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiykakiaeyypa0zaaebcaeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcyha8jabfa6agjabcycasiab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiykakcaleaacqwgrbwacqgh9aqpcqaixaqmaeaacqwgubgba0gaey4diunaaaa@4f3a@ in accordance to this makes the assumption that knock-down experiments are statistically independent from each other. hence, eq. can be written down as p=∏i=1m∑j=1n∏k=1npp mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebarcqgg8bafcqqhmogrcqggpaqkcqgh9aqpdaqewbqaamaaqahabawaaebcaeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcyha8jabfa6agjabcycasggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeyypa0jaemoaaomaeiykakiaemiuaalaeiikagiae8hude3aasbaasqaaiabdmgapbqabagccqgh9aqpcqwgqbgacqgg8bafcqqhmogrcqggpaqkasqaaiabdugarjabg2da9iabigdaxaqaaiabd6gaubqdcqghpis1aawcbagaemoaaomaeyypa0jaegymaedabagaemoba4ganiabgghildaaleaacqwgpbqacqgh9aqpcqaixaqmaeaacqwgtbqba0gaey4diunaaaa@6543@ the only thing missing is the definition of p. for this purpose we suppose the dik to be drawn from a mixture of a uniform distribution reflecting the null hypothesis and another distribution f reflecting the alternative hypothesis: p = γk + ·f γk ∈ under the alternative hypothesis there is a high density for small p-values and a strong decrease for increasing p-values. both distributions overlap with mixing coefficient γk. p can therefore be decomposed as: p={f1if Φ predicts an effect1otherwise mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcyha8jabfa6agjabcycasggaciab=h7axnaabaaaleaacqwgpbqaaeqaaogaeiykakiaeyypa0zaaiqabeaafaqaaegacaaabagaemozay2aasbaasqaaiabigdaxaqabagccqggoaakcqwgebardawgaawcbagaemyaakmaem4aasgabeaakiabcmcapaqaaiabbmgapjabbagamjabbccagiabfa6agjabbccagiabbchawjabbkhayjabbwgaljabbsgakjabbmgapjabbogajjabbsha0jabbohazjabbccagiabbggahjabb6gaujabbccagiabbwgaljabbagamjabbagamjabbwgaljabbogajjabbsha0bqaaiabigdaxaqaaiabb+gavjabbsha0jabbigaojabbwgaljabbkhayjabbeha3jabbmgapjabbohazjabbwgalbaaaiaawuhaaaaa@6f57@ the density function f reflects the strength of the knock-down effect on e-gene i under the alternative hypothesis. if it is greater than the alternative hypothesis would be accepted, and if it is smaller than rejected. in this work we assume f to be a mixture of a beta distribution and a small uniform component: f = πk + beta in practice we set πk = and tuned the parameter βk on the full distribution of raw p-values for knock-down experiment k such that f > if the benjamini-hochberg false discovery rate for dik was ≤ 10% and f ≤ otherwise. an alternative treatment using a fitting procedure with expectation maximization is described in our recent publication. regularization eq. allows one to specify a prior p on the network structure itself. this can be thought of as biasing the score of possible network hypotheses towards prior knowledge. it is crucial to understand that in principle in any inference scheme there exist two competing goals: belief in prior assumptions/prior knowledge versus belief in data. only trusting the data itself may lead to overfitting, whereas only trusting in prior assumptions does not give any new information and prevents learning. therefore, we need a trade-off between both goals. this technique is known as regularization in the machine learning literature. we have to take into account at this point that our assumptions may only be true up to a certain degree. hence, for each edge we should suppose a prior probability reflecting the degree of belief in its existence. in principle, this degree of belief can be very different for each edge. we summarize all prior edge probabilities in an n × n matrix Φ^ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuqhmogrgaqcaaaa@2e36@. making the assumption that all edge priors p are independent, i.e. p=∏i,jp mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqqhmogrcqggpaqkcqgh9aqpdaqeqbqaaiabdcfaqjabcicaoiabfa6agnaabaaaleaacqwgpbqacqwgqbgaaeqaaogaeiykakcaleaacqwgpbqacqggsaalcqwgqbgaaeqaniabg+givdaaaa@3ef4@ allows us to define the connection between Φij and Φ^ij mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuqhmogrgaqcamaabaaaleaacqwgpbqacqwgqbgaaeqaaaaa@311a@ for each edge separately. note that Φij ∈ { } depending on whether we set the edge i → j or not. hence, for each edge we have a certain difference |Φij−Φ^ij| mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaabdaqaaiabfa6agnaabaaaleaacqwgpbqacqwgqbgaaeqaaogaeyoei0iafuopdykbakaadawgaawcbagaemyaakmaemoaaogabeaaaogaay5bslaawia7aaaa@399b@ to our prior assumptions. the smaller this difference, the higher p should be. we can therefore model p as a laplacian distribution with parameter λ: p=λ2exp⁡ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgqbaucqggoaakcqqhmogrdawgaawcbagaemyaakmaemoaaogabeaakiabcyha8hgaciab=t7asjabcmcapiabg2da9maalaaabagae83udwgabagaegomaidaaigbcwgaljabciha4jabcchawnaabmaabagaeyoei0iae83udw2aaqwaaeaacqqhmogrdawgaawcbagaemyaakmaemoaaogabeaakiabgkhitiqbfa6agzaajawaasbaasqaaiabdmgapjabdqgaqbqabaaakiaawea7cagliwoaaiaawicacaglpaaaaaa@5029@ if we now write down the log-posterior of eq. log p ∝ log p + log p ∝log⁡p−λ∑i,j|Φij−Φ^ij| mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqghdistcyggsbabcqggvbwbcqggnbwzcqwgqbaucqggoaakcqwgebarcqgg8bafcqqhmogrcqggpaqkcqghsisliigacqwf7oabdaaeqbqaamaaemaabagaeuopdy0aasbaasqaaiabdmgapjabdqgaqbqabagccqghsislcuqhmogrgaqcamaabaaaleaacqwgpbqacqwgqbgaaeqaaagccaglhwuaayjcsdaaleaacqwgpbqacqggsaalcqwgqbgaaeqaniabgghildaaaa@4e82@ we see that λ specifies the trade-off between the model's fit to our data and our prior assumptions. an important special case of the latter would be Φ^ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuqhmogrgaqcaaaa@2e36@ = i.e. the matrix consisting only of zeros. the meaning of this prior would be to generally prefer sparse networks structures over dense ones. setting λ → ∞ corresponds to completely trusting the prior, whereas λ = leads to a maximum likelihood estimate, i.e. complete trust in data. in practice we would like to choose a λ balancing both goals. this leads to an instance of the classical model selection problem in statistical learning. one way of dealing with it is to tune λ such that the akaike information criterion aic = - log p + 2d becomes minimal. here d denotes the number of free parameters in the network structure Φopt optimizing. large scale network inference the inference framework does not answer the question how to come up with a candidate network topology, which we would like to score. markowetz et al. completely enumerate all possible topologies. this is, however, only suitable for small networks of up to s-genes. for s-genes there already exist more than, and for genes more than possible network topologies. in this context it should be noted that the scoring scheme cannot distinguish between two network hypotheses, if they only differ in transitive edges. this issue is known as prediction equivalence. hence, it only makes sense to consider the set of all transitively closed network hypotheses. however, restricting ourselves to this limited class of network structures does not generally solve the problem, since even then the number of networks to consider scales in a similar way with the number of s-genes. hence, we have to resort to heuristics. stochastic sampling a quite obvious idea to prevent the computational effort to enumerate all possible network hypothesis is to sample from the set of all transitively closed network graphs randomly. we decided to use simulated annealing here. sa is rather similar to markov chain monte carlo sampling, but additionally makes use of a so-called cooling scheme, which gradually decreases the neighborhood size of a given state in search space. sa has been successfully applied to many difficult optimization problems from various disciplines, including bioinformatics. in order to use sa, we have to define a state transition function t: s → s, which defines how to come from one graph to a modified one in search space. a special challenge in this context is that we need to guarantee that in principle all possible transitively closed network topologies can be reached by our function t. supposed we have functions add and del, which add and remove edges from a given transitively closed graph and produce a new one from this. we will restrict ourselves to the set of all transitively closed directed graphs here for reasons that will become clear soon. we now define add and del in a formal way as follows: definition 1 let tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ be the set of all transitively closed graphs dags with n nodes. suppose the nodes for g ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ are indexed in some way by natural numbers. we define add: tn×ℰi→tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqwgpbqaaeqaaogaeyokh4qae83ext1aasbaasqaaiabd6gaubqabaaaaa@43ff@ as a function, which takes a graph g and a pair of node indices ∈ ℰ1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegymaedabeaaaaa@38a0@ ⊂ ℕ × ℕ, inserts the edge between i, j and transitively closes g then. ℰ1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegymaedabeaaaaa@38a0@ is defined as the set of all pairs of node indices, for which there exists no edge in g and which after edge insertion do not violate the dag property. likewise, del: tn×ℰ2→tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqaiyagmaeqaaogaeyokh4qae83ext1aasbaasqaaiabd6gaubqabaaaaa@4396@ is a function, which takes a graph g and a pair of node indices ∈ ℰ2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegomaidabeaaaaa@38a2@ ⊂ ℕ × ℕ. g is transitively reduced, the edge between i, j deleted and then transitively closed again. with a transitive reduction g* of a graph g ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ we mean a graph with a minimal number of edges such that the transitive closure of g* is g. ℰ2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwfwesrdawgaawcbagaegomaidabeaaaaa@38a2@ is defined as the set of all pairs of node indices, for which there exists an edge in g. in contrast to general graphs the transitive reduction of a dag is unique, which is the reason for our restriction. this way we can guarantee that the del function is well defined and injective. this gives rise to the following lemma: lemma 2 the operations add: tn×ℰ1→t′n mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqaixaqmaeqaaogaeyokh4qaf83extlbauaadawgaawcbagaemoba4gabeaaaaa@43a0@and del: tn×ℰ2→t′n mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaakiabgena0kab=btifnaabaaaleaacqaiyagmaeqaaogaeyokh4qaf83extlbauaadawgaawcbagaemoba4gabeaaaaa@43a2@have the following inverse property to each other: for any g ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@del, i, j) = g, provided that edge i, j does not exist in g and does not violate the dag property. likewise, add, i, j) = g, provided that edge i, j exists in g. proof. the transitive reduction of a transitively closed dag is unique. hence, the del operation is a well defined injective function. additionally note that in the add operation we can never insert an edge, which lies in the transitive hull of g ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@, since otherwise it would have been there already. therefore, for any g ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@del, i, j) = g, provided that edge i, j does not exist in g and does not violate the dag property. likewise, add, i, j) = g, provided that edge i, j exists in g.   □ theorem 3 from any graph a ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@in the set of all transitively closed dags we can reach any other graph b ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@using add and del operations. proof. let ∈ tn mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaat0uy0hwztfgdpnwy1egaryqthrhal1wy0l2yhvdaiqaacqwftepvdawgaawcbagaemoba4gabeaaaaa@39d7@ be the dag without any edges. let us denote by x → y that from dag x we can get dag y using add and del operations on appropriate edges. note that we have a → and b → because we can use del to remove edges successively until no edges are left. because add and del are inverse operations to each other according to the lemma, we have → a, → b. hence, we get a → → b, which proves the theorem.   □ still, the sa approach suffers from a potential problem: both, the add and the del operation, at the bottom line perform a whole cascade of changes on the original graph. thus there may be harsh changes in the scoring function when applying such an operation to a given candidate network. this may make it difficult to come close to the optimal network hypothesis. module networks rather than looking for a complete network hypothesis at once the idea of the module network is to build up a graph from smaller subgraphs, called modules in the following. the module network is thus a divide and conquer approach: we first split the complete node set into smaller subgroups. this can be done by pam clustering on the p-value density profiles of the s-genes. the idea is that s-genes with a similar e-gene response profile should be close in the signaling path. the number of clusters for the pam clustering is chosen between and half of the number of s-genes such that the average silhouette index becomes maximal. the silhouette value for each point in a cluster is a measure of how similar that point is to points in its own cluster vs. points in other clusters, and ranges from - to +. it is defined as: s=min⁡j)−d¯wmax⁡,min⁡j) mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgtbwucqggoaakcqwgpbqacqggpaqkcqgh9aqpdawcaaqaaigbc2gatjabcmgapjabc6gaunaabaaaleaacqwgqbgaaeqaaogaeiikagiafmizaqmbaebadawgaawcbagaemoqaieabeaakiabcicaoiabdmgapjabcycasiabdqgaqjabcmcapiabcmcapiabgkhitiqbdsgakzaarawaasbaasqaaiabdefaxbqabagccqggoaakcqwgpbqacqggpaqkaeaacyggtbqbcqgghbqycqgg4baecqggoaakcuwgkbazgaqeamaabaaaleaacqwgxbwvaeqaaogaeiikagiaemyaakmaeiykakiaeiilawiagiyba0maeiyaakmaeioba42aasbaasqaaiabdqgaqbqabagccqggoaakcuwgkbazgaqeamaabaaaleaacqwgcbgqaeqaaogaeiikagiaemyaakmaeiilawiaemoaaomaeiykakiaeiykakcaaaaa@6349@ where d¯w mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgkbazgaqeamaabaaaleaacqwgxbwvaeqaaogaeiikagiaemyaakmaeiykakcaaa@328f@ is the average distance from the i-th point to the other points in its own cluster, and d¯b mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacuwgkbazgaqeamaabaaaleaacqwgcbgqaeqaaogaeiikagiaemyaakmaeiilawiaemoaaomaeiykakcaaa@34a2@ is the average distance from the i-th point to points in another cluster j. each cluster of s-genes now forms one module. these modules are eventually further subdivided into smaller submodules until each submodule contains only s-genes at most. this way we obtain a tree structure of modules, where each node has children. we begin with estimating the leaves in the module tree. as each leaf module can contain s-genes at maximum this can be done by enumerating all possible transitively closed network hypotheses and taking the highest scoring one. after the leaves in the module tree have been built, their connection is estimated. for this purpose we score all pairwise connections between any pair of s-genes from leaves l and l denoting by |l1| and |l2| the number of s-genes in both leaves, these are |l1|·|l2| tests altogether, because between any pair of s-genes we can either have no edge, an edge from n to n an edge from n to n or an edge in both directions. after the best connection between l and l has been estimated, the corresponding subgraph is transitively closed. after all connections between leaves belonging to the same submodule in the module tree have been established, we recursively continue with connecting submodules in the same fashion as we did for leaf modules until the topology for the total network is completed. generalized inference framework: proof of principle to show the correctness of our generalized inference framework, we conducted experiments on the drosophila dataset by boutros et al.. this dataset was also employed by markowetz et al. as a proof of principle with discretized data. the dataset consists of expression profiles from affymetrix microarrays: genes were stimulated by lipopolysaccharide for minutes and then knocked-down by rnai with replicates for each expression profile. additionally there were replicates of control experiments without lps and rnai and replicates of expression profiling with lps but without rnai. the dataset is available in a preprocessed form as a supplement of. we took the same genes showing a secondary effect as used in this publication and calculated p-values for differential gene expression between lps stimulated and knock-down conditions by fitting an empirical bayes model using the limma package in the r statistical computing environment. we enumerated all possible transitively closed network topologies and calculated their scores using. the scores of the top models and the best model are depicted in figure the score distribution of the top models is slightly different, because of our modified inference scheme. we had a closer look at the best models and found them to be identical to those shown in. the second best model differs from the best model only in the missing edge key → rel. the next two models are either missing the edge tak → rel or tak → key. the key feature is preserved in all of them: the signal runs through tak before splitting into two pathway branches, one containing mkk4/hep, the other both key and rel. this fits exactly to the findings of boutros et al.. large scale inference: evaluation on artificial networks to test our methods and to get better insights into the performance of our large scale inference methods, we generated data from artificial random networks. the sampling procedure for artificial networks is described in section "methods". we sampled networks with n = s-genes. for each number of s-genes we varied the number m =..., 4n of e-genes and the parameter β = describing the beta component of the f distribution. we compared the sa approach with the module network. we evaluated both methods in terms of average sensitivity and specificity over generated networks for each parameter combination. the initial temperature for the sa was set to and the maximum number of iterations to 100n. the initial network structure was always the graph with no edges. a logarithmic temperature cooling scheme according to was used. the results are shown in figure – figure in general all methods achieve a higher specificity than sensitivity, which is due to our "p-value" sampling strategy, and they show a high robustness against a varying number of e-genes. all in all the module network approach shows a superiority to the sa approach, especially in terms of sensitivity. using module networks the sensitivity and specificity for n = goes up to almost 100%. for n = the sensitivity lies around 80%, while the specificity reaches more than 90%. moreover, for all tested values of β the curves are relatively close together. we also compared the computation times for both approaches and found the module network to be substantially faster for n =. the average running time for network inference with n = nodes was only 4s with the module network on our amd dual core opteron ghz machine. in conclusion we think that the module network offers the most reliable and fast mechanism for large scale network inference among our tested approaches and is therefore taken as our inference method in the following section. application to rnai data from human er-α pathway we applied the module network to infer the complete topology for a network of silenced genes in the er-α pathway. the genes were selected from previous microarray studies in our department to be influenced by er status in breast cancer patients. each of the genes was silenced individually using two different sirnas, and the effect on gene expression was studied on whole genome cdna microarrays. the data were generated in our department, details on the data generation and preprocessing steps are described in section "methods". we found several known interdependencies between e- and s-genes as well as among s-genes by an intensive literature screen. the corresponding information was obtained from the ingenuity™ software and can be visualized in form of a interdependency graph. it represents some prior knowledge, which can be used for the network inference with our module networks algorithm via the regularization technique. we considered situations for the network inference: no prior knowledge, inclusion of known interdependencies between e- and s-genes: for known interdependencies we set p = and p = while otherwise we have a uniform prior p = 1n mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabigdaxaqaaiabd6gaubaaaaa@2f11@. in case of several interdependencies for one e-gene p is rescaled appropriately. additional inclusion of literature knowledge for interdependencies between s-genes with parameter λ chosen from the set {,..., 10-2} via the aic criterion. the prior for the network structure is chosen such that = if an interaction between s-gene i and s-gene j is known and = otherwise. that means missing a known interdependency is punished more than introducing an edge where nothing is known. to ensure the statistical stability of the inferred network we employed bootstrapping: we sampled m e-genes from the total set of e-genes times with replacement and each time ran the module network for topology induction. at the end we only considered edges, which were found in more than 50% of all bootstrap trials.  we proposed a method for reconstructing signaling pathways from secondary effects, which were observed on microarray after silencing genes of interest via rnai. our approach systematically extends and generalizes previous work by markowetz et al: an inference scheme was developed, which can deal with p-values for differential gene expression and does not rely on discretized data only. regularization was employed to incorporate prior assumptions on the network architecture into our framework. finally, new algorithms for large scale inference of signaling pathways were developed and evaluated in a systematic fashion on artificially created data. thereby, our module network, which recursively build up the complete topology from smaller pieces, revealed the best performance in terms of sensitivity and specificity. we used the module network to infer the signaling pathway for genes in the er-α pathway in human mcf- breast cancer cells and used a bootstrapping approach to ensure the statistical stability of the result. the induced edges in our inferred network were found with good consistency and could in many cases be also confirmed by the literature. future biological experiments are planned to validate our reconstructed network in a systematic way. in conclusion of our results we think that our approach offers a scale-able, reliable and fairly general way for large scale inference of signaling pathways from secondary effects and therefore provides researchers with a valuable tool to gain insight into complex cellular processes. the code for the module network inference method is available in the latest version of the -r-package nem, which can be obtained from the bioconductor homepage. 
in this paper we address this challenging problem by extending previous work by markowetz et al., who proposed a statistical framework to score networks hypotheses in a bayesian manner. using a bootstrapping approach this reconstruction can be found with good statistical stability. our proposed module network is employed to infer the signaling network between genes in the er-α pathway in human mcf- breast cancer cells. the code for the module network inference method is available in the latest version of the r-package nem, which can be obtained from the bioconductor homepage. third and most important, we propose methods to scale up the original approach, which is limited to around genes, to large scale networks. these secondary effects can be used to computationally reverse engineer features of the upstream signaling pathway. second, we show how prior assumptions on the network structure can be incorporated into the scoring scheme using regularization techniques. in combination with dna microarray technology this enables researchers to gain insights into signaling pathways by observing downstream effects of individual knock-downs on gene expression. our extensions go in three directions: first, we introduce a way to omit the data discretization step needed in the original framework via a calculation based on p-values instead.
 the advent of rna interference techniques enables the selective silencing of biologically interesting genes in an efficient way. in combination with dna microarray technology this enables researchers to gain insights into signaling pathways by observing downstream effects of individual knock-downs on gene expression. these secondary effects can be used to computationally reverse engineer features of the upstream signaling pathway.  in this paper we address this challenging problem by extending previous work by markowetz et al., who proposed a statistical framework to score networks hypotheses in a bayesian manner. our extensions go in three directions: first, we introduce a way to omit the data discretization step needed in the original framework via a calculation based on p-values instead. second, we show how prior assumptions on the network structure can be incorporated into the scoring scheme using regularization techniques. third and most important, we propose methods to scale up the original approach, which is limited to around genes, to large scale networks.  comparisons of these methods on artificial data are conducted. our proposed module network is employed to infer the signaling network between genes in the er-α pathway in human mcf- breast cancer cells. using a bootstrapping approach this reconstruction can be found with good statistical stability. the code for the module network inference method is available in the latest version of the r-package nem, which can be obtained from the bioconductor homepage. 