Guillaume Blin
email: gblin@univ-mlv.fr

Eric Blais
email: eblais@mcb.mcgill.ca

Danny Hermelin

Pierre Guillon
email: pguillon@univ-mlv.fr

Mathieu Blanchette
email: blanchem@mcb.mcgill.ca

Nadia El-Mabrouk

Gene Maps Linearization using Genomic Rearrangement Distances

A preliminary step to most comparative genomics studies is the annotation of chromosomes as ordered sequences of genes. Different genetic mapping techniques often give rise to different maps with unequal gene content and sets of unordered neighboring genes. Only partial orders can thus be obtained from combining such maps. However, once a total order O is known for a given genome, it can be used as a reference to order genes of a closely related species characterized by a partial order P . Our goal is to find a linearization of P that is as close as possible to O, in term of a given genomic distance. We first prove NP-completeness complexity results considering the breakpoint and the common interval distances. We then focus on the breakpoint distance and give a dynamic programming algorithm whose running time is exponential for general partial orders, but polynomial when the partial order is derived from a bounded number of genetic maps. A time-efficient greedy heuristic is then given for the general case

Introduction

Despite the increase in the number of sequencing projects, the choice of candidates for complete genome sequencing is usually limited to a few model organisms and species with major economical impact. For example, the rice genome is the only crop genome that has been completely sequenced. Other grasses of major agricultural importance such as maize and wheat are unlikely to be sequenced in the short term, due to their large size and highly repetitive composition. In this case, all we have are partial maps produced by recombination analysis, physical imaging and other mapping techniques that are inevitably missing some genes (or other markers) and fail to resolve the ordering of some sets of neighboring genes. Only partial orders can thus be obtained from combining such maps. The question is then to find an appropriate order for the unresolved sets of genes. This is important not only for genome annotation, but also for the study of evolutionary relationships between species.

Once total orders have been identified, the classical genome rearrangement approaches can be used to infer divergence histories in terms of global mutations such as inversions, transpositions and translocations [START_REF] Bergeron | Reversal distance without hurdles and fortresses[END_REF]Bourque et al., 2005b;[START_REF] El-Mabrouk | Sorting signed permutations by reversals and insertions/deletions of contiguous segments[END_REF]Hannenhalli and Pevzner, 1999;[START_REF] Pevzner | Human and mouse genomic sequences reveal extensive breakpoint reuse in mammalian evolution[END_REF][START_REF] Tang | Phylogenetic reconstruction from gene rearrangement data with unequal gene contents[END_REF].

In a recent study, Sankoff et. al. generalized the rearrangement by reversal problem to handle two partial orders [START_REF] Sankoff | Reversals of fortune[END_REF][START_REF] Zheng | Reversal distance for partially ordered genomes[END_REF]. The idea was to find two total orders with the minimal reversal distance. The problem has been conjectured NP-hard, and a branch-and-bound algorithm has been developed for this purpose. The difficulty of this problem is partly due to the fact that both compared genomes have partially resolved gene orders. However, once a total order is known for a given genome (for example a completely sequenced genome), it can be used as a reference to order markers of closely related species.

Given a reference genome characterized by a total order O and a related genome characterized by a partial order P , the goal is to find a permutation consistent with P minimizing a given genomic distance with respect to O. We consider the breakpoint distance, which is the simplest measure of gene order conservation usually used as a first attempt to solve any genome rearrangement problem. Moreover, half the breakpoint distance gives a lower bound for the inversion distance. We also consider a more general measure of synteny, the common interval distance, that has been widely studied in the last two years [START_REF] Blin | Conserved interval distance computation between nontrivial genomes[END_REF][START_REF] Figeac | Sorting by reversals with common intervals[END_REF][START_REF] Bérard | Conservation of combinatorial structures in evolution scenarios[END_REF]Bourque et al., 2005a;[START_REF] Blin | Inferring positional homologs with common intervals of sequences[END_REF].

After introducing the basic concepts in Section 2, we prove NP-complete results in Section 3. Section 4 then presents two exact dynamic programming algorithms for the breakpoint distance. The first algorithm applies to arbitrary partial orders. Its running time is exponential in the number of genes, but when the partial order is the intersection of a bounded number of genetic maps of bounded width, it runs in polynomial time. The second is a linear time algorithm that applies only to partial orders derived from a single genetic map. Section 5 then presents a fast and accurate heuristic for the general case. We finally report results on simulated data, and applications to grass genetic maps in Section 6.

A graph representation of gene maps

Hereafter, we refer to elementary units of a genetic map as genes, although they could in reality be any kind of markers. Moreover, as the transcriptional orientation of genes is usually missing from genetic maps, we consider unsigned genes.

Formally, a genetic map is represented as an ordered sequence of gene subsets or blocks B 1 , B 2 , . . . , B q , where for each 1 ≤ i ≤ q, genes belonging to block B i are incomparable among themselves, but precede those in blocks B i+1 , . . . , B q and succeed those in blocks B 1 , . . . , B i-1 . For example, in Figure 1.a, {4, 5} is a block, indicating that the order of genes 4 and 5 is left undecided. Maps M 1 , . . . , M m obtained from various protocols can be combined to form a more complex partial order P on the union set of the genes of all maps as follows: a gene a precedes a gene b in P if there exists a map M i where a precedes b.

However, combining maps can be a problem in itself, due to possible inconsistencies which would create precedence cycles (e.g. a precedes b in M 1 but b precedes a in M 2 ), or the presence of multiple loci (markers that are assigned to different positions in the same map).

These issues have been considered in previous studies [START_REF] Yap | A graph-theoretic approach to comparing and integrating genetic, physical and sequence-based maps[END_REF][START_REF] Zheng | Reversal distance for partially ordered genomes[END_REF][START_REF] Sankoff | Reversals of fortune[END_REF], and a software is available for combining genetic maps [START_REF] Jackson | Consensus genetic maps: a graph theory approach[END_REF]. In this paper, we assume that the partial order P is already known.

A partial order P is represented as a DAG (directed acyclic graph) (V P , E P ), where V P is the set of vertices (genes) and E P is the edge set representing the available order information (Figure 1 using data from "IBM2 neighbors 2004" [START_REF] Polacco | IBM neighbors: a consensus GeneticMap[END_REF], "Paterson 2003" [START_REF] Bowers | A high-density genetic recombination map of sequencetagged sites for Sorghum, as a framework for comparative structural and evolutionary genomics of tropical grains and grasses[END_REF] and "Klein 2004" [START_REF] Menz | A High-Density Genetic Map of Sorghum Bicolor (L.) Moench Based on 2926 Aflp, Rflp and Ssr Markers[END_REF]). E P is a minimum set of edges, in the sense that an edge can not be deduced by transitivity from others.

Preliminary definitions

Let P be a partial order represented by a DAG (V P , E P ). A vertex a is P -adjacent to a vertex b (denoted by a < P b) if there is an edge from a to b, and a precedes b (denoted

Problems

Formally we define the two following problems:

Minimum-Breakpoint Linearization (MBL) problem

Given: A partial order P and a permutation O on the set of genes {1, 2, . . . , n}, Find: A linearization of P into a permutation O so that bkpts(O, O ) is minimized.

Maximum-Common Interval Linearization (MCIL) problem

Given: A partial order P and a permutation O on the set of genes {1, 2, . . . , n}, Find: A linearization of P into a permutation O so that ICommon(O, O ) is maximized.

One may note that minimizing the number of breakpoints is equivalent to maximizing the number of O-adjacencies. Moreover, as an adjacency is a common interval of size 2, the MCIL problem is a generalization of the MBL problem. W.l.o.g, we assume from now on that O is the identity permutation (1, 2, . . . , n).

Remark 1 ( Signed genes and unequal gene content) All hardness results and algorithmic solutions developed in this paper hold for signed genes as well, using the classical

definition of breakpoints: a breakpoint of O with respect to O is a pair (a, b) such that a < O b, but neither a ≮ O b nor -b ≮ O -a.
As for common intervals, the same definition holds in the case of signed permutations. Moreover, from a theoretical point of view, all algorithmic solutions developed are also applicable to a partial order and a permutation with unequal gene content. However, as the goal is to find an appropriate order of a genetic map G using a reference genome H, only common genes of G and H are of interest.

In this section, we prove that the decision version of both the MBL and the MCIL problems is NP-complete. The properties of partial orders linearization that are presented in this section and that are required to prove the hardness result are also useful in the understanding of the two exact dynamic programming algorithms presented in the next section.

We propose a reduction from the NP-complete problem Maximum Independent Set [START_REF] Garey | Transforming cabbage into turnip (polynomial algorithm for sorting signed permutations by reversals)[END_REF]: given a graph G = (V, E) and an integer k, can one find an independent set of vertices of G -i.e. a set V ⊆ V such that no two vertices of V are connected by an edge in E -of cardinality greater than or equal to k ?

We initially note that the MBL and the MCIL problems are in NP since given a permutation O and a linearization O of P , one can compute the number of breakpoints in linear time and the number of common intervals in quadratic time.

For convenience, we define a reduction from a slightly different set of instances for the Maximum Independent Set problem: connected graphs. This can be done w.l.o.g. since the problem is still NP-complete in that case. Let G = (V, E) be a connected graph of n vertices. We define the permutation O and the partial order P as follows.

The permutation O is defined as a string

O = δ α 1 β 1 γ 1 α 2 β 2 γ 2 . . . α n β n γ n χ ,
and the partial order P as a DAG P = (V P , E P ) with

V P = {δ, α 1 , α 2 , . . . , α n , β 1 , β 2 , . . ., β n , γ 1 , γ 2 , . . . , γ n , χ, } and E P = {(γ 1 , δ), (δ, χ)} ∪ {(γ i+1 , γ i )|1 ≤ i < n} ∪ {(χ, α i ), (χ, β i )| 1 ≤ i ≤ n} ∪ {(β i , α j ), (β j , α i )|∀(v i , v j ) ∈ E} ∪ {(α i , ), (β i , )|1 ≤ i ≤ n}.
In the following, we will refer to any such construction as a PO-construction. An illustration of a PO-construction of a graph G with 6 vertices is illustrated in Figure 2.

First, let us present some interesting properties of any instance (O, P ) obtained by a PO-construction (Lemmas 1 and 2) relative to the breakpoints distance. Then, we will present other interesting properties (Lemmas 3 and 4) relative to common intervals. Finally, we will use those properties to prove that both the Minimum-Breakpoint Linearization and the Maximum-Common Interval Linearization problems are NP-complete. Lemma 1 Let G = (V, E) be a graph and P = (V P , E P ) be a partial order obtained from G by a PO-construction. There exists no linearization O of P where both α i < O β i and

α j < O β j , for any α i , α j , β i , β j ∈ V P such that (v i , v j ) ∈ E.
Proof. By contradiction, let us assume that there exists such a linearization O where α i < O β i and α j < O β j . Since (v i , v j ) ∈ E, we have (β i , α j ) ∈ E P and (β j , α i ) ∈ E P . Therefore, in any linearization of P -and consequently Oβ i O α j and β j O α i -which leads, by transitivity, to

β i O α i ; a contradiction. 2 tion O of P , bkpts(O, O ) = (3n + 2) -k where k is the number of pairs (α i , β i ) such that α i < O β i . Proof. By construction, in O, (i) δ < O α 1 < O β 1 < O γ 1 , (ii) ∀1 < i ≤ n, γ i-1 < O α i < O β i < O γ i and (iii) γ n < O χ < O . In any linearization O of P , (i) γ 1 < O δ, (ii) ∀1 < i ≤ n, γ i < O γ i-1 , (iii) δ < O χ and (iv) either α j < O or β j < O for a given 1 ≤ j ≤ n.
Therefore, in any linearization O of P , the only adjacencies that can be preserved are the ones of the form α i < O β i for some 1 ≤ i ≤ n. Let k be the number of pairs (α i , β i ) such that

α i < O β i . If k = 0 then no adjacencies at all are preserved, therefore bkpts(O, O ) = (3n+2). Consequently, if k > 0 then bkpts(O, O ) = (3n + 2) -k. 2
Lemmas 1 and 2 will be used afterwords in the proof of the NP-completeness of the MBL-problem. We now state some properties relative to the types and maximum numbers of common intervals between the permutation and any linearization of the partial order obtained by a PO-construction.

Lemma 3 Let G = (V, E) be a graph of n vertices, O and P = (V P , E P ) be respectively a permutation and a partial order obtained from G by a PO-construction. Given any linearization O of P , any non-trivial common interval between O and O is one of the following:

{[δ, χ], [δ, ], [α i , β i ]}, where 1 ≤ i ≤ n .
Proof. Let us first characterize any common interval of size greater or equal to 3. By construction, considering O as the reference, any interval of that size is of one of the following forms: (i) {α j , β j , γ j }, (ii)

{β i , γ i , α i+1 }, (iii) {γ i , α i+1 , β i+1 }, (iv) {γ n , χ, }, (v) {δ, α 1 , β 1 } for 1 ≤ j ≤ n and 1 ≤ i < n.
Since, by construction, in any linearization

O of P γ i O δ < O χ O {α i , β i } for 1 ≤ i ≤ n,
any common interval containing an α and a γ has to contain also χ and δ.

Therefore, there is no common interval of type (i), (ii) or (iii); the smallest common interval including simultaneously an α and a γ is indeed [δ, χ].

By construction, in any linearization O of P , γ n O δ < O χ. Thus, any common interval of type (iv) has to contain also δ. Therefore, there is no common interval of type (iv); the smallest common interval including simultaneously γ n , and

χ is indeed [δ, ].
Finally, in any linearization

O of P , δ < O χ O {α i , β i } for any 1 ≤ i ≤ n. Therefore,
any common interval of type (v) has to contain also χ. Therefore, there is no common interval of type (v); the smallest common interval including simultaneously δ and an α or a

β is indeed [δ, χ].
We just proved that any common interval of size greater or equal than 3 is either [δ, χ] or [δ, ]. Let us now characterize the common intervals of size 2. By construction any of these intervals is of one of the following forms: (i) {δ, α 1 }, (ii) {α i , β i }, (iii) {β i , γ i }, (iv) {γ j , α j+1 }, (v) {γ n , χ}, (vi) {χ, } for 1 ≤ i ≤ n and 1 ≤ j < n.

By construction, in any linearization

O of P , γ i O δ < O χ O {α i , β i } for any 1 ≤ i ≤ n.
Therefore, any common interval of type (i), (iii), (iv) or (v) has to contain also χ and δ. Therefore, there is no common interval of those types. Finally, in any linearization

O of P , χ O {α i , β i } O . Whereas in O, β n < O γ n < O χ < O .
Therefore, any common interval of type (vi) has to contain also δ and γ n .

We have thus proved that given any linearization O of P , any common interval of size at least two between O and O is one of the following:

{[δ, χ], [δ, ], [α i , β i ]} where 1 ≤ i ≤ n.
Lemma 4 Let G = (V, E) be a graph of n vertices, O and P = (V P , E P ) be respectively a permutation and a partial order obtained from G by a PO-construction. Given any lineariza-

tion O of P , ICommon(O, O ) = k + 3n + 5
where k is the number of pairs

(α i , β i ) such that α i < O β i or β i < O α i . Proof. Let (v i , v j ) ∈ E.
By construction, in P we have β j < P α i and

β i < P α j . Let O be a linearization of P s.t. α i < O β i . In O , [α i , β i ] is thus a common interval between O and O .
We will prove that [α j , β j ] cannot be a common interval between O and O .

Indeed, there are three cases: α j < O β j , β j < O α j or none of those two. By Lemma 1, the first case is not possible. Therefore, consider first that β j < O α j . Then, in O we have

β j O α i < O β i O α j . Consequently, [α j , β j ] cannot be a common interval. Similarly, consider that in O , β i < O α i . In O , [α i , β i ] is a common interval between O and O . If α j < O β j , in O we have β i O α j < O β j O α i ; [α i , β i ]
is not a common interval anymore, a contradiction. If β j < O α j then by construction we have β j O β i < O α i O α j . In both cases, [α j , β j ] cannot be a common interval. Therefore, by Lemma 3, in any linearization of P , there are k + 2 + 3(n + 1) common intervals where k is the number of pairs (α i , β i ) such

that α i < O β i or β i < O α i .
2

We now turn to the proof of the following theorems.

Theorem 1 A connected graph G = (V, E) admits an independent set of vertices V ⊆ V of cardinality greater than or equal to k if and only if there exists a linearization O of P such that bkpts(O, O ) ≤ (3n + 2)k, where O and P result from a PO-construction of G.

Proof. (⇒) Let V ⊆ V such that |V | ≥ k and V is an independent set. Let O be a linearization of P defined by O = P 1 δ χ P 2 P 3 P 4 where:

• P 1 is the linearization of the subset of vertices V 1 = {γ i |1 ≤ i ≤ n} such that ∀1 < i ≤ n, γ i < P 1 γ i-1 ; • P 2 is any linearization of the subset of vertices V 2 = {β i |v i ∈ V -V }; • P 3 is any linearization of the subset of vertices V 3 = {α i , β i |v i ∈ V } such that ∀v i ∈ V , α i < P 3 β i ; • P 4 is any linearization of the subset of vertices V 4 = {α i |v i ∈ V -V }.
For example, given the instances illustrated in Figure 2 and an independent set V = {1, 3, 4}, O is defined by:

O = γ 6 γ 5 γ 4 γ 3 γ 2 γ 1 δ χ β 2 β 5 β 6 α 1 β 1 α 3 β 3 α 4 β 4 α 2 α 5 α 6 .
By (⇐) Suppose we have a linearization O of P such that bkpts(O, O ) ≤ (3n + 2)k. Let V ⊆ V be the set of vertices such that:

∀(α i , β i ) such that α i < O β i , add v i to V
By Lemma 1, we can affirm that V is an independent set. Let us verify that |V | ≥ k. Let V ⊆ V be the set of vertices such that:

By

∀(α i , β i ) such that α i < O β i or β i < O α i , add v i to V
By Lemma 3, we can affirm that V is an independent set. Let us verify that |V | ≥ k.

By Lemma 4, ICommon(O, O ) = k + 5 + 3n where k is the number of pairs (α i , β i ) such that α i < O β i or β i < O α i . Therefore, we obtain |V | = k 2
4 Exact dynamic programming algorithms

Hereafter, we describe two exact dynamic programming algorithms for solving the MBL problem. The first algorithm works on an arbitrary partial order P , but has a running time that can be exponential in |V P |. However, we show that the algorithm's running time is polynomial in the more realistic case where P is built from a bounded set of genetic maps of bounded width. The second algorithm applies to the case where P is built from a single genetic map, and runs in linear time.

We begin with some preliminary definitions. Let A be a subset of vertices of V P . A is a border of P iff any pair of vertices of A are incomparable, and a maximal border iff any other vertex of V P is comparable to at least one vertex of A. We also define, for any subset

B ⊆ V P , f ront(B) = {x ∈ B :
x has no successor in B}. Finally, for any subset A ⊆ V P , we denote pred(A) = A ∪ {x ∈ V P : ∃y ∈ A s.t. x P y}.

A dynamic algorithm for arbitrary partial orders

Let A be a border. We denote by X A,i the maximum number of adjacencies that can be obtained from a linearization of pred(A) that is consistent with the partial order P , and that ends with vertex i (i.e. i is the rightmost vertex in the total ordering of pred(A)). It is easy to see that the number of adjacencies in the global optimal solution is max i∈F X F,i adjacencies, where F = f ront(V P ). The following theorem provides a recursive formula for the computation of X A,i .

Theorem 3 For any border A and any vertex i ∈ A,

X A,i = max j∈A X A ,j +        1 if |j -i| = 1 0 otherwise where A = f ront(pred(A) \ {i})) = (A \ {i}) ∪ {k | (k, i) ∈ E P and k ∈ pred(A \ {i})}
A recursive algorithm follows from the previous theorem. The recursion begins with A = F = f ront(V P ) and stops as soon as A is the empty set.

Computing each entry of the dynamic programming table only requires operations which can be done in linear time. If the partial order P admits b(P ) possible borders, the running time is O(b(P ) • |V P | 2 ). In the general case, the number of borders of P can be as much as 2 |V P | , if P consists of a single block of incomparable vertices. However, we are more interested in the case where P is obtained by combining a small number m of genetic maps,

where each map contains a maximum of q blocks and the size of each block is at most some Furthermore, as much as possible, we want to join consecutive elements in neighboring blocks as well. The set L i is used to keep track of which elements can be put last in the ordering of B i and therefore possibly be matched with an element in the block B i+1 . If the elements of B i are stored in an ordered list, then the recursive definition of Theorem 4 can be implemented in a recursive algorithm for which each iteration requires O(|B i | + |L i-1 |) time to run, for a total time complexity of O(n) in the case of a genetic map of n genes.

An efficient heuristic

Since our exact dynamic programming for the general problem has a worst-case running time that is exponential in the number of genes, a faster heuristic is required to solve large problem instances. In this section, a greedy heuristic is developed for general partial orders obtained from the intersection of an arbitrary number of maps. It aims to find a maximum number of O-adjacencies consistent with a partial order P . At each step, the partial order is updated by incorporating adjacencies of the longest O-adjacency path that can be part of a linearization of P . The algorithm does not necessarily end up with a total order. Rather, it stops as soon as no more adjacencies can be found. All linearizations of the obtained partial order are then equivalent in the sense that they all give rise to the same number of adjacencies.

A direct (resp. indirect) adjacency path of P is a sequence of vertices of form (i,

i + 1, i + 2, • • • i + k) (resp. (i + k, • • • i + 2, i + 1, i))
such that for any 0 ≤ j < k, either i + j < P i + j + 1 (resp. i + j + 1 < P i + j), or i + j and i + j + 1 are incomparable. For example, in Figure 3,

(1, 2, 3, 4) (resp. (11, 10, 9, 8, 7)) is a direct (resp. indirect) adjacency path. Notice that adjacencies of this indirect path can not belong to any linearization of P , as gene 5 should be located after 11 but before 7.

We say that an adjacency path p of P is valid iff there is a linearization O of P such that p is a subsequence of O . For example, (1, 2, 3, 4) is a valid path of the partial order in Figure 3. Lemma 5 gives the conditions for an adjacency path to be valid. We need a preliminary definition.

Definition 1 Given two vertices i and j, we say that i is compatible with j iff the two following conditions hold:

1. i and j are either incomparable or i P j;

2. Any vertex v verifying i P v P j belongs to the interval [i, j] (or [j, i] if j < i).

Lemma 5 A direct (resp. indirect) adjacency path of P from i to i + k (resp. from i + k to i) is valid if and only if, for any j 1 , j 2 such that 0 ≤ j 1 < j 2 ≤ k, i + j 1 is compatible with i + j 2 . (resp. i + j 2 is compatible with i + j 1 ).

A preliminary preprocessing of P = (V P , E P ) is required to efficiently compute successive adjacency paths.

1. Create the matrix M of size |V P | × |V P | verifying, for any i, j ∈ V P , M (i, j) = 1 iff i < P j and M (i, j) = 0 otherwise. 2. Compute the transitive closure of M , that is the matrix M T of size |V P |×|V P | verifying, for any i, j ∈ V P ,

M T (i, j) =    1 iff i < P j 2 iff i P j but i ≮ P j 0 otherwise
M T is computed from M using the Floyd-Warshall algorithm [START_REF] Floyd | Algorithm 97: Shortest path[END_REF].

After the preprocessing step, the following Steps 1 and 2 are iterated as long as P contains an adjacency path.

•

Step 1: Find a longest valid direct or indirect adjacency path (see details below).

• The algorithm for valid indirect paths is obtained by replacing the three first lines of Algorithm 4 by the following:

For i = |V | to 1 do If (i < P i -1)
or (i and i -1 are incomparable) then Add (i, i -1) to L;

Step 1 consists in running successively both algorithms for direct and indirect paths, and taking the longest resulting path. probability of possible O-adjacencies. We simulated twenty different instances for each triplet of parameters (n, p, q) with k ∈ {30, 50, 80, 100}, p ∈ {0.7, 0.9} and q ∈ {0.4, 0.6, 0.8}. We did not consider p values lower than 0.7, as the dynamic programming algorithm exponentialtime prevented us from testing such instances.

Both the heuristic and the dynamic programming algorithm were run on the dataset described above. We evaluated two criteria: the CPU time and the number of breakpoints (or similarly adjacencies) induced by the returned linearization.

Figure 5 shows that the running time of the dynamic programming algorithm grows exponentially with the width of the DAG (defined as the size of the DAG's largest border), while the heuristic is not affected by it (this was expected, as the time-complexity only depends on the DAG's size). Moreover, the greedy heuristic can easily handle partial orders consisting of thousands of genes, as illustrated in Figure 6.

We now evaluate the heuristic's optimality, namely the number of O-adjacencies resulting from the obtained linearization compared to the optimal solution (obtained with the dynamic programming algorithm). As illustrated by Table 1, the greedy algorithm always returns a linearization containing at least 90% of the adjacencies of the optimal solution, and usually close to 100%.

Illustration on grass genomes: The Gramene database contains a large variety of maps of different grass genomes such as rice, maize and oats, and provides tools for comparing individual maps. A visualization tool allows to identify regions of 'homeology' between species, that is a linear series of markers in one genome that maps to a similar series of loci on another genome. Integrating marker orders between different studies remains a challenge to geneticists. However, as total orders are already obtained for widely studied species such as rice, which has been completely sequenced, one can use this information to order markers on another species by using the adjacency maximization criterion. Extracting the linear orders of markers using the Gramene visualization tool remains unpractical for hundreds of markers, as no automatic tool is provided for this purpose. We therefore illustrate the method on maps that are small enough to be extracted manually.

Maize has been chosen instead of rice as it has shorter maps, though non-trivial, that can be represented graphically.

We used the "IBM2 Neighbors 2004" [START_REF] Polacco | IBM neighbors: a consensus GeneticMap[END_REF] map for chromosomes 5

(Figure 1) and 1 (Figure 7) of maize as a reference, and compared it with the "Paterson 2003" [START_REF] Bowers | A high-density genetic recombination map of sequencetagged sites for Sorghum, as a framework for comparative structural and evolutionary genomics of tropical grains and grasses[END_REF] and "Klein 2004" [START_REF] Menz | A High-Density Genetic Map of Sorghum Bicolor (L.) Moench Based on 2926 Aflp, Rflp and Ssr Markers[END_REF] maps of the chromosomes labeled C and LG-01, respectively, of sorghum. We extracted all markers of maize indicated as having a homolog in one of the databases of sorghum. All are found completely ordered in maize. This linear order is considered as the identity permutation. For markers of sorghum that are located on maize chromosome 5 (resp. 1), a total order maximizing the adjacency criterion is indicated in Figure 1.d (resp.

Conclusion

We have presented a detailed complexity result and algorithmic study for the problem of linearizing a partial order that is as close as possible to a given total order, in term of the breakpoint and common intervals distances. Applications on the grass genomes show that this may be helpful to order unresolved sets of markers of some species using the totally ordered maps of well studied species such as rice. However, preliminary to the application of our algorithms is generating the appropriate partial orders. For this purpose, an automated preprocessing of the Gramene comparative database would be required to output the considered genetic maps, and then combine them on a single partial order. The absence of such tools prevented us from presenting more suitable applications.

The next step of this work will be to generalize our approach to two (or more) partial orders, as previously considered in [START_REF] Sankoff | Reversals of fortune[END_REF][START_REF] Zheng | Reversal distance for partially ordered genomes[END_REF] for the reversal distance. As conjectured by Sankoff, an NP-complete result for this problem should be proved. A dynamic programming approach may also be developed for this case.

An adjacency of two genes being just a common interval of size 2, a simple extension of the greedy heuristic would be to order genes that remain unordered after maximizing adjacencies, by using the constraint of maximizing intervals of size 3, 4 and so on. An

Figure 1 :

 1 Figure 1: Data extracted from the comparison of maize and sorghum in the Gramene

Figure 2 :

 2 Figure 2: Example of a PO-construction. The graph (a) is a connected graph of 6 vertices.

  Lemma 2, we can affirm that bkpts(O, O ) = (3n + 2) -|V |. Since, by hypothesis, |V | ≥ k, we obtain bkpts(O, O ) ≤ (3n + 2)k.

  Lemma 2, bkpts(O, O ) = (3n + 2)k where k is the number of pairs (α i , β i ) such that α i < O β i . Therefore, we obtain |V | = k. 2 Theorem 2 A connected graph G = (V, E) admits an independent set of vertices V ⊆ V of cardinality greater than or equal to k if and only if there exists a linearization O of P such that ICommon(O, O ) ≥ k + 3n + 5, where O and P result from a PO-construction of G. Proof. (⇒) The proof is almost the same as for Theorem 1. Let V ⊆ V such that |V | ≥ k and V is an independent set. Let O be a linearization of P defined as in the proof of Theorem 1 (i.e. O = P 1 δ χ P 2 P 3 P 4 ). By Lemma 4, we can affirm that ICommon(O, O ) = |V | + 5 + 3n. Since, by hypothesis, |V | ≥ k, we obtain ICommon(O, O ) ≥ k + 5 + 3n. (⇐) Suppose we have a linearization O of P such that ICommon(O, O ) ≥ k + 5 + 3n.

Figure 3 :

 3 Figure 3: Dotted edges are all O-adjacencies that can, individually, be part of a linearization of P . A bi-directional edge represents the superposition of two edges, one in each direction. An adjacency path of P is a directed sequence of consecutive dotted edges.

AlgorithmFigure 4 :

 4 Figure4: Finding a longest valid adjacency path of P . L is the list of adjacency paths of size k, L j denotes the j th path of L, and L j [i] the i th vertex of L j .

  Step 2: Incorporate the new adjacencies in M T , and compute the transitive closure of M T . Algorithm 4 describes the search of the longest valid direct path. Valid direct paths are computed beginning with paths of size 2. For a fixed k, any path p = (i, i+1, • • • i+k) of size k is obtained from a concatenation of two valid consecutive paths p 1= (i, i + 1, • • • i + k -1) and p 2 = (i + 1, i + 2, • • • i + k) of size k -1.As p 1 and p 2 are valid paths, the path p is valid iff i is compatible with i + k.

Figure 5 :

 5 Figure 5: CPU time expended by (a) the dynamic programming algorithm and (b) the

Figure 6 :

 6 Figure 6: CPU time expended by the heuristic, for DAGs of 1000 vertices with different

FigureFigure 7 :

 7 Figure 7: (a) The partial order of markers in sorghum that are located and totally ordered

Table 1 :

 1 Percentage of O-adjacencies resulting from the heuristic's linearization compared to the optimal solution (obtained with the dynamic programming algorithm). Results are obtained by running the heuristic and dynamic programming algorithm on 10 different simulated DAGs for a given size and width.

		3	4	5	6	7	8	9	10	11
		30 100 100 98,15 97,41 94,33 96,18 93,18 95,54 100	100
	Genome size	50 100 98,04 96,43 97,62 95,25 98,61 100 86,96 95,26 94,94 80 100 98,21 97,90 87,54 96,79 93,89 100 95,83 98,33 100
		100 100 98,81 95,65 96,83 89,70 93,95 90,38 95,30 94,63 94,95

small number k. In this case, there are at most q m maximal borders in P . Furthermore, two elements that are in the same border cannot be in different blocks on a genetic map, so each maximal border is of size at most km, which allows 2 km possible subsets. Therefore, the total number of borders of P is bounded above by b(P ) ∈ O(q m • 2 km ). Since, in practice, only two or three different genetic maps are combined to form a partial order, the dynamic algorithm yields a practical and exact solution to the MBL problem.

A linear-time algorithm for single genetic map

When P is built from a single genetic map consisting of a list of blocks B 1 , B 2 , . . . , B q , a much faster linearization algorithm exists. Let X i be the maximum linearization score obtained in the partial subset

The maximum linearization score of P is thus equal to X q . Let L i represent the set of elements in B i that can be placed at the last position in a total ordering of B 1 ∪ • • • ∪ B i that achieves the score X i . Define the functions

the values X i and L i can be determined recursively as follows.

Theorem 4 Define X 0 = 0 and L 0 = {}. Then, for any

The intuition behind the recursive definition of X i is as follows: to get the maximum linearization score, we always want to join as many elements x, x + 1 within a same block.

Complexity: Computing the transitive closure of the adjacency matrix in the preprocessing phase, as well as in Step 2, is done using the Floyd-Warshall algorithm [START_REF] Floyd | Algorithm 97: Shortest path[END_REF] in

Experimental results

We first test the efficiency of the heuristic compared to the dynamic programming algorithm for general partial orders on simulated data, and then illustrate the method on grass maps obtained from Gramene (http://www.gramene.org/).

Simulated data:

We used simulated data to assess the performance of our greedy algorithm. We simulate DAGs of fixed size n that can be represented as a linear expression involving the operators '→' and ',' where P-adjacent genes are separated by a '→' and incomparable genes by a ','. Such a representation is similar to the one used in [START_REF] Lander | MAP-MAKER: an interactive computer package for constructing primary genetic linkage maps of experimental and natural populations[END_REF][START_REF] Yap | A graph-theoretic approach to comparing and integrating genetic, physical and sequence-based maps[END_REF]. For example, the DAG in Figure 1.c has the following string representation: , 15, 16, 17, 21} → {18, 19} → 20 DAGs are generated according to two parameters: the order rate p that determines the number of ',' in the expression, and the gene distribution rule q corresponding to the


Guillaume Blin
email: gblin@univ-mlv.fr

Florian Sikora
email: sikora@univ-mlv.fr

Stéphane Vialette
email: vialette@univ-mlv.fr

Querying Graphs in Protein-Protein Interactions Networks using Feedback Vertex Set

Keywords: Graph Query, Pattern-Matching, Dynamic Programming, Protein-Protein Interactions networks

 gave an algorithm (which is not implemented) to query graphs with a bounded treewidth in PPI networks (the treewidth of the query being involved in the time complexity). In this paper, we propose another algorithm for querying pattern in the shape of graphs, also based on dynamic programming and the color-coding technique. To transform graphs queries into trees without loss of informations, we use feedback vertex set coupled to a node duplication mecanism. Hence, our algorithm is FPT for querying graphs with a bounded size of their feedback vertex set. It gives an alternative to the treewidth parameter, which can be better or worst for a given query. We provide a python implementation which allows us to validate our implementation on real data. Especially, we retrieve some human queries in the shape of graphs into the fly PPI network.

I. INTRODUCTION

C ONTRARY to what was predicted years ago, the human genome project has highlighted that human complexity may not only rely on its genes (only 25 000 for human compared to the 30 000 and 45 000 for the mouse and the poplar respectively). This observation increased the interest in protein properties (e.g. their numbers, functions, complexity and interactions). Among other protein properties, the set of all their interactions for an organism, called Protein-Protein Interactions (PPI) networks, have recently attracted lot of interest. The number of reported interactions increases rapidly due to the use of various genome-scale screening techniques [START_REF] Gavin | Functional organization of the yeast proteome by systematic analysis of protein complexes[END_REF], [START_REF] Ho | Systematic identification of protein complexes in Saccharomyces cerevisae by mass spectrometry[END_REF], [START_REF] Uetz | A comprehensive analysis of protein-protein interactions in Saccharomyces cerevisae[END_REF]. Unfortunately, acquiring such valuable resources is prone to high noise rate [START_REF] Gavin | Functional organization of the yeast proteome by systematic analysis of protein complexes[END_REF], [START_REF] Reguly | Comprehensive curation and analysis of global interaction networks in saccharomyces cerevisiae[END_REF].

Comparative analysis of PPI tries to determine the extent to which protein networks are conserved among species. Indeed, it was observed that proteins functioning together in a pathway (i.e., a path in the interactions graph) or a structural complex (i.e., an assembling of strongly connected proteins) are likely to evolve in a correlated fashion and during evolution, all such functionally linked proteins tend to be either preserved or eliminated in a new species [START_REF] Pellegrini | Assigning protein functions by comparative genome analysis: protein phylogenetic profiles[END_REF].

In this article, we focus on the following related problem called GRAPH QUERY (formally defined later). Given a PPI network and a pattern with a graph topology, find a subnetwork of the PPI network that is as similar as possible to the pattern, in respect to the initial topology. Similarity is measured both in terms of sequence similarity and graph topology conservation.

Unfortunately, this problem is clearly equivalent to the NP-complete subgraph homeomorphism problem [START_REF] Garey | Computers and Intractability: a guide to the theory of NP-completeness[END_REF]. Recently, several techniques have been proposed to overcome the difficulty of this problem. By restricting the query to a path of length less than five, Kelley et al. [START_REF] Kelley | Conserved pathways within bacteria and yeast as revealed by global protein network alignment[END_REF] developped PathBlast, a software with a factorial time complexity which allows one consecutive mismatch. Later on, Shlomi et al. [START_REF] Shlomi | QPath: a method for querying pathways in a protein-protein interaction network[END_REF] proposed an alternative, called QPath, for querying paths in a PPI network which is based on the color-coding technique introduced by Alon, Yuster and Zwick [START_REF] Alon | Color coding[END_REF]. The use of this technique allows to define a fixedparameter tractable (FPT) algorithm parameterized by the size of the query. Recall that a parameterized problem is FPT if it can be determined in f (k)n O (1) time, where f is a function only depending on the parameter k, and n is the size of the input [START_REF] Downey | Parameterized Complexity[END_REF]. In addiction of being faster, QPath deals with longer paths (until size ten) and allows more flexibility by considering a bounded number of non-exact matches.

By restricting the query to a tree, Pinter et al. [START_REF] Pinter | Alignment of metabolic pathways[END_REF] proposed an algorithm that is restricted to forest PPI networks (i.e., collection of trees). Finally, Dost et al. [START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF] developed QNet, an algorithm to handle tree query in the general context of PPI networks. The authors also gave some theoretical results for querying graphs using the tree decomposition of the query.

Since QNet is the main reference in this field and is quite related to the work presented in this paper, let us present it briefly. QNet is an FPT algorithm for querying trees in a PPI network. The time complexity is

2 O(k) m ln( 1 ǫ ),
where k is the number of proteins in the query, m the number of edges of the PPI network and 1-ǫ the success probability (for any ǫ > 0). As QPath, QNet uses dynamic programming together with the color-coding technique. For querying graphs in a network, QNet uses, as a subroutine, an algorithm to query trees. To do so, they perform a tree decomposition (a formal definition of a tree decomposition can be found in [START_REF] Bodlaender | A tourist guide through treewidth[END_REF]). Roughly speaking, it is a transformation of a graph into a tree, a tree node (or a bag) can contain several graph nodes. There exists several algorithms to perform such a transformation. The treewidth of a graph is the minimum (among all decompositions) of the cardinality of the largest bag minus one. Computing the treewidth is, however, NP-Hard [START_REF] Arnborg | Complexity of finding embeddings in a k-tree[END_REF]. From this tree decomposition, the time

complexity of QNet is 2 O(k) n t+1 ln( 1 ǫ )) time,
where k is the size of the query, n is the size of the PPI network, t is the treewidth of the query, and 1ǫ is the success probability (for any ǫ > 0).

QNet is an algorithm for querying trees in a PPI network. A logical extension would be to query graphs. The authors of [START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF] provide a theoretical solution, without implementation and which depends on the treewidth of the query. We propose here an alternative solution, that uses the color-coding technique (Section II). We provide in Section III some experimental results. II. PADA1 AS AN ALTERNATIVE TO QNET In this section, we propose an alternative to QNet called PADA1 (Protein Alignment Dealing with grAphs). At a more general level, QNet and PADA1 use the very same approach: transform the query into a tree and find an occurrence of that tree in the PPI network by dynamic programming. However, whereas QNet uses tree decompositions, PADA1 combines feedback vertex sets together with nodes duplications (Algorithm GRAPH2TREE). Let note that independently, Cheng et al. use a similar technique to transform a graph into a tree in [START_REF] Cheng | Fast Alignments of Metabolic Networks[END_REF] in order to query a graph in a network. However, unlike our approach, in [START_REF] Cheng | Fast Alignments of Metabolic Networks[END_REF] the authors do not use nodes duplication and hence, it is not clear they can ensure that all the edges of the query are kept in the results (especially the ones linking the nodes belonging to the feedback vertex set). It is worth mentioning that, following the example of QPath and QNet, we will consider non-exact matches (i.e., allowing indels). Since we allow queries to be graphs, PADA1 is clearly an extension of QPath and a real alternative to QNet.

A. Transforming the query into a tree

Let us first present Algorithm GRAPH2TREE which transforms a graph G = (V, E) into a tree. Our transformation is lossless, and hence, one can reconstruct the graph starting from the tree. The main idea of Algorithm GRAPH2TREE is to transform the graph into a tree by iteratively finding a cycle C, duplicating a node of C, and finally breaking cycle C by deleting one of its edges (see Figure 1 for an illustration on how to break a cycle at vertex v 1 ). Central in our approach is thus the node duplication procedure (Algorithm DUPLICATE). For each u ∈ V , write d(u) for the set of all copies of vertex u including itself and N (u) for the set of all its neighbors. Let F denote the set of all nodes of G that have been duplicated at the end of Algorithm GRAPH2TREE, i.e., F = {v ∈ V : |d(v)| > 1}. The cardinality of F turns out to be an important parameter since, as we will prove soon, the overall time complexity of PADA1 mostly depends on |F | and not on the total number of duplications. Minimizing the cardinality of F is the well-known NP-complete FEEDBACK VERTEX SET problem [START_REF] Karp | Reducibility among combinatorial problems[END_REF]: Given a graph G, find a minimum cardinality subset of vertices with the property that removing of these vertices from the graph results in an acyclic graph.

Function GRAPH2TREE(G) 1 begin 2 for all u of V , d(u) ← {u}; 3 F V S ← FEEDBACKVERTEXSET(G);
We have implemented a "brute-force" algorithm for the FEEDBACK VERTEX SET problem. Once this set is computed, we duplicate a node of a cycle as long as there is a cycle in the graph. By definition of the feedback vertex set, there is at least one node in a cycle belonging to the feedback vertex set. The time complexity of the Algorithm GRAPH2TREE is dominated by the process of computing the feedback vertex set, which is in O(2 |V | × |E|) since there are 2 |V | potential subgraphs. Nevertheless this solution is practical since it is still running in seconds if |V | is smaller than twenty. Indeed, the overall complexity of PADA1 considerably limits the size of our graph query. However, one may also consider an efficient FPT algorithm such as the one of Guo et al. [START_REF] Guo | Compression-based fixed-parameter algorithms for feedback vertex set and edge bipartization[END_REF], using iterative compression, in addition of the quadratic kernalization of Thomassé [START_REF] Thomasse | A quadratic kernel for feedback vertex set[END_REF] in order to compute efficiently the feedback vertex set.

B. Tree matching

We now assume that the query has been transformed into a tree (with duplicated nodes) by Algorithm GRAPH2TREE, and hence we only consider tree queries from this point. We show that an occurrence of such a tree can be found in a PPI network by dynamic programming. Let us fix notations. PPI networks are represented by undirected edge weighted graphs G N = (V N , E N , w) ; each node of V N represents a protein and each weighted edge (v i , v j ) ∈ E N represents an interaction between two proteins. A query is given by a tree T Q = (V Q , E Q ) (output of Algorithm GRAPH2TREE on the graph query). The set V Q represents proteins while E Q represents interactions between these proteins. As in QNet, we do not give weight on the query. Indeed, in PPI networks, the weight represents probabilities of the interactions. There is no clue to use those probabilities in the query. In the following, we will consider that T Q is an ordered tree where q 1 , . . . , q nq denote the ordered n q children of any node q. As we will show afterwards, the result does not depend on this ordering.

Function FEEDBACKVERTEXSET(G = (V, E)) 1 begin 2 for (i = 0 ; i < |V | ; i + +) do 3 foreach subgraph G ′ = (V ′ , E ′ ) of G such that |V ′ | = |V | -i do 4 if G ′ is acyclic then 5 return V \ V ′ ;
Function DUPLICATE(G = (V, E), v a , v b , d) 1 begin 2 Let i ← |d(v b )|; 3 Let v b,i be a new node; 4 V ← V ∪ {v b,i }; 5 d(v b ) ← d(v b ) ∪ {v b,i } ; 6 E ← E -{(v a , v b )}; 7 E ← E ∪ {(v a , v b,i )};
Let h(p 1 , p 2 ) be a function that returns a similarity score between two proteins p 1 and p 2 . The similarity considered here will be computed according to amino-acid sequences similarity (using BLASTp [START_REF] Altschul | Basic local alignment search tool[END_REF]). In the following, given two nodes v 1 and v 2 of V Q (or V N ), we write h(v 1 , v 2 ) for the similarity between the two proteins corresponding to v 1 and v 2 . A node v 1 is considered to be homologous to a node v 2 if the corresponding similarity score h(v 1 , v 2 ) is above a given threshold. Biologically, one can assume that two homologous proteins have probably similar functions. Clearly, for every node v of F , all nodes in d(v) are homologous with the same protein.

An alignment of the query T Q and G N is defined as:

(i) a subgraph G A = (V A , E A , w) ⊆ G N = (V N , E N , w), such that V A ⊆ V N and E A ⊆ E N , and (ii) a mapping σ : V Q → V A ∪ {del}.
More precisely, the function σ is defined such that if σ(q) = v then q and v are homologous.

For a given alignment of T Q and G N , a node q of V Q is said to be deleted if σ(q) = del and matched otherwise. Moreover, any node v a of V A such that σ -1 (v a ) is undefined is said to be inserted. Note that, similarly to QNet, only nodes of degree two can be deleted (we can only contract paths). For practical applications, the number of insertions (resp. deletions) is limited to be at most N ins (resp. N del ), each involving a penalty score δ i (resp. δ d ).

The GRAPH QUERY problem can be thus defined as follow: Given a query T Q with duplicated nodes, a PPI network G N , a similarity function h, penalty scores δ i and δ d for insertions and deletions, find an alignment (G A , σ) between T Q and G N of maximum score. The score of an alignment is defined as the sum of (i) similarity scores of aligned nodes (i.e.,

v∈V A σ -1 (v) defined h(v, σ -1 (v))
), (ii) the sum of weights of edges in E A (i.e., e∈E A w(e)), (iii) a penalty score δ d for each node deletion (i.e., q∈V Q σ(q)=del δ d ), and (iv) a penalty score δ i for each node insertion (i.e.,

v∈V A σ(v) -1 undefined δ i ).
The general problem is NP-complete. However, it is Fixed Parameter Tractable in case the query is a tree by a combination of the color-coding technique [START_REF] Alon | Color coding[END_REF] and dynamic programming. This randomized technique allows to find simple paths of length k in a network in O(2 k ) time (to be compared to the O(n k ) time brute-force algorithm), where n is the number of proteins in the network [START_REF] Scott | Efficient algorithms for detecting signaling pathways in protein interaction networks[END_REF]. In [START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF], the authors of QNet adapted this technique for their query algorithm. Since one is looking for an alignment, each node of the query has to be considered once (and only once) in an incremental build of the alignment by dynamic programming. Thus, one has to maintain a list of the nodes already considered in the query. Therefore, on the whole, one has to consider all O(n k ) potential alignments, with

n = |V N | and k = |V Q |.
Using color-coding, one may decrease this complexity to O(2 k ). First, nodes of the network are colored randomly using k colors, where k = |V Q |. Then, looking for a colorful alignment (i.e., an alignment that contains each color once) leads to a solution, which is not necessarily optimal. Therefore, one only needs to maintain a list of the colors already used in the alignment, storable in a table of size in O(2 k ). In order to get an optimal solution, this process is repeated. More precisely, according to QNet [START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF], since a colorful alignment happens with probability k! k k ≃ e -k , the coloration step has to be done ln( 1 ǫ )e k times to obtain an optimal alignment with high probability (1ǫ, for any ǫ).

The QNet dynamic programming algorithm can be summarized as follows. By an incremental construction, for each (q i , q j ) ∈ E Q when one is considering q i ∈ V Q aligned with a node v i ∈ V N , check whether the score of the alignment is improved through: (i) a match of q j and any v j of V N such that q j and v j are homologous and (v i , v j ) ∈ E N , (ii) an insertion of a node v j of V N in the alignment graph G A , and (iii) a deletion of q j . This is done for a given coloration of the network, and repeated for each coloration.

Hereafter, we define an algorithm, inspired from QNet, which consider a query tree T Q , a PPI network G N and seeks for an alignment (G A , σ). It is worth noticing that for a given coloration, our algorithm, as QNet, is exact. To deal with duplicated nodes (cf. GRAPH2TREE algorithm), we pre-compute all possible assignment, called A, of the duplicated nodes σ(q ′ ) = v ∀q ′ ∈ d(q). We then compute for each assignment A the score of an alignment with respect to A. We denote BESTCONSTRAINTALIGNMENT this step. The difficulty lies in the construction of the best alignment by dynamic programming, with respect to A.

V Q of T Q . More precisely, ∀q ∈ F, ∀v ∈ V N , define Function PADA1 (T Q ,G N , h, threshold) 1 begin 2 BestG A ← ∅; BestScore ← -∞;
As done in QNet, we use a set S C of k + N ins colors (as needed by color-coding) which will be used when a node is matched or inserted. Moreover, in order to deal with potential duplicated nodes in T Q , we have to use another multiset S of colors (i.e., the colors in this set can appear more than once), rather than a classical set as in QNet. Indeed, every node in d(q) such that q ∈ F , must use the same color.

As a preprocess to PADA1, G N can be pruned, as shown in [START_REF] Bruckner | Topology-free querying of protein interaction networks[END_REF]. Let u ∈ G N be a protein which is not homologous with any protein of the query and v ∈ G N be a protein which is homologous with a protein of the query. Then, u can be too far from any v in terms of shortest path length to be inserted in the solution in regards to the maximum number of insertions (i.e., N ins ). According to this remark, u is kept in G N only if there are two proteins v 1 and v 2 , both homologous with a protein of the query, such that dist(u, v 1 ) + dist(u, v 2 ) ≤ N ins + 1, where dist(u, v) is the length of the shortest path between u and v. Otherwise, u can never been in a solution, and hence can be safely deleted from G N .

Once G N has been pruned, PADA1 can be launched for each valid connected component of G N . A component is said to be valid if it contains at least k -N del proteins which are homologous with a protein of the query, where k is the size of the query. Otherwise, a solution can never be found in this component, and hence there is no need to consider it. As stated in [START_REF] Bruckner | Topology-free querying of protein interaction networks[END_REF], there is in practice only 5% of the network proteins which are on average homologous with a query protein.

Algorithm 4 may be summarized as follow. Perform ln( 1 ǫ )e k random colorations of the PPI network G N to ensure optimality with a probability of at least 1ǫ. The coloration consists in defining a function c : V N → S C , where the color in S C is randomly chosen. Then, for each coloration, we build all possible valid assignments A of the duplicated nodes. An assignment A is valid if no two non homologous nodes are matched in A. For each such assignment A, we compute the best score of an alignment according to A with Algorithm BESTCONSTRAINTALIGNMENT. We keep the best score of these trials and obtain the corresponding alignment by a classic backtracking technique. The score of the assignment of the duplicated nodes is computed separatly as follows. Indeed, in order not to take into account the homology score of q more than once (i.e. |d(q)| times) in the overall score (namely ScoreG A ), we precompute the part of the score induced by the duplicated nodes (namely score(A)).

score(A)

= q∈F A(q) =del h(q, A(q)) + q∈F A(q)=del δ d
Let us now describe in details the BESTCONSTRAINTALIGNMENT step that returns the score of G A according to the precomputed assignment A.

ScoreG A ← max v∈V N W M (root, v, S, 1, A) + score(A)
The best alignment score is obtained by finding among all possibilities the best way to align the root of the query to any protein v of the network. Similarly to QNet, the root is selected arbitrarily, but it is always a node of degree one. Moreover, the score is computed only if the root and v are homologous. In this initial step, S represents the multiset of colors defined previously. As in QNet, for each query node q, let denote by q 1 , q 2 , . . . , q nq its n q children. To obtain the best alignment, we use three tables, namely W M , W I and W D , which are filled as follows.

If |S| ≤ 1, W M (q, v, S, j, A) ← -∞ Else, W M (q, v, S, j, A) ← max u:(u,v)∈E N S ′ ⊂S c(v)∈S ′ c(u)∈S-S ′ W M (q, v, S ′ , j -1, A) +                      (* Matching, child j *) W M (q j , u, S -S ′ , n q j , A) + w(u, v), (* Insertion, node u *) W I (q j , u, S -S ′ , A) + w(u, v), (* Deletion, child j *) W D (q j , v, S -S ′ , A)
In the computation of W M (q, v, S, j, A), we consider that q is already aligned with v. Thus, the value stored in W M (q, v, S, j, A) is the maximum score of the subtree rooted at q and considering only its j first children. This score corresponds to the sum of the score of the subtree rooted at q considering only its j -1 first children (i.e., W M (q, v, S ′ , j -1, A)) and the score for the best aligment of the j th child of q -denoted q j . Indeed, when considering q j one can either (1) match it with a neighbor u of v, (2) delete it or (3) insert a neighbor u of v in the alignment. In the dynamic programming equation, we denote by n q j the number of children of q j .

To obtain the optimal solution, each subset of the multiset S of colors has to be considered. In others words, we consider each subset of colors used for the first j -1 subtrees of q (the subset S ′ ), and therefore each subset of colors used for the j th subtree (the subset S -S ′ ). In case S is a singleton, the score is -∞ since there are at least j + 1 nodes (i.e., q and the nodes of the j subtrees rooted at q 1 , q 2 , . . . q j ) to add in the solution with only one color, while the deletion of nodes with degree one is forbidden.

W M (q, v, S, 0, A) ←    h(q, v) if |d(q)| = 1, -∞ if A(q) = v, 0 else (* |d(q)| > 1 and A(q) = v *)
Entries corresponding to W M (q, v, S, 0, A) are the base cases of the recursion. When q is not a duplicated node (i.e., |d(q)| = 1), the value is simply the similarity score with v given by h. Otherwise, the assignment of q has already been defined in A and taken into account in score(A), and has to be preserved; otherwise, we return -∞ to forbid such an alignment.

W I (q, v, S, A) ←            if A -1 (v) = ∅ and |S| > 1 max u:(u,v)∈E N c(u)∈S-{c(v)} W M (q, u, S -{c(v)}, n q , A) + w(u, v) + δ i , W I (q, u, S -{c(v)}, A) + w(u, v) + δ i , else -∞
When computing W I (q, v, S, A), the node v is considered to be inserted. Thus, the computation will continue with a neighbor u of v, performing either a match or another insertion.

In order to be coherent with A, if v has to be aligned with a duplicate node (i.e., A -1 (v) = ∅), inserting v is forbidden.

W D (q, v, S, A) ←                                          if degree(q) = 2 -∞ else if |d(q)| = 1 max u:(u,v)∈E N    W M (q 1 , u, S, n q 1 , A) + w(u, v) + δ d , W I (q 1 , u, S, A) + w(u, v) + δ d , W D (q 1 , v, S, A) + δ d else if A(q) = del max u:(u,v)∈E N    W M (q 1 , u, S, n q 1 , A) + w(u, v), W I (q 1 , u, S, A) + w(u, v), W D (q 1 , v, S, A) else -∞
Finally, when computing W D (q, v, S, A), the node v and the father of q are considered as aligned and the node q to be deleted. The alignment then continues from q 1 , the only son of q since deletion is only allowed for nodes with degree two. If q is a duplicated node (i.e., |d(q)| > 1), the deletion is only allowed when A(q) = del to respect the reference assignment A. Again, in this case, the cost of a deletion δ d is already counted in score(A).

Let us note that since in the dynamic programming, for each j s.t. 1 ≤ j ≤ n q , all set S ′ ⊂ S are considered, hence the result does not depend on the ordering of T Q .

Let us now analyze the complexity of PADA1. The whole complexity depends essentially on lines 5 to 12. Let us consider the complexity of one iteration (we have ln( 1ǫ )e k iterations). The random coloration can be done in O(n), where n = |V N |. There are n |F | possible assignments in the worst case (i.e., if all the proteins in F are homologous with the n proteins of the network). The complexity of BESTCONSTRAINTALIGNMENT is 2 O(k+N ins ) mN del as in QNet, where k is the size of the graph query and m = |E N |, since our modifications are essentially additional tests which can be done in constant time.

Let us note that the complexity of GRAPH2TREE is negligible compared to the overall complexity of Algorithm PADA1. Indeed, the complexity of Algorithm GRAPH2TREE only depends on the query size k, with k ≪ n. Therefore, on the whole, the complexity of PADA1 is O(n |F | 2 O(k+N ins ) mN del ln( 1 ǫ )) time for any desired success probability 1ǫ (with ǫ > 0). Observe that the time complexity does not depends on the total number of duplicated nodes (i.e., q∈F |d(q)|), but on the size of F .

III. EXPERIMENTAL RESULTS

According to the authors of QNet, one may query a PPI network by running an 2 O(k) n t+1 time algorithm ln( 1ǫ )e k times, where t is the treewidth of the query graph. Thus, the difference between the two algorithms is mainly related to the "t + 1 versus |F |" question (where |F | is the size of the set of families of duplicated nodes computed by Algorithm GRAPH2TREE). Let us recall that the size of |F | given by our algorithm is equal to the feedback vertex set size of the query graph. According to Bodlaender and Koster [START_REF] Bodlaender | Combinatorial optimization on graphs of bounded treewidth[END_REF], the treewidth of a graph is at most equal to the feedback vertex set plus one. However, it only gives an upper-bound. We have conducted some experimental tests to compare these two parameters in practice on random graphs. The method is as follows: for each different size of graph (the number of nodes varies from 4 to 16 while the number of edges varies from |V | to 2.|V |), we get the average treewidth and feedback vertex set values over 30000 connected graphs, randomly constructed with the NetworkX library (http://networkx.lanl.gov/). Treewidth is computed with the exact algorithm provided by http://www.treewidth.com/, while the size of the feedback vertex set is computed with our GRAPH2TREE algorithm. Results in Figure 3 suggest that parameter |F | is usually smaller for moderate size graphs (i.e., those query graphs for which PADA1 and QNET are still practicable). In summary, PADA1 is an alternative to the QNet algorithm, with a different parameter involved in the time complexity. One can determine which parameter is more suitable for a given query and then uses the appropriate algorithm.

In practice, our upper-bound is largely over estimated. Indeed, each element of F must be assigned to a different node of the network, and hence, there are less than n possibilities for each element of F . The worst number of runs of BESTCONSTRAINTALIGNMENT is n! (n-|F |)! , the number of combinations. Moreover, we only consider valid assignments and there are only few such assignments. Indeed, a protein is, on average, homologous to dozens of proteins, which is quite less than the number of proteins in a classical PPI network (e.g. n ≃ 5.000 for the yeast). For example, if |F | = 3 and if the protein represented by this unique element of F is homologous to ten proteins in the PPI network, then, the number of assignment will not be n 3 but only 10 3 . Here, the running time is largely reduced. Therefore, and not surprisingly, the BLAST threshold used to determine if a protein is homologous to another has a huge impact on the running time of the algorithm.

Finally, observe that in QNet, for a given treewidth, the query graph can be very different. For example, in the resulting tree decomposition of the graph, there is no limit on the number of bags of size t. Furthermore, in a given bag, the topology is arbitrary (e.g., a clique), potentially requiring an exhaustive enumeration upper-bounded by n t+1 . Therefore, the treewidth value does not indicate how many times an exhaustive enumeration has to be done.

We would have liked to compare in practice our algorithm to QNet, but, unfortunately, their version querying graphs is not yet implemented. Comparing our algorithm for simple trees queries with QNet would not make sense since PADA1 is not optimized for this special case. [START_REF] Kanehisa | The KEGG resource for deciphering the genome[END_REF], with three cycles. b) The alignment graph given by our algorithm in the fly PPI network. Dashed lines denotes the BLAST homology scores between the two proteins. Our algorithm retrieves a query graph in an other network. As in QNet [START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF], it seems to be that there is some conservation between these two species.

In order to validate our algorithm, we perform the experimental tests on real data proposed by QNet [START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF]. In our experiments, the data for the PPI network of the fly and the yeast have been obtained from the DIP database1 [START_REF] Xenarios | DIP, the Database of Interacting Proteins: a research tool for studying cellular networks of protein interactions[END_REF]. The yeast network contains 4 738 proteins and 15 147 interactions, whereas the fly network contains 7 481 proteins and 26 201 interactions.

The first experiment consists in retrieving trees. To do so, the authors of QNet extract randomly trees queries of size 5 to 9 from the yeast network and try to retrieve them in this network. Each query is modified with at most two insertions or deletions. We also have successfully retrieved these queries.

The second experiment was performed across species. The Mitogen-Activated Protein Kinase (MAPK) are a collection of signal transduction queries. According to [START_REF] Dent | MAPK pathways in radiation responses[END_REF], they have a critical function in the cellular response to extracellular stimuli. They are known to be conserved through different species. We obtained the human MAPK from the KEGG database [START_REF] Kanehisa | The KEGG resource for deciphering the genome[END_REF] and tried to retrieve them in the fly network as done in QNet. While QNet uses only trees, we were able to query graphs. The results were satisfying since we retrieved them, with few or without modifications. The Figure 4 shows a sample of our results on real data. This suggests a potential conservation of patterns across species.

The BLAST threshold have deep impact on the running time. Moreover, we could certainly speed-up the running time by stopping earlier the dynamic programming. Indeed, one can stop if there are more nodes to color than the number of available colors in a step of the dynamic programming. This trick implies to look for the available number of deletions and insertions. Thus, for the moment we only stop when there is only a single color available. Another improvement can be to switch the coloration step with choices of assignements A. A final speed-up possibility can be to use the Hüffner et al. technique [START_REF] Huffner | Algorithm Engineering For Color-Coding To Facilitate Signaling Pathway Detection[END_REF], which basically consists in increasing the number of colors used during the coloration step.

IV. CONCLUSION In this paper, we have tried to improve our understanding in PPI networks by developing a tool called PADA1 (available uppon request), to query graphs in PPI networks. This algorithm, which is FPT for query graphs with a constant size of their feedback vertex set, has a time complexity of O(n |F | 2 O(k+N ins ) mN del ln( 1 ǫ )), where n (resp. m) is the number of nodes (resp. edges) in the PPI network, k is the number of nodes in the query, N ins (resp. N del ) is the maximum number of insertions (resp. deletions) allowed, ǫ is any value > 0 such that 1ǫ is the desired success probability, and |F | is the minimum number of nodes which have to be duplicated to transform the query graph into a tree (solving the FEEDBACK VERTEX SET problem). This last parameter is the main difference with QNet of Dost et al. [START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF], which uses the treewidth of the query (unimplemented algorithm). Consequently, PADA1 is an alternative to QNet and one can determine which parameter is better considering a query. These algorithms both use the color coding technique and are both exact for a given coloration. We have performed some tests on real data and have retrieved known paths in the yeast PPI network. Moreover, we have retrieved known human paths in the fly PPI network.

The time complexity of our algorithm depends on the number of nodes which have to be duplicated in the graph query. This number is directly connected to the initial topology of the query graph. Obtaining more information about the topology of the queries and the average number of homologous for proteins in the query are of particular interest in this context. Future works include using these informations to predict average time complexity.

Knowing if GRAPH TOPOLOGICAL CONTAINMENT -determining if a graph G has a subgraph that is a subdivision of a parameter graph H -is FPT or W [START_REF] Blin | Querying Protein-Protein Interaction Networks[END_REF]-hard is still an open problem (conjectured to be FPT by Fellows [START_REF] Fellows | Parameterized complexity: new developments and research frontiers[END_REF]). In the context of proteins queries, one can ask if the query H appears in the network G with an unbounded number of insertions.

Finally, recently, a problem close to this one, called GRAPH MOTIF, has been settled by Lacroix, Fernandes and Sagot [START_REF] Lacroix | Motif search in graphs: application to metabolic networks[END_REF]. Roughly speaking, one is only concerned in find a connected occurence of the query, which is defined without a given topology, that is a set or a multiset of proteins. This new definition leads to new issues, and has already been investigated [START_REF] Fellows | Sharp tractability borderlines for finding connected motifs in vertex-colored graphs[END_REF], [START_REF] Betzler | Parameterized algorithms and hardness results for some graph motif problems[END_REF], [START_REF] Dondi | Weak pattern matching in colored graphs: Minimizing the number of connected components[END_REF], [START_REF]Maximum motif problem in vertex-colored graphs[END_REF], [START_REF] Blin | GraMoFoNe: a cytoscape plugin for querying motifs without topology in protein-protein interactions networks[END_REF], [START_REF] Schbath | Assessing the exceptionality of coloured motifs in networks[END_REF], [START_REF] Bruckner | Topology-free querying of protein interaction networks[END_REF].

4 while ∃ a cycle C in G do 5 7 DUPLICATE

 457 Let u be one vertex randomly chosen in C ∩ F V S;6 Let v be one vertex randomly chosen in C ∩ N (u);

10 Algorithm 2 :

 102 Compute the Feedback Vertex Set of a graph.

8 end 9 Algorithm 3 :

 893 Algorithm to duplicate a node when a cycle is detected.

Figure 1 .

 1 Figure 1. Steps when Duplicate(G, v3, v1, d) is called on graph a). b) A node v1,1 from v1 is created. c) The edge (v3,v1) is deleted: the cycle is then broken. d) The edge (v3,v1,1) is added. Finally, the resulting graph is acyclic, and d(v1) = {v1, v1,1}.

Figure 2 .

 2 Figure 2. a) The graph query with a cycle, before calling GRAPH2TREE algorithm. c) The query after calling GRAPH2TREE where q1 has been duplicated. Thus, q1 and q1,1 have to be aligned with the same node of the network. b) and d) denote the resulting graph alignment GA, subgraph of the network GN . The horizontal dashed lines denote a match between two proteins.

3 for 5 foreach valid assignment A do 6 ScoreG 7 if ScoreG A > BestScore then 8 Save coloration, A; 9 BestScore ← ScoreG A ; 16 Algorithm 4 :

 356789164 (i = 0; i < ln( 1 ǫ )e k ; i + +) do 4 randomly colorize G N with k + N ins colors; A ← BESTCONSTRAINTALIGNMENT(G N , T Q , A, h, threshold) + score(A); Sketch of the PADA1 algorithm to align a query graph to a network.

Figure 3 .

 3 Figure 3. Comparison between QNet (i.e., the treewidth+1 value) and PADA1 parameter (i.e., the size of the feedback vertex set). In QNet (resp. in PADA1), the treewidth (resp. the feedback vertex set) is computed over the query. Here the size of the query graph correspond to the number of nodes (there are usually between five and fifteen proteins in a classical query).

Figure 4 .

 4 Figure 4. A result sample of our algorithm. a) A MAPK human query, get from[START_REF] Kanehisa | The KEGG resource for deciphering the genome[END_REF], with three cycles. b) The alignment graph given by our algorithm in the fly PPI network. Dashed lines denotes the BLAST homology scores between the two proteins. Our algorithm retrieves a query graph in an other network. As in QNet[START_REF] Dost | QNet: A Tool for Querying Protein Interaction Networks[END_REF], it seems to be that there is some conservation between these two species.

http://dip.doe-mbi.ucla.edu/

ACKNOWLEDGEMENT

The authors would like to thank Banu Dost for providing us the QNet source code and their test data. We also thank the anonymous reviewers for their helpful comments and suggestions for improving the manuscript.

Stéphane Vialette a Researcher in the LIGM at Université Paris-Est -France. He received his PhD degree in 2001 from the University Paris 7 -France. His research interests are in computational biology, algorithmics and combinatorics.