Victor Epain

Rumen Andonov

Dominique Lavenier

Inverted Repeats Scaffolding for a Dedicated Chloroplast Genome Assembler

Keywords: 2012 ACM Subject Classification Applied computing Genome assembly, Inverted repeats, Contig, Scaffolding, De Bruijn graph, Assembly graph, Elementary longest path problem, Integer programming, NP-Complete problem Digital Object Identifier 10.4230/LIPIcs

This paper describes a novel assembly approach for chloroplast genomes. It contains two modular steps. In the first step, based on the hypothesis that chloroplasts genomes are over-represented compared to the nuclear genome in the plant's cell, we assemble contigs with a De Bruijn graph based approach using short reads with a high k-mer coverage. Connections between oriented contigs are also provided here.

The second step determines the order and the orientation of the contigs (scaffolding). Taking advantage of the knowledge that chloroplast genomes posses well studied circular structure, we develop a particular formulation of the scaffolding problem, called Nested Inverted Fragments Scaffolding. It aims to assemble highly conserved inverted repeats. We formulate it as an optimisation problem and we prove that it is NP-Complete. To solve the problem we propose and implement an integer linear programming formulation.

We evaluate our method on a set of real instances (a benchmark of 42 chloroplast genomes) and show that it obtains notable achievements with respect to the quality of the results. To further estimate the performance of our scaffolding module, we test it on huge artificially created instances. The results demonstrate an excellent behaviour of our integer formulation as even very large instances have been solved at the first Branch & Bounds node.

Introduction

Sequencing an organism consists in extracting DNA molecules contained in its cells and getting them as sequences of ATGC letters corresponding to nucleotides. Modern sequencing technologies are still unable to return one complete sequence for each DNA molecule but instead yield fragmented substrings of them called reads that need to be assembled to reconstruct the original genome. A major difficulty for the assembly is provoked by genomes repeats: a read from one region can overlap another read containing a repeat of the first one, while both reads are sequenced from distanced genomic locations. Moreover, reads are sequenced without distinction from the two complemented DNA strands. Thus, Inverted repeats (a DNA sequence is the reverse-complement of another one) can produce repeatinduced overlaps and can wrongly mix regions from both DNA strands.

XX:2

Inverted Repeats Scaffolding

This paper presents an approach for assembling the chloroplast genome, a plant organelle responsible for photosynthesis metabolism, which confers the ability to transform sunlight energy to chemical energy. The knowledge of chloroplast genomes allows evolutionary analyses [START_REF] Xiao-Ming | Inferring the evolutionary mechanism of the chloroplast genome size by comparing whole-chloroplast genome sequences in seed plants[END_REF], (meta)-barcoding [START_REF] De Vere | DNA Barcoding for Plants[END_REF], and is useful for biophotovoltaic process development [START_REF] Tschörtner | Biophotovoltaics: Green Power Generation From Sunlight and Water[END_REF].

Several specialised assemblers for chloroplast genomes have been already developed in [START_REF] Dierckxsens | NOVOPlasty: de novo assembly of organelle genomes from whole genome data[END_REF][START_REF] Jin | GetOrganelle: a fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] and have been compared in [START_REF] Freudenthal | A systematic comparison of chloroplast genome assembly tools[END_REF]. However, none of them explicitly formulates the specific subproblems that need to be solved. As a consequence, they apply essentially (meta)heuristics.

The chloroplasts possess a circular genome which is characterized by very specific but well-described repeated regions [START_REF]Genomics of Chloroplasts and Mitochondria[END_REF]. An instance of such a structure is given in Subfigure 1a. Deciphering these particular genome repeats is the first requirement for a successful chloroplast assembler. One more specific issue here is the separation of the plant reads from the chloroplast ones. Indeed, splitting them is not an easy task as there can be overlaps between both groups. As consequence, choosing repeat-induced or plant-chloroplast overlaps can lead to a misassembled genome.

Our approach can be cut in two main parts. The first part generates contigs that are expected to be part of the chloroplast. A multiplicity is attached to each contig. It provides an upper bound of the number of copies per contig. A list of potential links between contigs is also generated (c.f. Table 1 for an instance of input data: set of contigs and set of links). Finally, one contig is tagged as seed contig using proteome information from a close plant species. It is supposed to contain one reference gene present in a single copy.

In the second stage called scaffolding we finish the assembly, i.e. we order and orient the contigs yielded by the first step. The particularity of our approach is that it is based exclusively on the knowledge that most chloroplast genomes are circular and has two inverted repeats [START_REF]Genomics of Chloroplasts and Mitochondria[END_REF] (see for illustration Subfigure 1a). In contrast to some of the previous chloroplast genome scaffolders [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF], we need no additional distance-constraints encoding the insert-size information.

This step is modular and completely independent from the first one. Any assembly graph obtained by another assembler can serve here as an input data. The data are firstly represented in a structure that we call Multiplied Doubled Contigs Graph (MDCG). We search in this graph to finish the assembly by reconstructing inverted repeats using the observation that they represent nested couples of inverted regions. Thanks to a mathematical specification on the inverted repeats structure, we maximally reduce the assembly graph. Thus, we model the chloroplast scaffolding as building a sequence of oriented contigs, starting from and ending with the seed contig, which maximizes the number of nested and contiguous inverted fragments. We call the underlying optimisation problem Nested Inverted Fragments Scaffolding Problem (NIFSP) and prove that it is NP-Complete. We propose an Integer Linear Programming (ILP) formulation and show that the number of variables and constraints is polynomial in respect to the vertices and the edges number. Although our formulation is somehow similar to the one proposed in [START_REF] Gusfield | The RNA-Folding Problem[END_REF] for RNA folding, here we deal with contigs instead of nucleotide sequences, without any knowledge of their order, which obviously increases the underlying challenges. We use the Gurobi software [START_REF] Optimization | Gurobi Optimizer Reference Manual[END_REF] to solve the optimization problem, resulting in a small number of larger contigs that correspond to chloroplast's genomic regions. Finally, our tool yields a scaffolding enriched assembly graph, as illustrated Subfigure 1b. We perform numerical computations and evaluate the quality of the results using Quast evaluation tool [START_REF] Gurevich | QUAST: quality assessment tool for genome assemblies[END_REF]. For the contig indicated as u f,0 on the left-side green region, the contig in front of it (on the right-side green region) has the reverse orientation (so r). Because orientations are mutually exclusive, it is necessary to chose an occurrence that differs 0 (so 1). The couple u f,0 , ur,1 is defined as inverted fragments. Each inverted fragments couple is linked by blue dashed line. Inverted repeats can be modelled as sequences of nested inverted fragments. (b) Vertices correspond to genomic regions while edges indicate their neighbourhoods. Purple, green and red vertices correlate respectively with purple, green and red regions in Subfigure 1a. Two green vertices' label differ by their orientation 'f ' and 'r' (one is the reverse complement sequence of the other). They also differ by their occurrence number (0 and 1). It simply means that their associated genomic regions are simultaneously present in the genome sequence (i.e. the green vertices both participate in the solution). On the opposite, red vertices just differ by their orientation: it implies that only one genomic region from one of these two vertices appears in the genome (one genomic region is the reverse complement of the other).

2

Material & Method

Material & Method Overview

Input data are reads that come from short reads sequencing technology, and can be part of plant's or/and chloroplasts' genomes. We also use chloroplast proteome data from a near plant species, and a well known gene sequence that is contained with a very high probability in one of the chloroplast's unique genomic regions, defined as the seed gene. An overview of the method is illustrated in Figure 2. In the next section we briefly describe how chloroplast contigs with their links are extracted.

Data Filtering & Contigs and Links Generation

Our strategy supposes that the chloroplast genome is over-represented compared to the nuclear genome of the plant. We first perform a de Bruijn graph (DBG) assembly with Minia [3] from a subset of the sequencing dataset. Typically 3 millions of reads are randomly selected (this number can be adjusted by the user depending on the chloroplast coverage). A coverage is computed for each resulting contig. Contigs are then aligned to the chloroplast proteome of a close species. Those that present significant similarity are tagged as near-tothe-proteome contigs. The one that presents the highest similarity with a predefined gene (which is known to be present in a single copy), typically the matK gene, is set as the seed contig. Its coverage becomes the reference coverage. From these set of contigs, we then search for links between them. This is done by extending all the near-to-the-proteome contigs using overlapping read strategy. Contigs, Figure 2 Method overview.The method can be cut in two main parts. First step is a driven chloroplast assembly with reads that come from both plant's genome and chloroplasts' genome. This assembly is a De Bruijn graph assembly approach. It outputs contigs' successions and an estimated multiplicity for each contig. These results are structured in a Multiplied Doubled Contigs Graph (MDCG). The following scaffolding step is done thanks to an Integer Linear Programming (ILP) optimisation. Finally, we output a fewer number of larger contigs that correspond to chloroplast's genomic regions, and a scaffolding enriched assembly graph, as illustrated Subfigure 1b.

other than seed contigs, having connection with seed contigs are fished out and added to the contig set.

We define the contigs set C as the set containing near-to-the-proteome tagged contigs and those which are intermediate in paths between two tagged ones. For each contig we compute a multiplicity value as an approximated ratio between its estimated coverage and the seed contig's coverage (the reference coverage). Thus, a contig's multiplicity may be interpreted as an upper bound of the number of times the associated contig can be met in the final assembly.

To any contig is associated a nucleotide sequence which length is provided as input data. A contig can participate in the solution either in its forward orientation (so the sequence does not change), or in its reverse orientation (the sequence is read in inverse reading and each nucleotide is complemented). However, the orientation of a contig participating in the solution is unknown a priory.

Let N denote the set of natural numbers while N * = N \ {0}. As usual, B denotes the set of binaries.

Any contig c ∈ C is characterized by its identifier c id ∈ N and its multiplicity c mult ∈ N * (an upper bound of the number of occurrences of c). Contigs and their multiplicity are illustrated in Subtable 1a.

In the next step we to find paths between tagged contigs in the DBG. These paths are defined as links and we denote the corresponding set by L. When such a link traverses a limited number of non-tagged contigs, these contigs are selected to be in the final contigs set. Therefore, the links are succession relations between a couple of oriented contigs, i.e. L ⊂ (C × B) 2 . Each contig c is provided here with an orientation attribute orc or ∈ B (thus c becomes orc), where the binary value 0 stands for forward (f ), i.e. the original contig's sequence, while the binary value 1 stands for reverse (r), i.e. the reverse complement of the contig's sequence. Subtable 1b gives an illustration for links data.

Multiplied Doubled Contigs Graph (MDCG)

Building chloroplast genomic regions implies finding a sequence of linked oriented contigs, in the limit of multiplicity for the two orientations for each contig. Links between two multiplied (by their multiplicity) and oriented (forward or reverse) contigs can be structured in an Set C of contigs with their multiplicity (c mult ); (b) Set of links L. An identifier (orc id ) and an orientation (orcor) is provided here for each oriented contig orc serving as a link extremity. For the sake of saving memory, for each link (orc, ord) ∈ L, only one of the existing two links is reported in the table, the one with orc id < ord id . oriented graph. Thus, we define the Multiplied Doubled Contigs Graph, M DCG = (V, E), a directed graph where V is the set of vertices and E is the set of edges1 . Even if in practice we do not multiply and double contigs and links data in memory, each vertex v ∈ V is one occurrence of an oriented contig and each edge (u, v) ∈ E corresponds to a link between occurrences of two oriented contigs. The construction rules and some MDCG properties are as follows:

(i) Each contig c ∈ C is doubled according forward and reverse orientation. So ∀v ∈ V such as v id = c id , v has an orientation attribute v or equals to 0 or 1 respectively for forward or reverse orientation. We denote by v the reverse orientation of v (so that v or = 1 -v or , and the other attributes remain the same). For each vertex v ∈ V , its reverse v is in the graph too. By definition, vertex's identifier remain the same v id = v id , and likewise for vertex's occurrences, v occ = v occ . (ii) Moreover, each doubled contig is multiplied by its multiplicity. Thus, if contig c is multiplied k times, it generates a set of 2k vertices v 0 , . . . , v k-1 , v 0 , . . . , v k-1 where all these vertices have the same identifier c id . Let 0 ≤ v occ < c mult be an occurrence of a multiplied contig c. (iii) The set E of edges is constructed based on the links data L and applying the same multiplication rules as for the contigs. Furthermore, for any edge in E, its reverse is also added in E in order to validate the feature (u, v)

∈ E ⇐⇒ (u, v) = (v, u) ∈ E.
Furthermore, we add into MDCG the data that we define as inverted fragments and denote by InvF (its formal definition is given in section 2.5.2). Each inverted fragments is a non-oriented couple of vertices (i, j), such that they are two different occurrences of the same contig, but one is in forward orientation, while the other is in reverse orientation. These couples are candidates to be part of inverted repeats pairs.

Nested Inverted Fragments Scaffolding Problem

Thanks to the MDCG structure, building a sequence of linked oriented contigs is equivalent to finding an elementary path in MDCG. This is a necessary, but not a sufficient condition to correctly retrieve a circular genome with pairs of inverted repeats. Especially, no more that one of the vertices v ∈ V and its reverse v ∈ V can be in the path. The seed contig is illustrated as a big black dot in Subfigure 1a & 3 (seed vertex s ∈ V , in arbitrarily forward orientation). The circularity of the genome corresponds to a path that begins at s and ends in s.

Finding a path requires giving a position for each vertex participating in it. As illustrated in Subfigure 1a, two inverted fragments (i, j) and (k, l) in inverted repeats regions are nested i.e. if you draw in the corresponding circuit a line from i to j, and another from k to l, then the two lines do not intersect. In order to know if two inverted fragments intersect, we must compare the positions of the associated vertices.

Problem definition

Finally, the goal of the Nested Inverted Fragments Scaffolding problem (NIFSP) is to find in the corresponding MDCG an elementary path from s to s, that includes at most one of any of the vertices v ∈ V or its reverse v ∈ V , and that traverses as many as possible inverted fragments in the way that they are located on the path as contiguous nested pairs.

An illustration for such a path is given in Figure 3. To solve the NIFSP we describe below a linear programming formulation where vertices' position, their relative location, chosen edges and chosen inverted fragments are all integer or binary variables.

Modelling the Scaffolding Problem as a Path Search in a Graph: Integer Linear Programming Formulation

In the sequel, let N - v and N + v denote respectively the sets of predecessors and successors of a vertex v ∈ V .

In order to define a path in a graph, one usually needs to designate its start and its terminal vertex. We use here for this purpose the given seed contig. We denote its corresponding vertex by sc and we suppose that it is in its forward orientation. We also know that it is present in the graph with only one occurrence. Furthermore, we replace sc by two new vertices s and t, both in forward orientation and with only one occurrence, where s gets all outgoing sc edges, while t gets all incoming sc edges. Note that we add neither s nor t and that N - s = N + t = ∅, where ∅ denotes the empty set.

Path Definition

Equations ( 1)- [START_REF] François | Global Optimization for Scaffolding and Completing Genome Assemblies[END_REF] in this section follow the single-flow formulation of an s -t path that is used in the field [START_REF] Taccari | Integer programming formulations for the elementary shortest path problem[END_REF] and that we have adapted to the scaffolding problem specificities in our previous papers [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF][START_REF] François | Global Optimization for Scaffolding and Completing Genome Assemblies[END_REF]. Vertices s and t will be used as the source (start) and the target (terminal), respectively, of the path we are looking for. We introduce a binary variables x e to indicate if the corresponding edge e participates or not in the path:

∀e ∈ E : x e ∈ B. (1)
An unique path exits from s and enters in t, i.e.:

v∈N + s x sv = u∈N - t x ut = 1 (2) With any vertex v ∈ V \ {s, t} is associated a variable i v ∈ [0, 1] to encode whether v is, or isn't in the solution path. If v /
∈ {s, t} is an intermediate vertex in the path, then v possesses exactly one incoming/outgoing edge. If v is not in the path, then no edge inputs/exits v, i.e.

∀v ∈ V \ {s, t} : i v = u∈N - v x uv = w∈N + v x vw (3)
Note that although i v is declared as contiguous variable, it takes 0/1 value because of (3). This is a particular case of a more general theorem proven in [START_REF] François | Global Optimization for Scaffolding and Completing Genome Assemblies[END_REF].

Furthermore, for each vertex v, at most one of its orientations participates in the path:

∀v ∈ V \ {s, t} : i v + i v ≤ 1 (4)
We introduce a variable f e ∈ R to express the quantity of the flow circulating along the edge e ∈ E. No flow can use an edge e when x e = 0, i.e.

∀e ∈ E :

x e ≤ f e ≤ x e × M ( 5
)
where M = |V |. For the initial/intermediate vertex flow we set respectively:

v∈N + s f sv = 1 (6)
∀v ∈ V \ {s, t} :

w∈N + v f vw - u∈N - v f uv = i v (7)
Constraint ( 7) makes the flow increasing and forbids sub-tours. On the other hand, the exit flow value can be interpreted as a position/number of a vertex v on the genome in respect to the start s. Accordingly, we denote this value by pos v , i.e.

∀v ∈ V \ {t} : pos v = w∈N + v f vw (8)
Constraints ( 4) -( 7) define an elementary path from s to t. In the sequel, we add to this model non-intersecting inverted fragments constraints.

∈ N + i ⇐⇒ j ∈ N + l .
We aim to maximise the number of such adjacent inverted fragments in the solution. We do not consider dashed edges in AInvF because their extremities participate in no one of the InvF 's item. Note that edge (i, k) ∈ E (in bold) is a canonical -edge (l, j) ∈ E can be retrieved by adjrev((i, k)) = (l, j).

Inverted Fragments

Let R be the set of contigs with multiplicity greater than 1, i.e. R = {c ∈ C | c mult > 1}. An oriented couple consisting of a forward occurrence of a contig from R, and another occurrence of it, but in reverse orientation is called inverted fragment (c.f. Subfigure 4a for illustration). Denote InvF the set of inverted fragments:

InvF =    (i, j) ∈ V 2 , ∀c ∈ R (i id = j id = c id ) ∧ (i or = 1 -j or = 0) ∧ i occ = j occ -1 = 2 × k, 0 ≤ k < c mult 2   
Here the first set definition determines inverted fragments, while the second one describes the minimum number of items to build all distinct solutions. It is also necessary to describe a set of pairs of inverted fragments:

P InvF = (i, j), (k, l) ∈ InvF 2 (j id < k id ) ∨ (j id = k id ∧ j occ < k occ ) .
The condition on the identifier and the occurrences attributes assures that for each pair (p, q) ∈ P InvF , its permutation pair (q, p) / ∈ P InvF , and generates the minimum number of items for building all distinct solutions (c.f. Subfigure 4a for illustration). The set P InvF contains unordered pairs of inverted fragments that must be checked for admissibility in respect to constraint [START_REF] Jin | GetOrganelle: a fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF].

Let (i, j), (k, l) ∈ P InvF and consider the positions of these inverted fragments in the genome. In order to illustrate the options for their relative locations, and without loss of generality, we assume that pos i < pos j , pos k < pos l and that the interval pos i , pos j is smaller or equal to the interval pos k , pos l . Only one of the following three cases holds: pos i , pos j ∩ pos k , pos l = ∅ pos i , pos j ⊆ pos k , pos l pos i , pos j ∩ pos k , pos l = ∅ and pos i , pos j ⊆ pos k , pos l

The first two cases correspond to feasible solutions to our problem, while in the last one the inverted fragments intersect. The later case is not admissible and will be forbidden as shown below.

Denote by D alpha be the set of couples of vertices whose positions must be compared in order to detect intersection, and define it as follows:

D alpha = (i, k), (i, l), (j, k), (j, l) ∀ (i, j), (k, l) ∈ P InvF ∀(u, v) ∈ D alpha we introduce a binary variable α v
u to indicate the relative locations of these vertices. α v u equals 1 if pos u < pos v , otherwise 0. We formulate this definition as the below big-M constraint where M = |V |:

∀(u, v) ∈ D alpha : -α v u × M ≤ pos u -pos v ≤ (1 -α v u ) × M -1 (9) α v u ≤ pos u + pos v ( 10
)
Note that α v u = 1 -α u v , hence only one of the variables α v u and α u v needs to be maintained. Furthermore, let us consider two couples of inverted fragments (i, j), (k, l) ∈ P InvF . It is easy to check that amongst the 24 possible permutations between their respective positions, 8 correspond to intersections that are not feasible solutions and we need to exclude them from consideration. To any of these intersections we associate a binary variable. When the pair (i, j) is chosen to be the first argument, we introduce four binary variables inters n (ij, kl), n = 1, . . . , 4. For example, the case pos i < pos k < pos j < pos l will be associated with a binary variable inters 1 (ij, kl) that equals 1 if the previous three relations are true, 0 otherwise. This definition is related to (9) by the following constraints:

3 × inters 1 (ij, kl) ≤ α k i + (1 -α k j ) + α l j ≤ 2 + inters 1 (ij, kl) (11)
We need four additional binary variables inters n (kl, ij), n = 1, . . . , 4 to model the case when the pair (k, l) is the first argument.

Moreover, we use binary variables m(i, j), (i, j) ∈ InvF , to indicate if the inverted fragment (i, j) can be matched (chosen to be in the solution), or not. When the pairs (i, j) and (k, l) do intersect, no more than one of them can participate in the solution. We model this condition by the below constraint:

∀ (i, j), (k, l) ∈ P InvF m(i, j) + m(k, l) ≤ 2 - 4 n=1 inters n (ij, kl) + 4 n=1 inters n (kl, ij) (12)
The vertices of a chosen inverted fragment must belong to the path, i.e.

∀(u, v) ∈ InvF : m(u, v) ≤ i u ; m(u, v) ≤ i v . ( 13
)

Inverted Repeats Contiguity

Inverted fragments represent subsequences in the inverted repeats regions. In this section we study how to assemble them in a contiguous way, and how to make these assembled regions as long as possible. We achieve this goal by assembling adjacent inverted fragments. Those are fragments that are extremities of an edge in the considered MDCG. Since for any edge in E, its reverse is also in E, we keep only one of them (considered as canonical edge) in a set of edges called AInvF . However, only containing canonical edges implies to be be able to retrieve: the inverted fragments associated with any of the extremity of each (u, v) ∈ AInvF . This task is performed by the function invfrag

: {v ∈ V | v occ + (1 -v or ) < v mult } → InvF (c.
f. Subfigure 4a for illustration); the reverse edge for each (u, v) ∈ AInvF , task performed by the function adjrev :

AInvF → E (c.f. Subfigure 4b for illustration). We introduce a binary variable isadj(u, v), ∀(u, v) ∈ AInvF . It equals 1 if the two inverted fragments associated with the canonical edge (u, v) are contiguous in the solution, 0 otherwise. The contiguity of inverted fragments implies, ∀(u, v) ∈ AInvF : the canonical edge (u, v) and its adjacency reverse adjrev(u,v) must participate in the solution:

isadj(u, v) ≤ x uv ; isadj(u, v) ≤ x adjrev(u,v) (14)
the associated inverted fragments do not intersect:

isadj(u, v) ≤ m(invfrag(u)) ; isadj(u, v) ≤ m(invfrag(v)) (15)

Objective Function

The goal of the scaffolding problem considered here is to find an elementary s -t path in the M DCG that contains as many as possible contiguous non-intersecting inverted fragments. We obtain the below problem: Proof. By reduction from the longest path problem from vertex s to vertex t (LPSTP), known to be NP-Hard [START_REF] Schrijver | Combinatorial Optimization: Polyhedra and Efficiency[END_REF].

Consider an instance I ∈ LPSTP that is composed of a directed graph G = (V, E) and two vertices s ∈ V and t ∈ V . We shall build an instance transform function mdgf such that I ∈ LPSTP ⇐⇒ mdgf (I) ∈ NIFSP. As illustrated in Figure 5, function mdgf transforms the graph G from 5a to the graph G = (V , E ) in 5b, vertex s ∈ V to vertices s f ∈ V and t f ∈ V , while vertices i f and i r in 5b, are images of t in 5a.

All subgraphs G * in 5b are quasi-clones of G V \{s,t} in 5a. In the graphs

G f 0 = (V f 0 , E f 0 ) and G f 1 = (V f 1 , E f 1
) edges are oriented in the same direction as in G V \{s,t} . Vertices in V f 0 , V f 1 have forward orientation attribute, and respectively 0 and 1 as occurrence attribute. Vertices in V r0 , V r1 have reverse orientation, while their edges E r0 , E r1 are oriented in the opposite direction compared to E f 0 , E f 1 edges. The sets InvF , P InvF and AInvF are constructed based on G * graphs' data. We assume that inverted fragments are composed of vertices from V f 0 and V r1 uniquely, (i.e. ∀(u, v) ∈ InvF , u ∈ V f 0 , v ∈ V r1 ). This is depicted by blue dashed vertical lines from G f 0 to G r1 . Vertices i f and i r allow the transition from G f 0 and G f 1 to G r0 and G r1 .

It is straightforward to see that there exists a O(|V | + |E|) algorithm that computes this transform function.

As adjacent inverted fragments associate only vertices in V f 0 with those in V r1 , the path that maximizes the number of contiguous inverted fragments exits s f , goes through G f 0 to i f (or i r , it does not matter), and passes through G r1 to t f . Since G f 0 is a copy of G V \{s,t} , while G r1 is its reverse graph, there is a bijection between V f 0 and V r1 vertices sets. Hence maximising the number of contiguous inverted fragments is equivalent to finding the longest path

v f 0,1 • • • v f 0,n in G f 0 and its reverse v r1,n • • • v r1,1 in G r1
. These two paths have the same length, that equals the length of the longest path in G.

To conclude, as there is a linear time complexity transform function mdgf such that

I ∈ LPSTP ⇐⇒ mdgf (I) ∈ NIFSP, NIFSP is at least NP-Hard.
Proposition 2 (Solution verification is done in polynomial time). Given a path in MDCG and a set of chosen inverted fragments, there is an algorithm that verifies if these inverted fragments are nested along the path.

Proof. Let P ath be the path of oriented fragments in MDCG, and let InvF * be the set of chosen inverted fragments.

Algorithm 1 verifies if the solution is admissible. Note that the stack repeats has the following four methods: peek Returns the vertex on the top of the stack without removing it. If the stack is empty, it returns N ull.

add Adds the vertex on the top of the stack. remove Removes the vertex on the top of the stack. is_empty Returns T rue if and only if, the stack is empty.

Algorithm 1 is linear in the length of the path P ath.

From Propositions 1 and 2 we conclude that NIFSP is NP-Complete.

XX:12 Inverted Repeats Scaffolding

Algorithm 1 Given solution verification algorithm.

Require: M DCG, P ath and InvF * Ensure: Return T rue if, and only if the given solution is admissible.

1: repeats ← LIFO() repeats is a stack.

repeats.add(invfrag(v)[1 -v or ])

12: return repeats.is_empty()

Results

The first assembly step, from reads to contigs' sequences, multiplicities and successions was implemented in Python3. The described ILP formulation for the IR scaffolding was implemented in Python3 using the PuLP package where we run Gurobi solver with academic licence. All the instances have been executed on a Linux laptop computer (32GB RAM, Intel ® Core™ i7-10610U CPU @ 1.80GHz ×8).

Real Instances

Input Data

We tested the whole pipeline on reads that mix nuclear and chloroplasts genomes. All reads and reference genomes are the same as in the GetOrganelle paper [START_REF] Jin | GetOrganelle: a fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF]. A subset of 42 chloroplast genomes with inverted repeats and with a chloroplast coverage at least equal to 15 was chosen.

Quality of the Assemblies

The scaffolding module solved 37 instances out of 42. Further analyses revealed that all unsolved instances failed because of missing links. As a consequence, finding a circular path in MDCG was impossible. We were able to asses the quality of the solution concerning the 37 solved instances as the genomes for these instances are known. We used the well known assembly evaluation tool Quast [START_REF] Gurevich | QUAST: quality assessment tool for genome assemblies[END_REF]. The results corresponding to the main measures output by Quast are presented in Figure 6. Tables A.3-?? given in the appendix reports all details. Below we give our observations and explanations concerning the erroneous assemblies. All other genomes have been assembled almost perfectly (with a reference genome coverage higher than 99%, c.f. Subfigure 6b).

Aoysia citriodora genome reference is covered at 98.68%. The reason is that one contig is not connected to the main assembly graph component, and another one is missing. Boswellia sacra solution covers 98.51% of the reference because one contig is misassembled (the one corresponding to the long unique region). Inverted regions are successfully assembled. Dendrobium nobile solution (96.55%) has both contigs corresponding to the two unique regions misassembled, while the inverted repeats are successfully assembled. Haberlea rhodopensis solution covers 97.53% of the reference. Contig corresponding to IR region is shorter than it should be. It is probably due to underestimated multiplicities. Jasminum tortuosum solution (99.96%) has one misassembled contig corresponding to unique region. However, inverted repeats are successfully assembled. Pimenta dioica solution covers 98.01% the reference genome. Although inverted repeats are successfully assembled, the contig corresponding to LSC region is missing some subsequences.

Time of the Assemblies

For each instance, the first assembly step spends 3 minutes as mean time. Formatting its outputs to MDCG and writing the model with PuLP was less than one second. As NIFSP is proved NP-Complete, we measured the Gurobi solving time. Thus each instance has been run 10 times for the ILP solving step. Subfigure 6d shows the distribution of the LP relaxation and Branch & Bound (B&B) times sum.

Multimerism Revealed by Final Assembly Graph

For each feasible instance, the scaffolding module outputs 3 contigs: two of them correspond to two unique regions known as long and short single copy (LSC & SSC). The the third one corresponds to one of the inverted repeats region (IRa or IRb). In addition of these contigs, the final assembly graph that visualises the connections between the three genomic regions is provided. This assembly graph is enriched by reverse properties as illustrated in Subfigure 1b. Indeed, the graph in Subfigure 1b shows that there are two possible circular sequences of genomic regions: 0 -1 f,0 -2 f -1 r,1 and 0 -1 f,0 -2 r -1 r,1 . These two equivalent solutions illustrate multimerism phenomenon for chloroplast genomes, studied in [START_REF] Deng | The chloroplast genome exists in multimeric forms[END_REF]: two conformations can simultaneously exist in the same plant's cell. In fact, one of the unique genomic region between the two inverted repeats (e.g. the red region in Subfigure 1a & 1b) exists in forward orientation in some genomes, and in reverse orientation in others.

Artificial Data

In order to better analyse the limit and the behaviour of our linear model, we applied it to artificially generated data instances (contigs, links and multiplicities). Each instance has been run 5 times. Results are reported in Table 2. The average gap of the linear relaxation bound (U B) with respect to the optimal integer value (Opt), computed as 100 × (U B-Opt) U B (MIP gap), differs from zero for only one of the reported instances. Moreover, all the instances have been solved at the first B&B node. They confirm the quality of the integer linear formulation: the gap values stay extremely small even for the largest instance and all instances have been solved at the first B&B node.

Conclusion

In this paper we propose a dedicated chloroplast genome assembler. Although the first step is driven by a proteome from a close species, the second step is a scaffolding approach only based on chloroplast genomic knowledge. The first step provides contigs, an estimated multiplicity for each one and links between oriented contigs. Then the scaffolding strategy aims to built a special case of genomic regions: inverted repeats. Thus, it consists in finding an elementary circular path into MDCG that maximises the number of nested and contiguous inverted fragments. This optimisation problem defined as NIFSP is proven to be NP-Complete. We also design an integer programming formulation and apply the associated implementation to solve real genomic and artificial data instances.

Concerning real data instances, first results evaluated by Quast are very encouraging: the scaffolding module rarely outputs misassembled contigs and mismatches and indels errors rates is very low. Furthermore, the results on artificial data demonstrate that our formulation is very tight -the MIP gap value stays extremely small and we successfully solve large instances at the first B&B node.

Since chloroplast genomes are multimeric (regions between inverted repeats can be considered in their both orientations) this raises the question of their suitable representation. To answer this question our strategy here is to explicitly reveal this characteristic using an enriched assembly graph where each region corresponds to one node.

On the one hand, we propose here a scaffolding strategy for inverted repeats which do not suffer from any degeneration. However, a recent study suggests cases in which a gene can exist in a single occurrence instead of two whereas it belongs to inverted repeats [START_REF] Turmel | Divergent copies of the large inverted repeat in the chloroplast genomes of ulvophycean green algae[END_REF]. Therefore, one repeat is split into two parts, resulting in the loss of contiguity between the inverted fragments. Thus the contiguity constraint may appear restrictive in this context, as the solver may prefer to bypass the contig containing the gene in order to make contiguous the inverted fragments before and after it. We can easily solve this issue by removing the contiguity constraint.

On the other hand, not all chloroplast genomes posses inverted repeats structure. Some of them have a pair of directed repeats, others a combination of different repeat types. Our next step will be to mathematically formulate direct repeats and unique region reconstruction. The linear programming techniques that we have described in this article can be easily adapted to consider these additional cases.

A Detailed Results on Real Data

( a )Figure 1

 a1 Figure1Genomic knowledge and corresponding perfect assembly result. (a) Each arrow represents an oriented DNA sequence (contigs). Chloroplast's genome is circular thus it is a circular sequence of contigs, that begins at s and ends in s. Purple and red regions are unique, while green ones are inverted repeats (one is obtained by reversing and complementing the sequence of the other). For the contig indicated as u f,0 on the left-side green region, the contig in front of it (on the right-side green region) has the reverse orientation (so r). Because orientations are mutually exclusive, it is necessary to chose an occurrence that differs 0 (so 1). The couple u f,0 , ur,1 is defined as inverted fragments. Each inverted fragments couple is linked by blue dashed line. Inverted repeats can be modelled as sequences of nested inverted fragments. (b) Vertices correspond to genomic regions while edges indicate their neighbourhoods. Purple, green and red vertices correlate respectively with purple, green and red regions in Subfigure 1a. Two green vertices' label differ by their orientation 'f ' and 'r' (one is the reverse complement sequence of the other). They also differ by their occurrence number (0 and 1). It simply means that their associated genomic regions are simultaneously present in the genome sequence (i.e. the green vertices both participate in the solution). On the opposite, red vertices just differ by their orientation: it implies that only one genomic region from one of these two vertices appears in the genome (one genomic region is the reverse complement of the other).

( a )

 a Contigs data: c ∈ C c id c mult Links data: (orc, ord) ∈ L orc id orcor ord id ordor

Figure 3

 3 Figure 3 An instance of Multiplied Doubled Contigs Graph. The input contigs were multiplied by their multiplicity number, then doubled according to two DNA strands. The obtained graph possesses 42 vertices and 130 edges. Vertices candidate to participate in inverted repeats have one of their reverse oriented versions linked by a blue edge. They form a couple of vertices defined as inverted fragments. The solution (the assembled genome) is represented as a path in red. It passes through 8 adjacent inverted fragments that represent inverted repeats. It begins with the biggest vertex (a given starter, as illustrated in Subfigure 1a) and finishes in the same vertex since the chloroplast genomes are circular.

( a )Figure 4

 a4 Figure 4 Inverted fragments sets illustrations. (a) Two vertices i, j ∈ V 2 | i id = j id are inverted fragments if they are coloured identically. They belong to the InvF set. As an illustration of the invfrag function, note that invfrag(u f,0 ) = (u f,0 , ur,1), and invfrag(ur,3) = (u f,2 , ur,3). Inverted fragments linked by a dashed line represent pairs that belong to P InvF set (here three). Two inverted fragments (p, q) ∈ P InvF can be chosen in the solution only if they do not intersect. (b) (i, j) ∈ InvF and (k, l) ∈ InvF are two adjacent inverted fragments in MDCG. Remind that k ∈ N +

1 .

 1 can show that the number of constraints is O(|V | 2 + |E|) and the above formulation requires O(|V | 2 + |E|) variables. Nested Inverted Fragments Scaffolding Problem (NIFSP) is NP-Hard.

( a )Figure 5

 a5 Figure 5 Transformation of a directed graph G for LPSTP to a directed graph G for NIFSP. Edges in bold-green in both subfigures correspond to the solution path for LPSTP and NIFSP problems respectively. (a) As the longest path exits s and enters t, dashed edges do not participate in the solution. (b) Blue dashed line between vertices in G f 0 and G r1 visualise the reverse inverted fragments.

2 : 5 :

 25 if P ath[0] = s ∨ P ath[|P ath| -1] = t then 3: return F alse 4: for v ∈ P ath do if invfrag(v) ∈ InvF * then 6: u ← repeats.peek() add the associated vertex of v in their inverted fragment 11:

Figure 6

 6 Figure 6 Distribution of results' measures. Boxplots in violin plots show quartiles, where the whiskers correspond to 1.5 × IQR distances. (a) Histogram of misassemblies number occurrences. Here, only one instance outputs contigs with 2 misassemblies, and two instances output each one with one misassembly. (b) Distribution of solution genomes coverage (%) by output contigs. (c) Distribution of local errors number (per 100kb) on contigs output. While the first distribution shows indels errors type, the second shows mismatch errors type. (d) Distribution of Gurobi solving time (in seconds). Each dot represents the mean of the LP relaxation and B&B running times sum.

Table 1

 1 An instance of input data. It contains two types of data: (a)

Table 2

 2 Benchmark on artificial data. id: instance's identifier; |V | and |E| are respectively the number of vertices and edges in MDCG; time (in second): average LP relaxation time (above), and B&B time (below); opt: optimal value for the LP relaxation (above), the integer solution (below); %gap: gap value of the linear relaxation bound with respect to the optimal integer value; nodes: number of B&B nodes exploration; iter: number of iterations for the LP relaxation (above) and for the B&B phase (below). U 1, IR (so IR) and U 2 are the regions illustrated in Subfigure 1a, and respectively correspond to purple, green left/right-side, and red regions. |U 1|, |IR| and |U 2| report number of contigs (above) and sum of contig's multiplicity (below), for respectively the purple region, the green regions (one version) and the red region.

	id |V | |E| |U 1| |IR| |U 2|	time	opt % gap nodes	iter
	0	152 248	15 17	20 41	15 17	.07 .66	39 39	0	1	1690 3083
	1	304 512	30 33	40 85	30 33	1.43 3.55	79 79	0	1	5938 11109
	2	438 704	45 47	60 121	45 50	2.18 14.43	119 119	0	1	13992 36472
	3	592 976	60 65	80 168	60 62	11.22 25.3	160 159	.62	1	21672 42873
	4	742 1228	75 81	100 211	75 78	22.67 46.06	200.5 199	.75	1	32511 64963
	5	906 1524	90 99	120 256	90 97	89.13 202.24	242.5 239	1.44	1	53394 156103
	6 1030 1678	105 113	140 289	105 112	87.98 129.44	282.5 279	1.24	1	63680 122852
	7 1188 1974	120 130	160 338	120 125	313.59 373.97	321.5 319	.78	1	96122 161433
	8 1348 2242	135 147	180 377	135 149	235.91 3778.72	366.3798 359	2.01	1	152044 312150
	9 1504 2526	150 164	200 426	150 161	564.12 1012.3	404.3697 399	1.33	1	207623 493610

Table A . 3

 A3 Scaffolding results on chloroplast benchmark with reference genomes. Instance: species name; size: genome length in number of nucleotides (in base-pairs); |V | and |E| are respectively the number of vertices and edges; time: CPU solver time (in seconds); %gnm: percentage of genome length covered by the solution computed by our approach; NGA50: NG50 corrected of assembly errors; #mis: number of misassemblies in the proposed solution; |mis|: misassemblies total length (in bp); mism: number of alignment mismatches per 100 kbp; indels: number of alignment indels per 100 kbp.

	V. Epain, R. Andonov and D. Lavenier XX:20 Inverted Repeats Scaffolding										XX:19
	Table A.3, continued Table A.3, continued
	instance instance	size |V | |E| time %gnm NGA50 #mis |mis| mism indels size |V | |E| time %gnm NGA50 #mis |mis| mism indels
	Eriobotrya japonica Prunus dulcis	159156 34 58 157723 30 46 .02	100.0 87263.0 .02 99.81 85866.0	0	0	0	0	3.14	5.65 2.54	7.62
	Fragaria virginiana Scutellaria lateriflora	155577 48 72 152283 42 74 .04	99.87 85577.0 .07 99.99 84432.0	0	0	0	0	1.93	9.01 7.22	9.85
	Ginkgo biloba Theobroma cacao	156945 24 38 < .01 160619 36 52	.02	-99.97 89335.0 -	-	0	-	0	-	-7.47 11.83
	instance Gnetum gnemon Zea mays	size |V | |E| time %gnm NGA50 #mis |mis| mism indels 115022 20 32 .01 99.93 66153.0 0 0 5.22 5.22 140447 14 34 .01 99.98 82202.0 0 0 6.41 8.55
	Allium sativum Haberlea rhodopensis	153118 12 20 < .01 153099 76 120 .15 97.53 83933.0 99.99 82090.0	0	0	0	0 14.06 27.46 0.0	1.96
	Aloysia citriodora Hydrastis canadensis	154699 58 90 160000 14 26 .01	.07 100.0 87147.0 98.68 85159.0	0	0	0	0	12.45 27.51 0.0 2.5
	Althaea officinalis Illicium anisatum	159987 10 18 < .01 142723 18 26 < .01 100.0 100791.0 99.96 87918.0	0	0	0	0	3.5	0.63 7.01	0.0
	Amborella trichopoda Illicium floridanum	162686 16 24 < .01 143571 24 32 < .01 99.95 101328.0 99.92 90428.0	0	0	0	0 156.86 28.91 27.18 17.42
	Artemisia annua Illicium henryi	150952 38 54 143240 28 26 < .01	.01	99.99 82772.0 --	-	0	-	0	-	0.0	-	5.3
	Azolla filiculoides Jasminum tortuosum	147665 40 64 < .01 162080 74 122 .08 99.96 89117.0 -	-	-1 89992	-22.22 29.63 -	-
	Boswellia sacra Laurus nobilis	159228 38 46 < .01 152750 40 56 .02	98.51 40853.0 99.9 93540.0	0	1 84629 0	0.0	0.64 5.24	3.83
	Citrus limon Lupinus albus	160101 136 208 151921 44 128 .11	.2 99.99 82282.0 99.45 86976.0	0	0	0	0	32.03 30.15 9.22 11.85
	Coffea arabica Magnolia biondii	155188 122 194 160002 32 56 .03	.36 99.99 88127.0 99.07 84008.0	0	0	0	0	59.19 44.88 0.0 2.5
	Dendrobium nobile Magnolia denudata	152018 42 50 < .01 160089 14 26 .01 99.99 88149.0 96.55 33709.0	0	2 105240 1392.16 358.1 0 0.0 0.62
	Digitalis lanata Magnolia officinalis	153108 46 74 160136 52 76 .03	.05 99.99 88169.0 99.77 83262.0	0	0	0	0	3.75	6.55 4.37	5.24
	Dioscorea villosa Mitragyna speciosa	153974 64 100 155600 68 104	.1	.09 99.56 85060.0 99.95 83877.0	0	0	0	0 13.56 25.18 0.0	5.2
	Echinacea atrorubens Nicotiana tabacum	151912 24 22 < .01 155943 8 16 < .01 99.99 86706.0 -	-	0	-	0	-	0.0	-	0.0	-
	Echinacea purpurea Physcomitrella patens	151913 154 240 122890 36 40 < .01	.19	99.95 83480.0 --	-	0	-	0	-	7.24 15.81 -
	Echinacea sanguinea Pimenta dioica	151926 34 44 158984 126 178 .23	.02 98.01 84223.0 99.99 83674.0	0	0	0	0	11.19 6.42 19.25	7.24
	Echinacea speciosa Piper auritum	151860 26 42 159909 50 86	.1	.02 99.95 87529.0 99.99 83531.0	0	0	0	0 15.02 15.64 3.29	8.56
	Eleutherococcus senticosus Piper nigrum	156863 22 34 161523 50 82 .07	.01 99.86 89114.0 99.99 86842.0	0	0	0	0 20.46 19.84 8.29	5.74
					(the table continues on the next page) (the table continues on the next page)

The graph modeling described here has been yet applied in our previous scaffolding papers[START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF][START_REF] François | Global Optimization for Scaffolding and Completing Genome Assemblies[END_REF]. However, the notation MDCG is used here firstly.