the input data or to output the result. In the read assembly stage, especially in
the OLC approach, the fragments correspond to the reads and the links are the
overlaps between two reads. In the scaffolding stage, the fragments represent the
contigs and links can come from paired-end read or long read alignments against
the contigs.

The literature describes three graph structures: a digraph, a bigraph and an
ungraph. While the first one highlights visually the double-strand sequencing, its
weakness is that it requires two vertices for each read (Kececioglu, 1991). Myers
(1995) is the first to employ a bigraph to store overlaps. The key idea is to aggregate
the two orientations of a read into only one vertex, a link and its reverse into
only one edge. Finally, the ungraph structure associates one vertex for each of
the two extremities of a fragment to simplify the graph traversal (Huson et al.,
2002). The ungraph handles two sets of edges. A fragment-edge connects the
two vertices corresponding to the two extremities of a fragment (the tail and the
head). Two vertices representing the extremities of a fragment are connected by a
fragment-edge. The second one is the multiset of link-edges. A link-edge represents
both a link and its reverse.

Although some authors have provided conceptual descriptions of these graph
structures, the latter have never been compared before from an implementation
stat point of view. In this thesis, we have proposed an implementation design based
on adjacency lists for each graph structure. At the opposite of a compressed sparse
row or column, the adjacency lists are more suitable for adding or deleting vertices
or edges in the graphs. We have also described transformation processes to pass
from one implementation to another. We then have visualised the graph structures
and their implementations in a map.

We have retained three implementations: *DGS* and *DGF* for the multidigraph
and *BGU* for the multibigraph. After we have compared their memory consumption,
we theoretically measure the cost function for each of their algorithm on elementary
operations, such as iterating over the neighbours, adding or deleting a vertex or
an edge. We come to the conclusion that if memory is the critical issue, then
the *BGU* implementation should be preferred. Otherwise, we recommend *DGS*
and *DGF* for iterating over the neighbours and for deleting a vertex or an edge.
*BGU* is preferable for adding vertices or edges. To conclude, *DGF* proves to be
well-balanced and ideally tailored for dynamic graph operations. *DGF* is available
in a *Python3* package named *revsymg*<fnref n="1" /> and is easily installable via *PyPI*.

<footnote n="1">¹https://pypi.org/project/revsymg/</footnote>
# IV GLOBAL EXACT OPTIMISATIONS FOR CHLOROPLAST STRUCTURAL HAPLOTYPE SCAFFOLDING

Bebo & Chucho. (2008). Tres Palabras [Song]. On *Juntos Para Siempre*. Sony Music | Latin

<figure><img src="image_77.png" /></figure>

## In this chapter

<table><tr><td>1</td><td>Introduction</td><td>88</td></tr><tr><td>1.1</td><td>Chloroplast genome specificities.</td><td>88</td></tr><tr><td>1.2</td><td>State-of-the-art.</td><td>89</td></tr><tr><td>1.3</td><td>Our approach.</td><td>89</td></tr><tr><td>2</td><td>Input data and notation</td><td>90</td></tr><tr><td>2.1</td><td>Set of contigs C</td><td>90</td></tr><tr><td>2.2</td><td>Set of links L</td><td>91</td></tr><tr><td>2.3</td><td>Mathematically defining genomic regions</td><td>91</td></tr><tr><td>3</td><td>Chloroplast scaffolding problem formulations</td><td>92</td></tr><tr><td>4</td><td>Graph and repeated fragment sets</td><td>96</td></tr><tr><td>4.1</td><td>Graph structure</td><td>97</td></tr><tr><td>4.2</td><td>Repeated fragment sets.</td><td>98</td></tr><tr><td>5</td><td>Integer Linear Programming (ILP) formulation</td><td>102</td></tr><tr><td>5.1</td><td>Circuit constraints</td><td>103</td></tr><tr><td>5.2</td><td>Repeated regions constraints</td><td>104</td></tr><tr><td>5.3</td><td>Fixing regions constraints.</td><td>107</td></tr><tr><td>5.4</td><td>Speed-up constraints</td><td>108</td></tr><tr><td>5.5</td><td>Scaffolding problems ILP</td><td>108</td></tr></table>
<table><tr><td>6</td><td>Hierarchical problem succession</td><td>109</td></tr><tr><td>7</td><td>From an ILP solution to a genome structure</td><td>110</td></tr><tr><td>8</td><td>Multiple genome forms</td><td>114</td></tr><tr><td>9</td><td>NP-completeness</td><td>117</td></tr><tr><td>10</td><td>Numerical results</td><td>121</td></tr><tr><td>10.1</td><td>Complexity validation on artificial data</td><td>122</td></tr><tr><td>10.1.1</td><td>Perfect artificial data.</td><td>122</td></tr><tr><td>10.1.2</td><td>Noisy artificial data</td><td>122</td></tr><tr><td>10.2</td><td>Synthetic chloroplast input data.</td><td>123</td></tr><tr><td>10.2.1</td><td>Input data generation</td><td>123</td></tr><tr><td>10.2.2</td><td>The evaluation's metrics</td><td>125</td></tr><tr><td>10.2.3</td><td>Initial version</td><td>126</td></tr><tr><td>10.2.4</td><td>Modified version.</td><td>129</td></tr><tr><td>11</td><td>Conclusion</td><td>130</td></tr><tr><td>12</td><td>Discussion and perspectives</td><td>131</td></tr></table>

In this chapter, we focus on the scaffolding problem specific to chloroplast genomes. We first recall the specificities of these genomes as described in the introduction (Chapter I, Section 2), highlight the state-of-the-art in chloroplast genome assembly (Chapter II, Section 3) and present our approach.

# 1 Introduction

## 1.1 Chloroplast genome specificities

Chloroplasts are plants' organelles derived from the integration of a cyanobacterium in an eukaryotic host. They conduct photosynthesis, a process to convert light energy into chemical energy. Over the evolution time, the chloroplast genome has reduced in length and loosened in terms of complexity (Xiao-Ming et al., 2017). As a result, chloroplast genomes possess few repeats that are usually identical in nucleotide sequences. One the most studied forms of chloroplast genome is a circular quadripartite DNA molecule. It consists of four regions: two identical (or highly similar) nucleotide subsequences, separated by two single-copies (Long Single-Copy, LSC, and Short Single-Copy, SSC) (Palmer, 1985; Bock and Knoop, 2012). There are two types of repeats: (i) the Direct Repeat (DR), where the sequences are highly similar; (ii) the Inverted Repeat (IR), where one sequence is
the reverse-complement of the other. Figure 5 illustrates the common chloroplast genome architectures.

Furthermore, each chloroplast has multiple copies of its genome, and the molecular forms of the copies differ (*structural haplotypes* leading to heteroplasmy, and multigenomic structures – not discussed here, Palmer 1983; Bendich 2004). This phenomenon can be induced by flip-flop inversion: one subsequence is reverse-complemented (*reversed*) during the DNA replication. This inversion is provoked by the existence of facing inverted repeats on either side of the reversed subsequence.

## 1.2 State-of-the-art

Although there are chloroplast genome assemblers and scaffolders, they do not fully exploit the chloroplast genome's specificities. Some of them are pipelines of generic methods applied on cleaned input dataset (Ankenbrand et al., 2018), or based on locally approaches as seed-and-extend algorithms (Coissac et al., 2016; Dierckxsens et al., 2017). Jin et al. (2020), in *GetOrganelle*, statistically compute the contigs' multiplicities by minimising the squared distance between them and the mapping coverage by the reads.

Concerning the handling of the distinct genome forms, *GetOrganelle* returns several solutions and explores in post-process the corresponding architectures. In Andonov et al. (2019) the flip-flop inversion breakpoints are detected in a post-scaffolding-process, and new optimal solutions can be constructed in polynomial time.

## 1.3 Our approach

We raise the following two questions: (i) how to mathematically model chloroplast genomic biological knowledge? (ii) How to reveal the structural haplotypes through the scaffolding problem formulation?

To formalise the scaffolding by integrating the structural haplotypes in its core, we postulate on the particularities of the chloroplast genome: (i) repeats are pairs of regions; (ii) the two regions of a repeat have identical (or reversed) nucleotide sequence; (iii) structural haplotypes can be seen as permutations in a sequence of oriented contigs.

We first describe the input data for our approach and provide mathematical definitions (Section 2). Based on the above three assumptions, we propose a new formulation for scaffolding chloroplast genomes without requiring any distances (Section 3). Our approach is region-driven and focuses on retrieving the repeats first. We model the optimisation problem on a digraph where we apply several ILP (Sections 4 and 5). We then detail how we combine the ILP solutions (Section 6).
The ILP's solutions and the digraph correspond to an oriented contig sequence representing one genome's form (Section 7). We partition this sequence into genomic regions we express in a region graph. This graph enables us to return multiple genome forms (Section 8).

We prove the decision version of the chloroplast genome scaffolding problem to be *NP-complete* (Section 9). However, we exactly solve the problem by profiting from the small size of the chloroplast instances and providing some numerical results (Section 10). We finally conclude (Sections 11 and 12).

## 2 Input data and notation

<caption>**Table 11 – Toy example of input data.**</caption>

(a) Set of contigs $C$. (b) Set of links $\mathcal{L}$. Each row represents a link. For the sake of space, for no one of the links in the table, its reverse is given, although it belongs to $\mathcal{L}$.

<table><caption>(a) Contig set $C$</caption><thead><tr><th>contig</th><th>mult</th><th>wex</th></tr></thead><tbody><tr><td>a</td><td>1</td><td>0.70</td></tr><tr><td>b</td><td>2</td><td>0.83</td></tr><tr><td>c</td><td>2</td><td>0.17</td></tr><tr><td>d</td><td>1</td><td>0.43</td></tr></tbody></table>

<table><caption>(b) Link set $\mathcal{L}$</caption><thead><tr><th>contig</th><th>orient</th><th>contig</th><th>orient</th></tr></thead><tbody><tr><td>a</td><td>f</td><td>c</td><td>r</td></tr><tr><td>a</td><td>r</td><td>c</td><td>r</td></tr><tr><td>b</td><td>r</td><td>c</td><td>f</td></tr><tr><td>b</td><td>f</td><td>d</td><td>f</td></tr><tr><td>b</td><td>f</td><td>d</td><td>r</td></tr></tbody></table>

### 2.1 Set of contigs $C$

Contigs are words in the nucleotides DNA alphabet $\Sigma_{nuc}^{+}$. A contig can occur in the genome up to an integer called *multiplicity*. Function *mult*: $C \to \mathbb{N}_{>0}$ provides its value. Any of contig's *occurrence* can appear in one of two possible reverse-complementary and mutually exclusive *orientations*: forward ($f = 0$) and reverse ($r = 1$). Each contig is provided with an *existence-weight* in $\mathbb{R}_{\ge 0}$ given by the function *wex*. The weight is proportional to the number of times a contig aligns with chloroplast sequences from a given set (from related or unrelated species). Finally, one contig in this set is defined as the *starter* ($s$) that must uniquely participate in the genome (*mult*($s$) = 1). The starter is a contig whose sequence matches a sequence shared by most chloroplast genomes in a single-copy. Table 11a gives an example of contig set.
## 2.2 Set of links $\mathcal{L}$

Each link is an ordered pair of oriented contigs. We denote by $\mathcal{L} \subset (\mathcal{C} \times \{f, r\})^2$ the link set<fnref n="1" />. The nature of the double-strands DNA requires that $\forall (c, d) \in \mathcal{L}, (\bar{d}, \bar{c}) \in \mathcal{L}$, where $\bar{c}$ and $\bar{d}$ denote the oriented contigs c and d in their reverse orientation, respectively (note that $\bar{c} = c$). The links between two oriented contigs c and d are valid for all occurrences of c and d respecting the same orientations. Table 11b gives an example of link set.

## 2.3 Mathematically defining genomic regions

We aim to order oriented occurrences of the contigs based on their links. Each genome form corresponds to a sequence of oriented contigs. Not all contigs or their occurrences are included. Indeed, the contig set may contain contigs belonging to the plant genome or other organelles. The link set may also contain artefact links. Definition 2.1 provides the properties the sequence of oriented contigs must respect.

### ► Definition 2.1: Sequence of oriented contigs

Let $SOC = (c_0, c_1, \dots, c_{n-1})$ be a sequence of oriented contigs:

$$ - \forall i \in [0, n-1], (c_i, c_{i+1}) \in \mathcal{L}; $$

$$ - \forall c \in \mathcal{C}, \sum_{c_i \in SOC|c_i=c} 1 \leq mult(c). $$

Based on the biological knowledge, we address the dedicated chloroplast scaffolding problem as a region-driven scaffolding, such that specific regions fit into a circular structure. We identify three types of regions to scaffold: the Inverted Repeat (IR), the Direct Repeat (DR) and the Single-Copy (SC).

### ► Definition 2.2: Region

A region $r = (c_0, c_1, \dots, c_{n-1})$ is a sequence of oriented contigs. Each region is oriented. Let $\bar{r} = (\overline{c_{n-1}}, \dots, \overline{c_1}, \overline{c_0})$ be the reverse region of $r$. It is composed of the oriented contigs of $r$, considered in their reverse orientation, and given in the reversed order. According to the reverse symmetry in the links, $\bar{r}$ also respects Definition 2.1.

<footnote n="1">¹In Chapter II, Section 1, $\mathcal{L}$ is a multiset while in this chapter, $\mathcal{L}$ is a set.</footnote>
► **Definition 2.3: Direct repeat (DR)**

A DR is a couple of regions ($dr_i, dr_j$) where $dr_i = dr_j$.

► **Definition 2.4: Inverted repeat (IR)**

An IR is a couple of regions ($ir_k, ir_l$) where $ir_l = \overline{ir_k}$.

► **Definition 2.5: Repeat**

A repeat is the generic term to denote either DR or IR. The length *replen*(*R*) of a repeat *R* = (*r*<sub>*i*</sub>, *r*<sub>*j*</sub>) equals the lengths of *r*<sub>*i*</sub> and *r*<sub>*j*</sub> (*replen*(*R*) = |*r*<sub>*i*</sub>| + |*r*<sub>*j*</sub>|).

► **Definition 2.6: Single-copy (SC)**

A SC is a region that is not part of a repeat.

► **Definition 2.7: Region weight**

The weight $rwex(r)$ of a region $r$ is defined as $rwex(r) = \sum_{c \in r} wex(c)$.

A chloroplast genome consists of a sequence of oriented regions. A genome form is a result of iterative transformations of an initial one. Section 8 introduces the region graph to model multiple genome forms (sequences of oriented contigs).

► **Definition 2.8: Sequence of oriented regions**

Consider a sequence SOR = (*r*<sub>0</sub>, *r*<sub>1</sub>, . . . , *r*<sub>*n*−1</sub>):

$$
- \forall i \in [0, n-1], (r_i[|r_i| - 1]^a, r_{i+1}[0]) \in \mathcal{L};
$$

$$
- \forall c \in \mathcal{C}, \sum_{r \in \text{SOR}} \sum_{c_i \in r | c_i = c} 1 \leq \text{mult}(c).
$$

3 Chloroplast scaffolding problem formulations

Solving the repeats is the most challenging task. A formulation that does not
restrict the occurrences can lead to misassemblies where the results are longer than

<footnote>aThe notation $x[i]$ denotes the $i^{th}$ element in $x$.</footnote>
the solution genomes. Therefore, the use of an occurrence is limited to conformity with the biological knowledge of genome forms. A contig should participate in the sequence only if it enables the formation of pairs of repeated regions. In this case, we would be inclined to assemble the minimum number of repeats.

However, in the case of repeat degeneration (e.g. two subsequences inside the two regions of an identified repeat differ, note that some IR losses have been reported in the chloroplast genomes of green algae - Turmel et al. 2017) finding the minimum number of repeats is not an appropriate model. Figure 31 illustrates the impact of degenerations on quadripartite structures. Indeed, in Figure 31b, we cannot guaranty to find both $IR_1$ and $IR_2$, but perhaps only one of them. For each repeat type, we address this issue by maximising the cumulative repeat lengths only if their regions respect a specific order.

<figure><img src="image_78.png" /><figcaption>■ **Figure 31 – Repeat degeneration and region orders.**</figcaption></figure>

With time, a repeat may degrade so that its occurrences differ. (a) and (c) show two common quadripartite structures with an IR and a DR, respectively. (b) and (d) highlight the impact of a degeneration on their structures. We give the below region orders according to the LSC arrow's direction. In (b), the degeneration results in two IRs: $IR_1 = (i, j)$ and $IR_2 = (k, l)$, such that $i$ is before $k$ and $l$ precedes $j$. In (d) it results in two DRs: $DR_1 = (i, j)$ and $IR_2 = (k, l)$, such that $i$ is before $k$ and $j$ precedes $l$.

► **Definition 3.1: Chloroplast scaffolding problem CHSP**

Given a set of contigs with their multiplicities and their weights, a starting contig and a link set. The aim is to obtain a circular sequence of oriented regions maximising the cumulative repeat lengths and minimising the number of repeats, with single-copies of maximum-weight.

For instance, let Cases (A) and (B) be two distinct and feasible sequences of oriented contigs:

(A) $(..., a, b, c, d, ..., a, b, c, d, ...)$ has one DR $(i, j)$, where $i = j = (a, b, c, d)$;
(B) $(... , a, b, ..., c, d, ..., a, b, ..., c, d, ...)$ has two DRs $(k, l)$ and $(m, n)$, where $k = l = (a, b)$ and $m = n = (c, d)$.

For Cases (A) and (B) the cumulative lengths are the same ($replen((i, j)) = replen((k, l)) + replen((m, n)) = 8$). However, Case (A) has one less repeat, which we prefer.

CHSP involves three subproblems, each one associated with a particular type of region: (i) DRP for the direct repeats (Definition 3.2); (ii) IRP for the inverted repeats (Definition 3.3) and (iii) SCP for the single-copies (Definition 3.4). We tackle CHSP in a hierarchical succession of DRP, IRP and SCP (Definition 3.6). Any intermediate problem must preserve the regions found by its predecessors (Definition 3.5).

DRP and IRP constrain the number of occurrences to the structure of pairs of repetitions. Indeed, each repeat type defines a valid repeat structure. The problems consist in maximising the cumulative length of the minimum number of repeats.

### ► Definition 3.2: Chloroplast direct repeat scaffolding problem DRP

Consider a set of contigs, their multiplicities, a starting contig and a link set. Find a circular sequence of oriented regions SOR, such that:

- it maximises the cumulative length of the minimum number of DRs, joined by regions of any kind;

- for any couple of DRs $(i, j)$ and $(k, l)$ found in SOR such that their respective positions in SOR given by function $\sigma$ respect $\sigma(i) < \sigma(j)$, $\sigma(k) < \sigma(l)$, and $\sigma(i) < \sigma(k)$, then:

$$ - [[\sigma(i), \sigma(j)]] \cap [[\sigma(k), \sigma(l)]] = \emptyset; $$

$$ - \text{ or } [[\sigma(k), \sigma(j)]] \subset [[\sigma(i), \sigma(l)]] $$

### ► Definition 3.3: Chloroplast inverted repeat scaffolding problem IRP

Consider a set of contigs, their multiplicities, a starting contig and a link set. Find a circular sequence of oriented regions SOR, such that:

- it maximises the cumulative length of the minimum number of IRs, joined by regions of any kind;

- for any couple of IRs $(i, j)$ and $(k, l)$ found in SOR such that their respective positions in SOR given by function $\sigma$ respect $\sigma(i) < \sigma(j)$, $\sigma(k) < \sigma(l)$ and $\sigma(i) < \sigma(k)$, then:
$$- [[\sigma(i), \sigma(j)]] \cap [[\sigma(k), \sigma(l)]] = \emptyset;$$

$$- \text{ or } [[\sigma(k), \sigma(l)]] \subset [[\sigma(i), \sigma(j)]].$$

Figure 32 provides examples of valid oriented contig positioning for each common genome structure (Figure 5). Although Figure 36 illustrates the authorised and forbidden positions for the latter defined repeated fragments, it is also applicable for the DRP and IRP regions’ position cases.

### ► Definition 3.4: Chloroplast single-copy scaffolding problem SCP

Consider a set of contigs, their multiplicities, their weights, a starting contig and a link set. Find a circular sequence of oriented regions such that all the single-copies maximise their weights.

Note that in Definition 3.4, if they are no repeats, the problem is reduced to find the maximum-weighted circuit of oriented contigs.

### ► Definition 3.5: Chloroplast scaffolding problem succession

DRP, IRP and SCP (Definitions 3.2 to 3.4) can also take as input a set of fixed regions that must be preserved in the resulting sequence of oriented regions.

Our hierarchical approach prioritises the scaffolding of the repeats previously to the SC regions. Indeed, scaffolding the repeats is the most difficult task as it can lead to misassemblies (wrongly chosen links). Hence, a contig should only be used as many times as possible if its occurrences enable the scaffolding of repeats that most closely represent the architecture of the chloroplast genome, as illustrated in Figures 5 and 32. We assume that the weights concern events that are less relevant comparing to the genomes architecture. Also, the weights on the contigs are less relevant than if they were weights on the links.

### ► Definition 3.6: Hierarchical problem succession

The form of each solution of CHSP satisfies one of the two problem successions: DRP-IRP-SCP ($h_1$) and IRP-DRP-SCP ($h_2$).

The next question is how to prioritise DRP and IRP? We propose resolving the order by comparing the scores defined in Section 5.5: if DRP score is better than this of IRP, then the retained succession will be DRP-IRP-SCP, otherwise it will be IRP-DRP-SCP. In the equality case, we discriminate at a further step of the hierarchical successions. The process is detailed in Section 6.
Finally, each hierarchical problem succession produces a circular sequence of oriented regions. From the obtained sequence it is possible to extract a set of ordered pairs of oriented regions. This procedure allows the building of several circular sequences of oriented regions of the same length. Each of them represents one possible chloroplast genome form. This all-equivalent-form process is described in Section 8.

<figure><img src="image_79.png" /><figcaption>■ **Figure 32 – Chloroplast repeat scaffolding.**</figcaption></figure>

Each subfigure is a common chloroplast genome structure with its associated order of oriented contigs (coloured arrows). The green and the blue sequences of arrows are IR and DR, respectively. The purple and the red ones are single-copy regions. Contig *s* is the starter, and the right side black arrow determines the contigs' order. Contigs *a*₀, *a*₁, *b*₀, *b*₁, *c*₀, *c*₁ and *d*₀, *d*₁ are two occurrences of four contigs *a*, *b*, *c* and *d*, respectively. Each coloured dashed line links two occurrences of the same contig. (a) The order of the occurrences is reversed, and their arrows are oppositely oriented. Visually, an IR produces parallel dashed lines. (b) The order and the orientation of the occurrences is preserved, revealing a DR. (c) A chloroplast genome can contain the two types of repeats. Here, we will retain the hierarchical problem succession *IRP-DRP-SCP* (*h*₂) since the IR contains more contigs than the DR.

4 Graph and repeated fragment sets

In order to efficiently handle the multiplicities of the contigs, hence of the links, we need to build adapted data structures. On the one hand, finding a sequence of oriented contigs, when the links correspond to ordered pairs of oriented contigs, justifies the use of a directed graph to represent the oriented contigs and their links. Section 4.1 defines such a directed graph structure. On the other hand, scaffolding the repeats requires choosing pairs of contigs occurring several times in the oriented contig sequence. Section 4.2 defines the sets of such repeat contig candidates.
## 4.1 Graph structure

Here we describe a directed graph suitable for further algorithms and the mathematical formulation of the scaffolding problems from Definitions 3.2 to 3.4.

► **Definition 4.1: Multiplied Doubled Contig Graph – MDCG**

Given a set of contigs $\mathcal{C}$, their multiplicities and the link set $\mathcal{L}$, the multiplied doubled contig graph $\mathcal{MDCG} = (V, E, vwex)$ is defined such thatª:

- $V = \{v_{f,0}, \dots, v_{f,n-1}, v_{r,0}, \dots, v_{r,n-1} \mid c \in \mathcal{C}, n = \text{mult}(c)\}$ is the set of all the forward and reverse occurrences of all the contigs ($|V| = 2 \sum_{c \in \mathcal{C}} \text{mult}(c)$). The vertices are associated with four functions:

- *contig*: $V \to \mathcal{C}$ provides the contig associated with a vertex;

- *vor*: $V \to \{f, r\}$ provides the orientation of the contig;

- *vocc*: $V \to \mathbb{N}_{>0}$ provides the occurrence number of the contig;

- *vwex*: $V \to \mathbb{R}_{\ge 0}$ provides the weight of each vertex such that $\forall v \in V, vwex(v) = wex(contig(v))$.

- $E = \{(u, v) \in V^2 \mid ((contig(u), vor(u)), (contig(v), vor(v))) \in \mathcal{L}\}$ is the set of multiplied links ($|E| = \sum_{(c,d) \in \mathcal{L}} \text{mult}(c)\text{mult}(d)$).

Figure 33 illustrates the $\mathcal{MDCG}$ representing the example data given in Table 11.

<figure><img src="image_80.png" /></figure>

■ **Figure 33 – MDCG example.**

Each vertex is associated with an occurrence of an oriented contig, and each contig is represented by an even number of vertices. For example, vertex labelled $c_{r,1}$ means that it comes from contig $c = contig(c_{r,1})$, in its reverse orientation ($vor(c_{r,1}) = r$), and in its second occurrence ($vocc(c_{r,1}) = 1$). The colours are the same as the ones in Figure 5a to relate the input data with the IR architecture. The bold red edges draw a circuit corresponding to an IR scaffolding where $a_{f,0}$ is the starter.

<footnote>ªMDCG is a digraph, not a multidigraph as DG.</footnote>
## 4.2 Repeated fragment sets

A repeat is a couple containing two identical (or reverse for IR) sequences of oriented contigs (Definitions 2.3 and 2.4). Therefore, a repeat consists of couples containing two identical (or reverse) contigs. In the context of MDCG, this leads to the concept of repeated fragments.

### ► Definition 4.2: Repeated fragment

A repeated fragment is an unordered pair of vertices such that one of the corresponding oriented contig belongs to the first region of a repeat, while the other belongs to the second region. The vertices are associated with the same contig but their occurrences differ, i.e. for each repeat $(r_i, r_j)$, $\exists u, v \in V, c \in C$ where $contig(u) = contig(v) = c$ and $vocc(u) \neq vocc(v)$ such that $(c, vor(u)) \in r_i$ and $(c, vor(v)) \in r_j$.

For example, $(c_{f,0}, c_{r,1})$ is a repeated fragment for the IR in Figure 33. We then precise the set of repeated fragments for each repeat type. Denote by $\mathcal{R} = \{c \in C | mult(c) > 1\}$ the set of contigs candidate to be part of repeats. For the sake of clarity, for each vertex $v \in V$, we note $ctg_v = contig(v)$, $or_v = vor(v)$, $occ_v = vocc(v)$, $wex_v = vwex(v)$ and $mult_v = mult(contig(v))$. We also assume there is an arbitrary strict total order on $C$, i.e. $\forall c, d \in C, c \neq d \iff c < d \lor c > d$.

### ► Definition 4.3: Set of direct fragments DirF

A direct fragment $(u, v) \in V^2$ is a repeated fragment, such that $u$ and $v$ have the same orientation.

$$ DirF = \bigcup_{c \in \mathcal{R}} \left\{ \begin{array}{l} (i, j) \in V^2 \text{ s.t.} \\ ctg_i = ctg_j = c \\ \land or_i = or_j \in \{f, r\} \\ \land occ_i = occ_j - 1 = 2k \\ 0 \le k < \left\lfloor \frac{mult(c)}{2} \right\rfloor \end{array} \right\} $$

### ► Definition 4.4: Set of inverted fragments InvF

An inverted fragment $(u, v) \in V^2$ is a repeated fragment, such that the orienta-
<figure><img src="image_81.png" /><figcaption>(a) Direct fragments and their pairs</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>(b) Inverted fragments and their pairs</figcaption></figure>

**Figure 34** – Repeated fragment sets illustration for two contigs *c* and *d*. In the two subfigures, *mult*(*c*) = 4 and *mult*(*d*) = 2, so that *contig*(*u*<sub>or,occ</sub>) = *c* and *contig*(*v*<sub>or,occ</sub>) = *d*. Two vertices coming from the same contig are respectively direct/inverted fragments if they are in the same coloured box, and so they belong to *DirF/InvF*. A tight grey line connects two direct/inverted fragments if their pair belong to *PDirF/PInvF*. **(a)** |*DirF*| = 6, and, e.g. (*u*<sub>f,0</sub>, *u*<sub>f,1</sub>) ∈ *DirF* so *dirfrag*(*u*<sub>f,0</sub>) = *dirfrag*(*u*<sub>f,1</sub>) = (*u*<sub>f,0</sub>, *u*<sub>f,1</sub>). |*PDirF*| = 12, and, e.g. ((*u*<sub>f,0</sub>, *u*<sub>f,1</sub>), (*u*<sub>r,2</sub>, *u*<sub>r,3</sub>)) ∈ *PDirF*. **(b)** |*InvF*| = 3, and, e.g. (*u*<sub>f,2</sub>, *u*<sub>r,3</sub>) ∈ *InvF* so *invfrag*(*u*<sub>f,2</sub>) = *invfrag*(*u*<sub>r,3</sub>) = (*u*<sub>f,2</sub>, *u*<sub>r,3</sub>). |*PInvF*| = 3, and, e.g. ((*u*<sub>f,2</sub>, *u*<sub>r,3</sub>), (*v*<sub>f,0</sub>, *v*<sub>r,1</sub>)) ∈ *PInvF*.

tions of *u* and *v* differ.

$$
\text{InvF} = \bigcup_{c \in \mathcal{R}} \left\{ \begin{array}{l} (i, j) \in V^2 \text{ s.t.} \\ \quad ctg_i = ctg_j = c \\ \quad \wedge or_i = f \wedge or_j = r \\ \quad \wedge occ_i = occ_j - 1 = 2k \\ \quad 0 \le k < \left\lfloor \frac{\text{mult}(c)}{2} \right\rfloor \end{array} \right\}
$$

Figure 34 illustrates *DirF* and *InvF* sets. In addition, we add two functions to retrieve the repeated fragments from the vertices in $\Theta(1)$: *dirfrag*: $V' \subset V \rightarrow DirF$ and *invfrag*: $V' \subset V \rightarrow InvF$ (abstracted with the *repfrag* function, c.f. Appendix 1 for their definitions).

Furthermore, Definitions 3.2 and 3.3 constrain the order between the repeated
fragments. Hence, they respectively require comparing pairs of direct/inverted
fragments, that must be defined:

► **Definition 4.5: Set of pairs of direct fragments**

$$ PDirF = \left\{ \begin{array}{l} \left( (i, j), (k, l) \right) \in DirF^2 \text{ s.t.} \\ \quad ctg_j < ctg_k \\ \quad \vee \\ \quad ctg_j = ctg_k \wedge occ_j < occ_k \end{array} \right\} $$

► **Definition 4.6: Set of pairs of inverted fragments**

$$ PInvF = \left\{ \begin{array}{l} \left( (i, j), (k, l) \right) \in InvF^2 \text{ s.t.} \\ \quad ctg_j < ctg_k \\ \quad \vee \\ \quad ctg_j = ctg_k \wedge occ_j < occ_k \end{array} \right\} $$

Figure 34 illustrates *PDirF* and *PInvF* sets. The constraints defining *DirF*, *InvF*, *ADirF*, *AInvF*, *PDirF* and *PInvF* are explained in details in Appendix 2 where we proof they are the smallest sets enabling to find all the distinct solutions.

Furthermore, a repeat is a couple of regions (Definition 2.5), themselves defined as oriented contig sequences (Definition 2.2). We need to define the edges connecting two repeated fragments.

► **Definition 4.7: Set of adjacent repeated fragments**

$$ \text{An adjacent repeated fragment is an edge } (u, v) \in E \text{ such that } u \text{ and } v \text{ participate in two distinct repeated fragments.} $$

► **Definition 4.8: Set of adjacent direct fragments**

$$ \text{An adjacent direct fragment is an edge between two direct fragments. Let} $$
<figure><img src="image_83.png" /><figcaption>■ **Figure 35 – Adjacent repeated fragment sets examples.**</figcaption></figure>

The two subfigures represent the multiplied link (and its reverse) $((c, f), (d, f)) \in \mathcal{L}$, where $mult(c) = 2$ and $mult(d) = 2$, so that $contig(u_{or,occ}) = c$ and $contig(v_{or,occ}) = d$. Two vertices of the same colour visualise a repeated fragment. Bold edges (canonical) are the ones that belong to the adjacent repeated fragments sets. The functions $diradj/invadj$ enable to retrieve the normal edges with the bold ones, and vice-versa. Dashed edges do not participate in $ADirF/AInvF$. Remember that $\forall (u, v) \in E, (\bar{v}, \bar{u}) \in E$. **(a) $diradj(u_{f,0}, v_{f,0}) = (u_{f,1}, v_{f,1})$; (b) $invadj(u_{f,0}, v_{f,0}) = (v_{r,1}, u_{r,1})$**.

$ADirF$ be the set of adjacent direct fragments:

$$ ADirF = \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ ctg_u \neq ctg_v \\ \wedge occ_u = 2k \\ 0 \le k < \left\lfloor \frac{mult_u}{2} \right\rfloor \\ \wedge occ_v = 2k' \\ 0 \le k' < \left\lfloor \frac{mult_v}{2} \right\rfloor \end{array} \right\} \cup \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ ctg_u = ctg_v \\ \wedge (or_u = f \lor or_v = f) \\ \wedge occ_u = 2k \\ \wedge occ_v = 2k' \\ 0 \le k < k' < \left\lfloor \frac{mult_u}{2} \right\rfloor \end{array} \right\} $$

► **Definition 4.9: Set of adjacent inverted fragments**

An adjacent inverted fragment is an edge between two inverted fragments. Let
$AInvF$ be the set of adjacent inverted fragments:

$$AInvF = \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ \text{ctg}_u < \text{ctg}_v \\ \wedge \text{occ}_u = 2k + \text{or}_u \\ 0 \le k < \left\lfloor \frac{\text{mult}_u}{2} \right\rfloor \\ \wedge \text{occ}_v = 2k' + \text{or}_v \\ 0 \le k' < \left\lfloor \frac{\text{mult}_v}{2} \right\rfloor \end{array} \right\} \cup \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ \text{ctg}_u = \text{ctg}_v \\ \wedge (\text{or}_u = f \lor \text{or}_v = f) \\ \wedge \text{occ}_u - \text{or}_u = 2k \\ \wedge \text{occ}_v - \text{or}_v = 2k' \\ 0 \le k < k' < \left\lfloor \frac{\text{mult}_u}{2} \right\rfloor \end{array} \right\}$$

Edges in $ADirF$ and $AInvF$ play the role of *canonical edges* between two adjacent repeated fragments, see Figure 35. In addition, we add two functions to retrieve the adjacent repeated fragments from the edges in $\Theta(1)$: *diradj*: $E \to E$ and *invadj*: $E \to E$ (abstracted with the *repadj* function, c.f. Appendix 1 for their definitions).

# 5 Integer Linear Programming (ILP) formulation

Modelling *DRP*, *IRP* and *SCP* from Definitions 3.2 to 3.4 requires finding a valid circuit in *MDCG*.

## ► Definition 5.1: Valid circuit in MDCG

Given a graph $MDCG = (V, E)$ and a starting vertex $s$, where $ctg_s$ is the starting contig, $or_s = f$ and $occ_s = 0$. A circuit $cp$ in $MDCG$ is valid if:

- it starts and ends with $s$;

- $\forall v \in cp, \bar{v} \notin cp$, where $ctg_v = ctg_{\bar{v}}, or_v = 1 - or_{\bar{v}}$ and $occ_v = occ_{\bar{v}}$;

- consecutive vertices $u$ and $v$ in $cp$ are connected by an edge $(u, v) \in E$.

First we describe common constraint blocks in Sections 5.1 to 5.3 for ILPs formulations, and then we give the *DRP*, *IRP* and *SCP* scaffolding problems ILP in Section 5.5.

Let $M = \sum_{c \in C} \text{mult}(c)$ be a constant, $N_v^-$ and $N_v^+$ be the sets of predecessors and successors of vertex $v \in V$, respectively.
## 5.1 Circuit constraints

The following set of constraints defines a valid circuit of oriented contig in MDCG, and is defined with a flow formulation as in Andonov et al. (2019) instead of using Miller-Tucker-Zemlin constraints to avoid cycles (Miller et al., 1960).

### Binary variables

— $x_e$ encodes if the edge $e \in E$ participates in the circuit.

### Continuous variables

— $i_v \in [0, 1]$ encodes if the vertex $v \in V \setminus \{s, \bar{s}\}$ participates in the circuit. Although it is a continuous variable, it acts as a binary one as proven in (François et al., 2018).

— $f_e \in \mathbb{R}_{\ge 0}$ is the positive flow on the participating edge $e \in E$ in the circuit (zero otherwise).

Constraint C1 defines the flow. The circuit must start and end with the starter in its forward orientation (Constraints C2 to C5). If a vertex participates, its reverse cannot (Constraint C6). Defining a circuit is equivalent to requiring an edge to exit a vertex if it has an incoming one (Constraint C7). Constraint C8 forces the flow to be monotonically increasing. This property avoids cycles.

### CCircuit constraints

$$x_e \le f_e \le Mx_e \qquad \forall e \in E \qquad (C1)$$

$$\sum_{v \in N_s^+} x_{sv} = \sum_{v \in N_s^-} x_{vs} = 1 \qquad (C2)$$

$$\sum_{v \in N_s^+} f_{sv} = 1 \qquad (C3)$$

$$x_{\bar{v}s} = 0 \qquad \forall v \in N_s^- \qquad (C4)$$

$$x_{sv} = 0 \qquad \forall v \in N_s^+ \qquad (C5)$$
$$ \forall v \in V \setminus \{s, \bar{s}\} : \\ i_v + i_{\bar{v}} \le 1 \qquad (C6) $$

$$ \sum_{u \in N_v^-} x_{uv} \le i_v \le \sum_{w \in N_v^+} x_{vw} \qquad (C7) $$

$$ \sum_{w \in N_v^+} f_{vw} - \sum_{u \in N_v^-} f_{uv} = i_v \qquad (C8) $$

$$ x_e \in \{0, 1\} \quad \forall e \in E $$

$$ i_v \in [0, 1] \quad \forall v \in V \setminus \{s, \bar{s}\} $$

$$ f_e \in \mathbb{R}_{\ge 0} \quad \forall e \in E $$

## 5.2 Repeated regions constraints

The following constraints are general to define ILPs for DRP and IRP. Definitions 3.2 and 3.3 define the repeated regions according to the positions of the oriented contig in them. It follows that some order of the vertices in the pairs of repeated fragments are allowed, and some others are forbidden. We decide to write the constraints for the forbidden cases because they are fewer than the allowed ones. To model the forbidden orders between 4 vertices, we compare the positions between two.

Specifically for IRP, modelling the forbidden orders echoes the approach for the RNA folding problem (Gusfield, 2019), except that the positions of the RNA's nucleotides are known.

According to Definitions 3.2 and 3.3, and given PDirF and PInvF (Definitions 4.5 and 4.6), denote by ForbidDR and ForbidIR the sets of forbidden quartet vertices for the DRs and IRs, respectively:

$$ ForbidDR = \left\{ \begin{array}{l} (i, k, l, j) \mid i, j \in p, i \neq j \land k, l \in q, k \neq l \\ (k, i, j, l) \mid \forall (p, q) \in PDirF \end{array} \right\} $$

$$ ForbidIR = \left\{ \begin{array}{l} (i, k, j, l) \mid i, j \in p, i \neq j \land k, l \in q, k \neq l \\ (k, i, l, j) \mid \forall (p, q) \in PInvF \end{array} \right\} $$

Figure 36 illustrates the authorised and forbidden positions for DRP and IRP.

To know if we are in the forbidden cases described in these two sets, we propose to compare the vertices two-by-two. Denote by AlphaDR and AlphaIR the sets containing the couples of vertices to be compared to determine the forbidden cases
<figure><img src="image_84.png" /><figcaption>■ Figure 36 – Non-exhaustive illustrations for authorised and forbidden order cases for two repeated fragments $((i, j), (k, l)) \in PRepF$.</figcaption></figure>

(a) and (b) Authorised and forbidden orders for $PDirF$; (c) and (d) Authorised and forbidden orders for $PInvF$.

respectively associated with ForbidDR and ForbidIR sets, such that:

$$ AlphaDR = \left\{ \begin{array}{l} (i, j), (k, l), (i, k), \\ (j, l), (i, l), (j, k) \\ \text{s.t. } ((i, j), (k, l)) \in PDirF \end{array} \right\} $$

$$ AlphaIR = \left\{ \begin{array}{l} (i, k), (i, l), (j, k), (j, l) \\ \text{s.t. } ((i, j), (k, l)) \in PInvF \end{array} \right\} $$

In the following, the sets of repeated fragments and these for the forbidden orders are abstracted to generalise DRP and IRP. Table 12 gives the correspondence of the sets depending on the problem to solve.

<table><caption>■ Table 12 – ILP sets and functions corresponding table.</caption><thead><tr><th>ILP</th><th>RepF</th><th>PRepF</th><th>ARepF</th><th>Forbid</th><th>Alpha</th><th>repfrag</th><th>repadj</th></tr></thead><tbody><tr><td>DRP</td><td>DirF</td><td>PDirF</td><td>ADirF</td><td>ForbidDR</td><td>AlphaDR</td><td>dirfrag</td><td>diradj</td></tr><tr><td>IRP</td><td>InvF</td><td>PInvF</td><td>AInvF</td><td>ForbidIR</td><td>AlphaIR</td><td>invfrag</td><td>invadj</td></tr></tbody></table>

Binary variables

— $m_{ij}$ encodes if the repeated fragment $(i, j) \in RepF$ is a part of a repeat.
— isadj$_e$ encodes if two repeated fragments connected by the edge $e \in ARepF$ (and repadj($e$) $\in E$) are adjacent in the circuit.

— forbid$_{ijkl}$ encodes whether we are in the forbidden vertices order $(i, j, k, l) \in Forbid$.

— $\alpha_{uv}$ encodes whether the vertex $u$ is before the vertex $v$ in the circuit. Since $\alpha_{uv} = 1 - \alpha_{vu}$, even if $(v, u) \notin Alpha$, for clarity we write $\alpha_{vu}$ instead of $1 - \alpha_{uv}$.

## Continuous variables

— $i_v \in [0, 1]$ encodes if the vertex $v \in V \setminus \{s, \bar{s}\}$ participates in the circuit, and acts as a binary variable.

— $f_e \in \mathbb{R}_{\ge 0}$ is the positive flow on the participating edge $e \in E$ in the circuit (zero otherwise). We use the exiting flow to define the position $pos(v)$ of a vertex $v \in V$, i.e. $pos(v) = \sum_{w \in N_v^+} f_w$.

The vertices of participating repeated fragments must be in the circuit (Constraints C9 and C10). Constraints C11 to C13 implement with linear constraints the $\alpha_{uv}$ definition. Constraints C14 to C16 implement the forbid$_{ijkl}$ definition. Constraints C17 to C20 implement the isadj$_e$ definition.

## CRepeat constraints

Add the set of constraints CCircuit

$$ \forall (i, j) \in RepF : \qquad m_{ij} \le i_i \qquad (C9) $$

$$ m_{ij} \le i_j \qquad (C10) $$

$$ \forall (u, v) \in Alpha : \qquad pos(v) - pos(u) \le M\alpha_{uv} \qquad (C11) $$

$$ pos(u) - pos(v) \le M(1 - \alpha_{uv}) \qquad (C12) $$

$$ pos(u) + pos(v) \ge \alpha_{uv} \qquad (C13) $$
$$ \forall (i, j, k, l) \in \text{Forbid}: \qquad 3\text{forbid}_{ijkl} \le \alpha_{ij} + \alpha_{jk} + \alpha_{kl} \qquad (C14) $$

$$ 2 + \text{forbid}_{ijkl} \ge \alpha_{ij} + \alpha_{jk} + \alpha_{kl} \qquad (C15) $$

$$ \forall (p, q) \in \text{PRepF}: \qquad m_p + m_q \le 2 - \sum_{\substack{(i,j,k,l) \\ \in \text{Forbid} \\ \text{s.t. } (p,q)}} \text{forbid}_{ijkl} \qquad (C16) $$

$$ \forall (u, v) \in \text{ARepF}: \qquad \text{isadj}_{uv} \le x_{uv} \qquad (C17) $$

$$ \text{isadj}_{uv} \le x_{\text{repadj}(u,v)} \qquad (C18) $$

$$ \text{isadj}_{uv} \le m_{\text{repfrag}(u)} \qquad (C19) $$

$$ \text{isadj}_{uv} \le m_{\text{repfrag}(v)} \qquad (C20) $$

$$ m_p \in \{0, 1\} \quad \forall p \in \text{RepF} $$

$$ \text{isadj}_e \in \{0, 1\} \quad \forall e \in \text{ARepF} $$

$$ \text{forbid}_{ijkl} \in \{0, 1\} \quad \forall (i, j, k, l) \in \text{Forbid} $$

$$ \alpha_{uv} \in \{0, 1\} \quad \forall (u, v) \in \text{Alpha} $$

## 5.3 Fixing regions constraints

When repeats are previously scaffolded, the involved regions are fixed as input for the next problems. Let $ADirF^*$, $DirF^*$, $AInvF^*$ and $InvF^*$ respectively be the sets of (adjacent) direct and (adjacent) inverted fragments composing the direct and inverted repeats that have been scaffolded.

### CFixRegions constraints

$$ \forall (u, v) \in \text{DirF}^* \cup \text{InvF}^*: \qquad i_u = 1 \qquad (C21) $$

$$ i_v = 1 \qquad (C22) $$

$$ \forall (u, v) \in \text{ADirF}^* \cup \text{AInvF}^*: \qquad x_{uv} = 1 \qquad (C23) $$
$$x_{diradj(u,v)} = 1 \quad \forall (u, v) \in ADirF^* \qquad (C24)$$

$$x_{invadj(u,v)} = 1 \quad \forall (u, v) \in AInvF^* \qquad (C25)$$

## 5.4 Speed-up constraints

Constraints C26 and C27 prevent the solver to loop on strictly equivalent solutions, e.g. solutions that differ according to a permutation of the occurrences. Denote by *ConsOcc* the set of occurrence-consecutive vertices, such that:

$$ConsOcc = \left\{ (u, v) \in V^2 \mid \begin{array}{l} ctg_u = ctg_v \\ \land or_u = or_v = f \\ \land occ_u = occ_v - 1 \end{array} \right\}$$

Also, denote by *ConsRepF* the set of consecutive repeated fragments, such that:

$$ConsRepF = \left\{ \begin{array}{l} ((i, j), (k, l)) \in RepF^2 \text{ s.t.} \\ ctg_i = ctg_j = ctg_k = ctg_l \\ \land or_i = or_k \land or_j = or_l \\ \land occ_i = occ_k - 2 \\ \land occ_j = occ_l - 2 \end{array} \right\}$$

$$i_v + i_{\bar{v}} \le i_u + i_{\bar{u}} \qquad \forall (u, v) \in ConsOcc \qquad (C26)$$

$$m_q \le m_p \qquad \forall (p, q) \in ConsRepF \qquad (C27)$$

## 5.5 Scaffolding problems ILP

Finally, it is possible to define the ILP formulations for the *DRP*, *IRP* and *SCP* scaffolding problems as a union of the constraints described before.

For *DRP* and *IRP* the ILP formulations are the same, and it is sufficient to choose the sets *RepF*, *PRepF*, *ARepF*, *Alpha*, *Forbid* and *ConsRepF* according to the repeats the problems scaffold. We aim to maximise the cumulative length of the minimum number of repeats. The objective value corresponds to:

$$\sum_{r \in Repeats} replen(r) - |Repeats| \\ = \sum_{p \in RepF} 2m_p - \left( \sum_{p \in RepF} m_p - \sum_{e \in ARepF} isadj_e \right) \\ = \sum_{p \in RepF} m_p + \sum_{e \in ARepF} isadj_e$$

where *Repeats* is the set of repeats.
DRP and IRP models

$$ \max \sum_{p \in \text{RepF}} m_p + \sum_{e \in \text{ARepF}} isadj_e $$

s.t. CCircuit

CRepeat

(C26)

(C27) if no repeats to fix

CFixRegions otherwise

Traditionally, SCP finds the maximum weighted circuit.

SCP model

$$ \max \sum_{v \in V \setminus \{s, \bar{s}\}} wex_v i_v $$

s.t. CCircuit

(C26)

CFixRegions if repeats to fix

Both for DRP and IRP, the number of variables and constraints are in $O(|V|^2 + |E|)$. The number of variables and constraints for SCP are in $O(|V| + |E|)$.

# 6 Hierarchical problem succession

Finally, here we describe how we combine the DRP, IRP and SCP scaffolding problems. As described in Definition 3.6, two problem combinations are opposed. The combinations are kept depending on the value of the problems' objective functions.

## ► Definition 6.1: Hierarchical problem succession solutions

Denote by $h_1, h_2$ the two hierarchical problem successions DRP-IRP-SCP and IRP-DRP-SCP. For each $h \in \{h_1, h_2\}$, denote by $Z_h^* \in \mathbb{R}_{\ge 0}^3$ the vector containing the values of the objective functions for each problem in the order of the problem succession corresponding to $h$.
By $S$ we denote the set of optimal problem successions, such that:

$$S = \left\{ s \mid \forall k \in [0, 2], Z_s^*[k] = \max_{h \in \{h_1, h_2\}} Z_h^*[k] \right\}$$

Note that $0 \le |S| \le 2$, and Definition 6.1 is stable for any problem with an objective value equal to zero. For example, $Z_{h_1}^*[1] = 0$ means that there is no inverted repeat. For an easier interpretation of the architecture, we adopt a problem code combination, summarised in Table 13.

<table><caption>Table 13 – Problem code combinations</caption><thead><tr><th></th><th>Z<sub>h</sub><sup>*</sup></th><th>h<sub>1</sub></th><th>h<sub>2</sub></th></tr><tr><th></th><th>[0]</th><th>[1]</th><th>[2]</th><th>DRP-IRP-SCP</th><th>IRP-DRP-SCP</th></tr></thead><tbody><tr><th>0</th><td>0</td><td>-</td><td>sc</td><td>sc</td><td>sc</td></tr><tr><th>&gt; 0</th><td>0</td><td>-</td><td>dr-sc</td><td>ir-sc</td><td>ir-sc</td></tr><tr><th>&gt; 0</th><td>&gt; 0</td><td>-</td><td>dr-ir-sc</td><td>ir-dr-sc</td><td>ir-dr-sc</td></tr></tbody></table>

## 7 From an ILP solution to a genome structure

From a solution found by the best hierarchical problem succession, we extract the corresponding genome architecture. Let $m, n \in \mathbb{N}$ be the number of repeats (pair of regions) and the number of single-copies, respectively. A genome contains $2m + n$ regions. Two items are sufficient to describe a genome with its regions:

- COR = $(r_0, r_1, \dots, r_{m+n-1})$ contains the forward regions, i.e. it is a $(m+n)$-uplet of $m+n$ oriented contig sequences.

- SOR = $((rid_0, ror_0), \dots, (rid_{2m+n-1}, ror_{2m+n-1})) \in (\mathbb{N} \times \{f, r\})^{2m+n}$ is a linearised circular sequence of oriented regions. For each $i \in [0, 2m+n[$, if $ror_i = f$ then SOR[i] represents the forward region COR[rid_i], else if $ror_i = r$ then SOR[i] represents the reverse $COR[rid_i]$.

Figure 37 illustrates the regions extracted by Algorithm 25 in the toy example.

The key idea of such an extraction algorithm is to start from the starting vertex $s$ and walk over the chosen edges. During the walk, for each vertex we need to identify the type of its region, and then add the vertex to the corresponding region (to the corresponding region identifier $rid_r$). Algorithm 22 aims to determine the
<figure><img src="image_85.png" /><figcaption>■ **Figure 37 – Extracting the genome architecture in MDCG from a CHSP solution.**</figcaption></figure>

The illustrated CHSP solution consists of the red circuit $(a_{f,0}, c_{r,1}, b_{f,0}, d_{f,0}, b_{r,1}, c_{f,0})$ and of the two choosen inverted fragments $InvF^* = \{(c_{f,0}, c_{r,1}), (b_{f,0}, b_{r,1})\}$ (vertex pairs linked by the green dashed lines). Algorithm 25 returns four regions: $m = 1$ inverted repeat (the two green arrows pointing downwards) and $n = 2$ single-copies (the purple and the red arrows). $COR = (a_f, c_r, b_f, d_f)$ and $SOR = (0_f, 1_f, 2_f, 1_r)$, where $x_y = (x, y)$ for clarity.

region type for a given vertex. The sequence of oriented regions *SOR* begins by the starter's region that is necessarily a single-copy (because *mult(contig(s))* = 1). The first oriented contig may not be the starter. Algorithm 23 gives the *initial vertex* associated to the first oriented contig of the starter's region. During the walk in the solution circuit from the initial vertex, a new region begins each time the region type changes. When the current vertex participates in a repeat, we must check if the next one participates in the same repeat, although the region type may not change (Algorithm 24). At the end, Algorithm 25 builds *COR* and *SOR* from an ILP solution.

To build the repeated regions, we use a First In First Out (FIFO) for the DRs, and a Last In First Out (LIFO) for the IRs. Each queue *rep_queue*, is associated with three methods:

*rep_queuePUT(x)* append *x* to the FIFO/LIFO;

*rep_queueIS_EMPTY* returns true if the queue is empty;

*rep_queuePEEK* returns the first/last value in the FIFO/LIFO;

*rep_queuePOP* deletes the first/last value in the FIFO/LIFO and returns it.

In the following, the given time complexities are under the assumption that the belonging test "$is x \in X?"$ for an object *x* in a set *X* is in $\Theta(1)$. Algorithm 22 is in $\Theta(1)$. Algorithm 23 is in $O(|SC_s|)$, where $|SC_s|$ is the number of contigs composing the single-copy region that contains the starting vertex. Algorithm 24 is in $\Theta(1)$, and so Algorithm 25 is in $O(|V| + |E|)$.
► Algorithm 22: Get the region type for a given vertex

**Require:** A vertex $v$, the chosen direct fragments $DirF^*$ and the chosen inverted fragments $InvF^*$.

**Ensure:** Returns the region type of the vertex.

1: **function** GETRegionCode($v$)

2: if $dirfrag(v) \in DirF^*$ then

3:   **return** DR

4: if $invfrag(v) \in InvF^*$ then

5:   **return** IR

6:   **return** SC

► Algorithm 23: Get the initial vertex of the single-copy region containing
the starting vertex

**Require:** The starting vertex $s$, the set of fixed variable values $x_e^*, e \in E$, the chosen direct fragments $DirF^*$ and the chosen inverted fragments $InvF^*$.

**Ensure:** Returns the first vertex of the single-copy region containing the starting vertex.

1: **function** INITIALoinITEX()

2: $v \leftarrow s$

3: $u \leftarrow u \mid x_{us}^* = 1$

4: reg\_code $\leftarrow$ GET\_REGION\_CODE($u$)

5: while $u \neq s \land reg\_code = SC$ do

6:   $v \leftarrow u$

7: $u \leftarrow u \mid x_{uv}^* = 1$

8: reg\_code $\leftarrow$ GET\_REGION\_CODE($u$)

9: if $region\_code = SC$ then ▷ Special case of a unique single-copy region

10:   **return** $s$

11: return $v$

► Algorithm 24: Is the repeat contiguous?

**Require:** Two vertices $u, v \in V$, the set of fixed variable values $x_e^*, e \in E$, the previous and current region codes *prev_region_code* and *region_code*.

**Ensure:** Returns true if the repeat is contiguous, else false.

1: **function** REPEAT_IS_CONTINGUOUS($u, v, prev_region_code, region_code$)

2: if $prev_region_code \neq region_code$ then

3:   **return** False
4: **if** region_code = DR **then**

5: **return** x\*<sub>diradj(u,v)</sub> = 1

6: **if** region_code = IR **then**

7: **return** x\*<sub>invadj(u,v)</sub> = 1

8: **return** False

► Algorithm 25: ILP solution to regions

**Require:** Graph *MDCG*, the starting vertex s, the set of fixed variable values
$x_e^*, e \in E$, the chosen direct fragments *DirF*\* and the chosen inverted
fragments *InvF*\*.

**Ensure:** Returns the oriented region sequence *SOR*, and the oriented contig sequence for each forward oriented region *COR*.

1: **function** ILP_SOLUTION_TO Regions()

2: prev_region_code ← SC

3: SOR ← [0<sub>f</sub>] ▷ Oriented region sequence, initialised by the single-copy in
forward orientation that contains the starting vertex

4: COR ← [[]] ▷ Oriented contig sequence for each forward region, initialised
         for the first region

5: region_index ← 0 ▷ First region index

6: init_v ← INITIAL VERTEX()

7: u ← init_v; v ← init_v

8: rep_queue ← HASHTABLE() ▷ For each repeat index is associated a queue
         of vertices (FIFO for DR, LIFO for IR)

9: repf_rep ← HASHTABLE() ▷ Each repeated fragment is associated with
its repeat index

10: repeat

11: region_code ← GET REGION CODE(v)

12: **if** region_code = SC **then**

13: **if** prev_region_code ≠ region_code **then** ▷ New single-copy

14: region_index ← |COR|

15: SOR. APPEND(region_index<sub>f</sub>)

16: COR. APPEND([])

17: COR[region_index].APPEND(ctg<sub>v</sub>, or<sub>v</sub>)

18: **else**

19: repf ← dirfrag(v) if region_code = DR else invfrag(v)

20: **if** repf ∉ repf_rep **then** ▷ Fist region of the repeat

21: **if** ¬REPEAT_IS_CONTIGUOUS(u, v, prev_region_code,
region_code) **then**
22: region_index ← |COR|
23: SOR. APPEND(region_index_f)
24: COR. APPEND([])
25: rep_queue[region_index] ← FIFO() if region_code = DR
else LIFO()
26: COR[region_index].APPEND(ctg_v, or_v)
27: rep_queue. PUT(dfrag(v) if region_code =
DR else invfrag(v))
28: repf_rep[repf] ← region_index
29: else ▷ Second region of the repeat
30: if ¬REPEAT_IS_CONTINGUOUS (u, v, prev_region_code,
region_code) then
31: region_index ← repf_rep[repf]
32: SOR. APPEND(region_index_f if region_code =
DR else region_index_r)
33: assert v = rep_queue[region_index].POP()
34: prev_region_code ← region_code
35: u ← v
36: v ← v | x_{uv}^* = 1
37: until v ≠ init_v
38: return SOR, COR

8 Multiple genome forms

The sequence of oriented regions *SOR* represents one chloroplast genome form.
Recall that especially with the LSC-IR-SSC-IR architecture (Figure 5a), the SSC
can be reversed during the DNA replication phase. In the following, our goal is to
retrieve other forms from the one obtained by the hierarchical problem succession.

Towards this goal, we introduce a specific assembly graph: the *region graph*.
The discovery of multiple genome forms is associated with the search of Eulerian
circuits in this graph. Figure 38 illustrates the procedure for the toy example.

► **Definition 8.1: Region graph** *RegGraph*

Given the (*m* + *n*)-uplet *COR* of forward regions and the sequence *SOR* of
oriented regions, *RegGraph* = (*Vreg*, *Ereg*, *Φ*) denotes a multidigraph named
the region graph, such that:

$$
- Vreg = \{v_{0,f}, \dots, v_{m+n-1,f}, v_{0,r}, \dots, v_{m+n-1,r} \mid |COR| = m+n\} \text{ is the}
$$
<figure><img src="image_86.png" /><figcaption>■ **Figure 38 – Region graph for the toy example.**</figcaption></figure>

The toy example's solution given in Figure 37 results in the graph visualised in Figure 33 will give a region graph that shows an LSC-IR-SSC-IR architecture as the one given in Figure 5a. Two oriented region sequences (genome forms) are obtained by finding the Eulerian circuits: $0_f \to 1_f \to 2_f \to 1_r$ and $0_f \to 1_f \to 2_r \to 1_r$. They correspond to two structural haplotypes commonly found in the chloroplast cells, and described in Figure 6. (a) Each arrow represents a region. Each link connects two arrows' extremity. Entering the tail/head and exiting the head/tail of an arrow corresponds to choosing the region in its forward/reverse orientation. (b) The same information is visualised. Each vertex is a region with a fixed orientation, and each edge connects two oriented regions.

set of oriented regions ($|Vreg| = 2|COR|$). For each vertex, bijective function $vreg: Vreg \hookrightarrow \{r \in COR\} \times \{f, r\}$ provides the oriented region it represents.

- *Ereg* is the multiset of links between two oriented regions in *SOR* (including between the last and the first regions). Note that $\forall e \in Ereg, \bar{e} \in \bar{Ereg}$ denotes its reverse where $\Phi(e) = (u, v), \Phi(\bar{e}) = (\bar{v}, \bar{u})$ ($|Ereg| = 2|SOR|$).

- $\Phi: Ereg \to \{(u, v) | u, v \in Vreg\}$ is the incident function, such that for two consecutive oriented regions $r_i$ and $r_j$ in *SOR* (including the last and the first ones), $\exists! e \in Ereg | \Phi(e) = (vreg^{-1}(r_i), vreg^{-1}(r_j))$.

Figure 38 illustrates the region graph for the toy example's solution given in Figure 37. Based on the above graph, we can find different sequences of oriented regions. Each region, independently of its orientation, in each sequence, must participate the same number of times.
► **Definition 8.2: Eulerian circuit in RegGraph**

A circuit in RegGraph = (Vreg, Ereg, Φ) is defined as Eulerian when:

- it begins from and ends with vertex $0_f$ representing the region containing the starter in forward orientation, i.e. $vreg(0_f) = COR[0]$;

- it passes through exactly one of the two versions of each edge ($e \in Ereg$ otherwise $\bar{e} \in Ereg$).

► **Proposition 8.1: An eulerian circuit is a valid oriented contig sequence**
**for RegGraph**

An Eulerian circuit in RegGraph = (Vreg, Ereg, Φ) provides a valid sequence of oriented contigs (Definition 2.1).

▷ Proof

Let RegGraph = (Vreg, Ereg, Φ) be a region graph. Denote by euc = (v₀, v₁, ..., v2m+n-1) an Eulerian circuit in RegGraph. For each two consecutive oriented regions vi, vj in euc, there exists an edge e ∈ Ereg such that Φ(e) = (vi, vj). According to Definition 8.1, vreg(vi) = rᵢ and vreg(vj) = rⱼ are also consecutive in the oriented region sequence that has originally built RegGraph, otherwise in its reverse. Thus, according to Algorithm 25, there is an edge (u, v) ∈ E in MDCG such that (ctgᵢ, orᵢ) and (ctgⱼ, orⱼ) equal rᵢ[|rᵢ| - 1] and rⱼ[0], respectively.

<figure><img src="image_87.png" /></figure>

We can easily verify that the number of Eulerian circuits is bounded by $O(2^{m'})$, where $m' \le m$ is the number of inverted repeats.

Now, given a region graph RegGraph, finding all the Eulerian circuits is equivalent to retrieve all the possible chloroplast genome forms. Each Eulerian circuit traverses exactly the same regions, but not necessary in the same orientations. Figure 38 gives the resulting region graph obtained from the input data given in Table 11.

We accept all the Eulerian circuits, although they may contradict the repeated region interval constraints given in Definitions 3.2 and 3.3. For example, the oriented region sequence (0<sub>f</sub>, 1<sub>f</sub>, 2<sub>f</sub>, 3<sub>f</sub>, 1<sub>f</sub>, 3<sub>r</sub>) respects the definitions, where (1<sub>f</sub>, 1<sub>f</sub>) is a DR and (3<sub>f</sub>, 3<sub>r</sub>) is an IR. One of the Eulerian circuit produces the oriented region sequence (0<sub>f</sub>, 1<sub>f</sub>, 2<sub>f</sub>, 3<sub>f</sub>, 1<sub>r</sub>, 3<sub>r</sub>). Region 1<sub>f</sub> now evolves in a new IR (1<sub>f</sub>, 1<sub>r</sub>). The order between the regions of the two IRs contradicts Definition 3.3.
9 NP-completeness

To prove the *NP*-completeness of one of the two hierarchical problem successions (decision version), it is sufficient to focus on only one scaffolding problem for the repeat. Here, we will focus on the decision version of *IRP* (*DIRP*).

► **Definition 9.1: IRP decision problem (DIRP)**

Given a set of contigs *C*, their multiplicities, a set of links *L*, a starting contig *s*, two integers *k*, *m*' ∈ N, is there a valid sequence of oriented regions for *IRP* with Σ<sub>ir∈IR</sub> replen(<i>ir</i>) ≥ *k* and |*IR*| ≤ *m*', where *IR* is the solution set of inverted repeats?

► **Proposition 9.1:*** *DIRP is in NP*

Given the input of *DIRP*, a sequence of oriented regions *SOR*, the sequence of oriented contigs for each region *COR*, two integers *k*, *m'* ∈ N. There is a polynomial time algorithm that checks if the given solution is valid and if its cumulative repeat length equals at least *k* and the number of repeats equals at most *m*'.

▷ Proof

Algorithm 26 verifies if the sequence of oriented regions is valid. It requires two traversals of SOR: (i) to identify which regions in the sequence form IRs; (ii) to verify if the order between the IRs is valid (thanks to the use of a LIFO). It also checks if the associated cumulative repeat length equals at least *k* and the number of repeats equals at most *m*'.

A LIFO *ir_lifo* is associated with four methods:

*ir_lifoPUT(x)* append *x* to the LIFO;

*ir_lifoIS_EMPTY* returns true if the LIFO is empty;

*ir_lifoPEEK* returns the last value in the LIFO;

*ir_lifoPOP* deletes the last value in the LIFO and returns it.

Algorithm 27 verifies if the corresponding sequence of oriented contigs is valid. It first retrieves the total sequence of oriented contig *SOC* from the sequence of oriented regions *SOR* and the sequence of each forward regions *COR*. It traverses the oriented contigs in *SOC*, verify if each two consecutive
oriented contigs are linked in the link set, and count the number of times a contig occur. At the end, it checks if each contig does not appear more than its multiplicity number of times.

It is straightforward to see that Algorithm 26 and Algorithm 27 are linear according the size of SOR and COR. Remember we assume that the belonging test "is $x \in X?"" for an object $x$ in a set $X$ is in $\Theta(1)$.

Δ

► Algorithm 26: Verify the validity of the sequence of oriented regions

**Require:** Input of *DIRP*, a sequence of oriented regions *SOR*, the sequence of oriented contigs for each region *COR*, an integer $k \in \mathbb{N}$.

**Ensure:** Returns True if the sequence of oriented regions is valid.

1: **function** `IS_SOR_VALID()`
2:     `IR` ← `{}`
3:     `count_reg` ← `HASHTABLE()`
4:     `count_reg` ← `Keep the orientations of a region in a hash table`
5:     `Identify the inverted repeats`
6:     **for** `(rid, ror) \in SOR` **do**
7:         **if** `rid` ∉ `count_reg` **then**
8:             **else**
9:                 **if** `count_reg[rid] = 2` **then**
10:                 **return** False
11:             **if** `ror` ≠ `count_reg[rid][0]` **then**
12:                 `IR ADD(rid)`
13:             `count_reg APPEND(ror)`
14:         **then**
15:     **if** $\sum_{ir \in IR} |COR[ir]| < k$ or $|IR| > m'$ **then**
16:         **return** False
17:     **then**
18:     **for** `(rid, ror) \in SOR` **do**
19:         **if** `rid` ∈ `IR` **then**
20:             **else**
21:                 **if** `ir_lifo.IS_EMPTY()` ∨ `rid` ≠ `ir_lifo.PEEK()` **then**
22:                 `ir_lifo PUT(rid)`
23:             **else**
24:                 `ir_lifo.POP()`
25:     **return** `ir_lifo.IS_EMPTY()`
► **Algorithm 27: Verify the validity of the sequence of oriented contigs**

**Require:** Input of *DIRP*, a sequence of oriented regions *SOR*, the sequence of oriented contigs for each region *COR*.

**Ensure:** Returns True if the sequence of oriented contigs is valid.

1: **function** `IS_SOC_VALID()`

2: Transform *SOR* in the sequence of oriented contigs *SOC* using *COR*.

3: `contig_count` ← `HASHTABLE()`

4: **if** *SOC*[0] ≠ (*s*, *f*) **then** ► *The circuit must start with the starting contig in forward orientation*

5: **return** False

6: ► *Verify the links and count the occurrences* ◄

7: (*prev_c*, *prev_or*) ← *SOC*[|*SOC*| − 1]

8: ► *Count the contigs and verify the links* ◄

9: **for** (*c*, *or*) ∈ *SOC* **do**

10: **if** ((*prev_c*, *prev_or*), (*c*, *or*)) ∉ *L* **then**

11: **return** False

12: **if** *c* ∉ `contig_count` **then**

13: `contig_count[c]` ← 1

14: **else**

15: `contig_count[c]` ← `contig_count[c]` + 1

16: (*prev_c*, *prev_or*) ← (*c*, *or*)

17: ► *Verify the multiplicity* ◄

18: **for** *c* ∈ `contig_count` **do**

19: **if** `contig_count[c]` > `mult(c)` **then**

20: **return** False

► **Proposition 9.2: NP-hardness of *IRP* decision problem**

*DIRP* is *NP*-hard.

▷ Proof

By reduction from the longest path decision problem from vertex *s* to vertex *t* (*LPSTP*), known to be *NP-complete* (Schrijver, 2003).

Consider an instance $I \in LPSTP$ composed of a directed graph $G = (V, E)$, two vertices $s, t \in V$ and an integer $k \in \mathbb{N}$ (the hypothetical number of vertices between $s$ and $t$ in the longest path). We shall build an instance transform function $tf$ such that $I \in LPSTP \iff tf(I) \in DIRP$. Function
<figure><img src="image_88.png" /><figcaption>(a) Directed graph G</figcaption></figure>

<figure><img src="image_89.png" /><figcaption>(b) Directed graph $G'$</figcaption></figure>

■ **Figure 39** – From a digraph *G* for *LPSTP* to a digraph *G'* for *IRP*.

Bold red edges in both sub-figures correspond to the solution path for LPSTP and IRP problems, respectively. (a) $G_{V\{s,t\}}$ is the subgraph induced by the vertex set $V \setminus \{s,t\}$. As the longest path exits *s* and enters *t*, dashed edges do not participate in the solution. (b) Green dashed line between vertices in $G'_{f_0}$ and $G'_{r_1}$ visualise the inverted fragments.

$tf$ transforms graph $G$ to graph $G' = (V', E')$, vertex $s \in V$ to vertices $s_f, s'_f \in V', t \in V$ to $t_f, t_r \in V', k$ to $k' = 2k$, and fix parameter $m' = 1$. Figure 39 illustrates the transformation.

All four subgraphs $G'_{or,i}$ in Figure 39b can be seen as copies of $G_{V\{s,t\}} = (V \setminus \{s,t\}, E_{V\{s,t\}})$ (the subgraph induced by the vertex set $V \setminus \{s,t\}$) in Figure 39a. For each $or \in \{f,r\}$, for each $i \in \{0,1\}$, $G'_{or,i} = (V'_{or,i}, E'_{or,i})$ such that:

- there is a bijective function $vtrans_{or,i}: V\{s,t\} \hookrightarrow V_{or,i}$, where $\forall v \in V'_{or,i}$, $vor(v) = or$ and $vocc(v) = i$;

- there is a bijective function $etrans_{or,i}: E_{V\{s,t\}} \hookrightarrow E'_{or,i}$ where $\forall (u,v) \in E_{V\{s,t\}}$, $\exists! (u',v') \in E'_{or,i}$ such that:

$$vtrans_{or,i}(u), vtrans_{or,i}(v) = \begin{cases} u', v' & \text{if } or = f \\ v', u' & \text{if } or = r \end{cases}$$

There is a bijective function $vsttrans_{st}: \{s,t\} \hookrightarrow \{(s_f, s'_f), (t_f, t_r)\}$. There is a function $esttrans: \{(s,w) \in E\} \cup \{(u,t) \in E\} \to E'$ such that:
— $\forall (s, w) \in E$:

$$ esttrans(s, w) = \{(s_f, vtrans_{f,i}(w)) \in V'^2 | i \in \{0, 1\} \} \\ \cup \{(vtrans_{r,i}(w), s'_f) \in V'^2 | i \in \{0, 1\} \} $$

— $\forall (u, t) \in E$:

$$ esttrans(u, t) = \{(vtrans_{f,i}(u), t_f), (t_f, vtrans_{r,i}(u)) \in V'^2 | i \in \{0, 1\} \} \\ \cup \{(vtrans_{f,i}(u), t_r), (t_r, vtrans_{r,i}(u)) \in V'^2 | i \in \{0, 1\} \} $$

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

The sets $InvF$, $PInvF$ and $AInvF$ are built based on $G'_{or,i}$ graphs. Inverted fragments are visualised by green dashed vertical lines in Figure 39b, where $InvF = \{(i, j) \in V_{f,0} \times V_{r,1} | vtrans_{f,0}(i) = vtrans_{r,1}(j)\}$.

As the adjacent inverted fragments associate only vertices in $V_{f,0}$ with those in $V_{r,1}$, the path that maximises the number of contiguous inverted fragments exits $s_f$, goes through $G'_{f,0}$ to $t_f$ (or $t_r$, it does not matter), and passes through $G'_{r,1}$ to $s'_f$.

Since $G'_{f,0}$ is a copy of $G_{V\setminus\{s,t\}}$, while $G'_{r,1}$ is its reverse graph, there is a bijection between $V'_{f,0}$ and $V'_{r,1}$ vertices sets.

The way $G'$ is built implies only one IR is assembled (so parameter $m' = 1$ is respected). The length of this IR ($k' = 2k$) gives the hypothetical length of the longest path in $LPSTP$.

To conclude, as there is a linear time complexity transform function $tf$ such that $I \in LPSTP \iff tf(I) \in DIRP$, $DIRP$ is at least $NP$-hard.

Δ

From Propositions 9.1 and 9.2 we conclude that $DIRP$ is $NP$-complete.

# 10 Numerical results

We develop `khloraascaf`<fnref n="2" />, a python package that computes the scaffolding of chloroplast contigs. It can either use `Gurobi` solver or `CBC`. All the following runs have been executed on a Linux laptop computer (32GB RAM, Intel® Core™ i7-10610U CPU @ 1.80GHz ×8). Each time, we select `Gurobi` solver.

<footnote n="2">²https://khloraa-scaffolding.readthedocs.io/en/latest</footnote>
## 10.1 Complexity validation on artificial data

khloraascaf is also accessible as an API, that permits in this section to study the combinatorial behaviour of *IRP*.

We demonstrate in Section 9 that *DIRP* is in the general case *NP*-complete. Furthermore, the heteroplasmy for the chloroplast genome is very often caused by the presence of an inverted repeat, that reverses the region(s) between it.

Thus, we artificially build a contig set with the associated attributes, and a link set, such that the genome architecture behind corresponds to the following circular sequence of oriented regions: $SC1 - IR - SC2 - \overline{IR}$. In the following, we run what corresponds to *IRP* computation in khloraascaf on two types of growing generated data: perfect and noisy artificial ones. To emphasise the effect of the inverted repeats in the complexity of *IRP*, we fix the length of the single-copies, i.e. $|SC1| = |SC2| = |SC| = 20$, and we incrementally raise the length of the inverted repeat $|IR| = 20k$ for $k \in [1, 10]$.

### 10.1.1 Perfect artificial data

The data generated for this section correspond to the smallest set of contigs, links, and the smallest multiplicities to make sure that *IRP* is feasible. The multiplied doubled contig graph associated with these perfect data has exactly the same topology as the one illustrated in Figure 33. For instance, testing perfect artificial data acts as a control for further tests. Table 14 gives some Gurobi metrics<fnref n="3" />.

Observe that the gap is equal to 0% and the problem is solved either during the presolving or the linear relaxation. Indeed, for the class of graph that contains only the perfect artificial data, there exist a polynomial algorithm. However, the relaxation time seems to fit an exponential distribution as well as for the B&B time, as shown in Figure 40, even though the distributions should be treated with caution because of the limited number of points.

### 10.1.2 Noisy artificial data

Here we test the behaviour of the solver when we add noise to the perfect artificial data. In that case, for each generated contig, the multiplicity has 25% chance to be overestimated by one (that increases *InvF*, *PInvF* and possibly *AInvF* sets). Similarly, for each contig, there is 25% of chance to create a new link to another randomly chosen contig. This can generate more loops, and can increase

<footnote n="3">³To reproduce the results, please refer to https://khloraascaf-results.readthedocs.io/en/latest/benchmark_4/.</footnote>
<figure><img src="image_90.png" /><figcaption>■ **Figure 40 – Solver running time distributions for perfect artificial data.**</figcaption></figure>

Points are measured times, the red curves correspond to the best $ae^{bx}$ function applied on IR length axis.

the AInvF set<fnref n="4" />.

As expected, now the gap is not ensured to be null, and some instances are not solved at the presolving or at the linear relaxation steps.

These numerical results corroborate the *NP*-complete demonstration for a more general class of graphs.

## 10.2 Synthetic chloroplast input data

In this section, we aim to validate experimentally the relevance of our scaffolding problem definition by running khloraascaf on synthetic data.

### 10.2.1 Input data generation

Here we briefly describe our protocol for input data generation<fnref n="5" />. 200 chloroplast genomes (selected in CpGDB<fnref n="6" />) were downloaded from the NCBI<fnref n="7" />. For each of them,

<footnote n="4">⁴To reproduce the results, please refer to https://khloraascaf-results.readthedocs.io/en/latest/benchmark_5/.</footnote>

<footnote n="5">⁵For more details, please refer to https://khloraascaf-results.readthedocs.io/en/latest/benchmark_3/</footnote>

<footnote n="6">⁶CpGDB: A Comprehensive Database of Chloroplast Genomes http://www.gndu.ac.in/CpGDB/index.aspx</footnote>

<footnote n="7">⁷NCBI: The National Center for Biotechnology Information advances science and health by providing access to biomedical and genomic information https://www.ncbi.nlm.nih.gov/</footnote>
■ Table 14 – Gurobi solver metrics on perfect artificial growing data.

$|V|$, $|E|$, $|SC|$, $|IR|$ and $|\mathcal{L}|$ respectively stand for the number of vertices, edges, contigs in each single-copy region, contigs in each region of the inverted repeat, and links; **Time**: the presolve time plus the relaxation time (above) and the B&B time (below); **Opt.**: the linear relaxation bound $UB$ (above) and the integer optimal value $Opt$ (below); **% Gap**: the MIP gap equals $100 \times \frac{(UB-Opt)}{UB}$; **Nodes**: number of explored B&B nodes; **Iter.**: number of iterations for the LP relaxation (above) and for the B&B phase (below).

<table><thead><tr><th>V</th><th>E</th><th>SC</th><th>IR</th><th>L</th><th>Time</th><th>Opt.</th><th>% Gap</th><th>Nodes</th><th>Iter.</th></tr></thead><tbody><tr><td>160</td><td>244</td><td>20</td><td>20</td><td>122</td><td>0.25</td><td>39.00</td><td>0.00</td><td>1</td><td>1242</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.48</td><td>39.00</td><td></td><td></td><td>2944</td></tr><tr><td>240</td><td>404</td><td>20</td><td>40</td><td>162</td><td>1.23</td><td>79.00</td><td>0.00</td><td>1</td><td>4690</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>2.54</td><td>79.00</td><td></td><td></td><td>10 946</td></tr><tr><td>320</td><td>564</td><td>20</td><td>60</td><td>202</td><td>4.27</td><td>119.00</td><td>0.00</td><td>1</td><td>8945</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>10.45</td><td>119.00</td><td></td><td></td><td>23 864</td></tr><tr><td>400</td><td>724</td><td>20</td><td>80</td><td>242</td><td>12.15</td><td>159.00</td><td>0.00</td><td>1</td><td>15 196</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>26.29</td><td>159.00</td><td></td><td></td><td>39 562</td></tr><tr><td>480</td><td>884</td><td>20</td><td>100</td><td>282</td><td>13.46</td><td>199.00</td><td>0.00</td><td>1</td><td>23 109</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>66.28</td><td>199.00</td><td></td><td></td><td>68 740</td></tr><tr><td>560</td><td>1044</td><td>20</td><td>120</td><td>322</td><td>22.91</td><td>239.00</td><td>0.00</td><td>1</td><td>44 048</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>110.98</td><td>239.00</td><td></td><td></td><td>97 646</td></tr><tr><td>640</td><td>1204</td><td>20</td><td>140</td><td>362</td><td>35.68</td><td>279.00</td><td>0.00</td><td>1</td><td>46 071</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>80.26</td><td>279.00</td><td></td><td></td><td>125 084</td></tr><tr><td>720</td><td>1364</td><td>20</td><td>160</td><td>402</td><td>56.89</td><td>319.00</td><td>0.00</td><td>1</td><td>75 371</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>376.74</td><td>319.00</td><td></td><td></td><td>256 831</td></tr><tr><td>800</td><td>1524</td><td>20</td><td>180</td><td>442</td><td>321.52</td><td>359.00</td><td>0.00</td><td>1</td><td>71 971</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>634.01</td><td>359.00</td><td></td><td></td><td>196 905</td></tr><tr><td>880</td><td>1684</td><td>20</td><td>200</td><td>482</td><td>488.74</td><td>399.00</td><td>0.00</td><td>1</td><td>88 157</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>1458.66</td><td>399.00</td><td></td><td></td><td>236 406</td></tr></tbody></table>

a set of reads was generated. The contigs were generated with Minia (Chikhi and Rizk, 2012), a de Bruijn graph assembly approach. The links correspond to k-mer paths in the resulting compacted de Bruijn graph (cDBG) that connect two contigs. Finally, 31 instances were selected for which various difficulties have been detected, e.g. extra-links in the link set or combination of repeats.

The starter is the contig for which the matK gene, usually found in a single-copy region, maps on.

To obtain the multiplicity of a given contig $c$, we sum the length of the alignments of the reads mapping on $c$ (we denote by $MapR_c$ the set containing the reads that map on $c$). This sum is defined as the coverage $cov_c$ of the contig $c$ by the reads:

$$cov_c = \sum_{read \in MapR_c} |align_c^{read}|$$
■ Table 15 – Gurobi solver metrics on noisy artificial growing data.

The column descriptions are the same as the one in Table 14. Except that because of the noise on multiplicities, the sum of contig multiplicity for each region can change: this is the value below the number of contig in columns |SC1|, |SC2| and |IR|. Similarly, because of the noise on the number of links, the below value for $\mathcal{L}$ corresponds to the number of noisy links.

<table><thead><tr><th>|V|</th><th>|E|</th><th>SC1</th><th>SC2</th><th>IR</th><th>l</th><th>Time</th><th>Opt.</th><th>% Gap</th><th>Nodes</th><th>Iter.</th></tr></thead><tbody><tr><td rowspan="2">186</td><td rowspan="2">372</td><td>20</td><td>20</td><td>20</td><td>152</td><td>1.01</td><td>39.00</td><td>0.00</td><td rowspan="2">1</td><td>2862</td></tr><tr><td>26</td><td>24</td><td>43</td><td>30</td><td>5.16</td><td>39.00</td><td></td><td>6888</td></tr><tr><td rowspan="2">280</td><td rowspan="2">688</td><td>20</td><td>20</td><td>40</td><td>212</td><td>6.89</td><td>79.00</td><td>0.00</td><td rowspan="2">1</td><td>6727</td></tr><tr><td>28</td><td>23</td><td>89</td><td>50</td><td>12.26</td><td>79.00</td><td></td><td>19930</td></tr><tr><td rowspan="2">366</td><td rowspan="2">946</td><td>20</td><td>20</td><td>60</td><td>262</td><td>42.52</td><td>123.50</td><td>3.64</td><td rowspan="2">1</td><td>16821</td></tr><tr><td>25</td><td>26</td><td>132</td><td>60</td><td>79.28</td><td>119.00</td><td></td><td>50246</td></tr><tr><td rowspan="2">452</td><td rowspan="2">1208</td><td>20</td><td>20</td><td>80</td><td>320</td><td>90.06</td><td>161.50</td><td>1.55</td><td rowspan="2">1</td><td>27822</td></tr><tr><td>22</td><td>23</td><td>181</td><td>78</td><td>295.22</td><td>159.00</td><td></td><td>129174</td></tr><tr><td rowspan="2">556</td><td rowspan="2">1366</td><td>20</td><td>20</td><td>100</td><td>322</td><td>196.04</td><td>199.00</td><td>0.00</td><td rowspan="2">1</td><td>42292</td></tr><tr><td>23</td><td>24</td><td>231</td><td>40</td><td>244.41</td><td>199.00</td><td></td><td>92009</td></tr><tr><td rowspan="2">662</td><td rowspan="2">1804</td><td>20</td><td>20</td><td>120</td><td>412</td><td>1007.59</td><td>242.50</td><td>1.44</td><td rowspan="2">1</td><td>84639</td></tr><tr><td>29</td><td>26</td><td>276</td><td>90</td><td>1434.16</td><td>239.00</td><td></td><td>210798</td></tr><tr><td rowspan="2">736</td><td rowspan="2">1946</td><td>20</td><td>20</td><td>140</td><td>454</td><td>1108.09</td><td>283.00</td><td>1.41</td><td rowspan="2">1</td><td>228198</td></tr><tr><td>24</td><td>26</td><td>318</td><td>92</td><td>3540.79</td><td>279.00</td><td></td><td>691619</td></tr><tr><td rowspan="2">822</td><td rowspan="2">2212</td><td>20</td><td>20</td><td>160</td><td>514</td><td>1118.76</td><td>323.00</td><td>1.24</td><td rowspan="2">1</td><td>86592</td></tr><tr><td>26</td><td>27</td><td>358</td><td>112</td><td>2449.55</td><td>319.00</td><td></td><td>287146</td></tr><tr><td rowspan="2">902</td><td rowspan="2">2362</td><td>20</td><td>20</td><td>180</td><td>542</td><td>1591.18</td><td>363.00</td><td>1.10</td><td rowspan="2">1</td><td>91936</td></tr><tr><td>26</td><td>27</td><td>398</td><td>100</td><td>2576.60</td><td>359.00</td><td></td><td>269958</td></tr><tr><td rowspan="2">996</td><td rowspan="2">2656</td><td>20</td><td>20</td><td>200</td><td>602</td><td>2294.85</td><td>404.00</td><td>1.24</td><td rowspan="2">1</td><td>116315</td></tr><tr><td>26</td><td>26</td><td>446</td><td>120</td><td>3747.02</td><td>399.00</td><td></td><td>351501</td></tr></tbody></table>

Where $align_c^{read}$ is the sequence representing the alignment of the read *read* on the contig *c*. Its length equals the number of nucleotides of *read* that match on *c* (identity or substitution). Then the multiplicity $mult(c)$ of *c* is obtained by normalising its coverage $cov_c$ by this of the starter *s*, $cov_s$: $mult(c) = \max(1, [cov_c/cov_s - 0.1])$. As the multiplicity is an upper-bound of the usage of contig, we prefer to round up the normalisation only if the decimal part is greater than 0.1.

The existence-weight $wex(c)$ for a contig *c* is computed by counting the number of nucleotides of *c* that are covered by at least a gene of protein from a chloroplast near-species, normalised by the length of *c*.

### 10.2.2 The evaluation's metrics

For each synthetic instance, we know the sequence of the oriented contigs and the sequence of the oriented regions we search for. In the sequel we test our scaffolding
approach and evaluate the obtained region graph. For each instance, for each
optimisation problem combination, we provide the following metrics:

- the total number of eulerian circuits in the region graph (genome forms);

- how many of them coincide with the sequence of oriented contigs we search for;

- how many of them coincide with the sequence of oriented regions we search for.

Evaluating the sequence of the oriented contigs is stringent. On the other hand, although a result can be evaluated as a false one, we can still retrieve the sequence of the chloroplast genome by applying an alternative sequence. As a consequence, for each instance we use Quast (Gurevich et al., 2013) to evaluate the sequences associated to each genome form. As the genome reference is known, Quast tries to find the minimum number of differences (relocation, inversion, indels) between the reference and the sequence we provide. Three metrics are chosen to evaluate the best genome form for each problem succession:

- the genome fraction of the reference;

- the number of misassemblies;

- the number of local misassemblies.

For more detailed descriptions of Quast metrics, you can refer to Appendix 3.1.

It would be expected that Quast metrics illustrate wrong assembly for the
instances for which the sequence of the contigs, and a fortiori this of the regions,
are not retrieved. Analogously, the instance that truthfully retrieves the sequences
would have good Quast metric. However, these assertions may be contradicted
because of the contig and the link sequences generation.

In Section 10.2.3, we provide the two metric sets when khloraascaf is applied
to the original synthetic data, while in Section 10.2.4 the metrics are reported for
a subset of modified synthetic data.

10.2.3 Initial version

Table 16 provides all the metrics defined above for the 31 instances. The instances are solved very quickly (solver times < 4.5 sec., Table 19). khloraascaf successfully founds the sequence of the oriented contigs in 20 of them and retrieves the sequence of the oriented regions in 28 of them<fnref n="8" />. Three categories of failures are identified.

<footnote n="8">⁸For more details on the result for the initial version of the synthetic data, please refer to https://khloraascaf-results.readthedocs.io/en/latest/benchmark_3/v1_order_analysis/ and https://khloraascaf-results.readthedocs.io/en/latest/benchmark_3/v1_quast_analysis/.</footnote>
**Wrong starting contig and multiplicities estimations** This is the category for which our approach is dependent and sensible. In the presented version, we have used a given starting contig, and a wrong one can lead to reduce all the multiplicities. This is the case for Begonia_pulchrifolia and Lamprocapnos_spectabilis for which the starters are contigs that normally participate into an IR, that contradicts our assumption that the starter participates in an SC.

Independently of a right starting contig, the multiplicity computation is sensible from the noise on the contig coverage by the reads. Agathis_dammara and Pelargonium_nanum both suffer from only one contig under-estimated multiplicity.

*IRP's objective: maximising the cumulative length of the minimum number of repeats*
The sequence of oriented regions for the Cucumix_hystrix's reference contains the following sub-sequences: (..., *IR*, *IR'*, ..., *IR'* , *SC*, *IR*'). As *IRP* also aims to minimise the number of repeats, it results in the merge of *IR*–*IR'* and consequently does not retrieve *SC*, normally inserted between *IR'* and *IR*.

**Model's robustness** While the sequences of oriented contigs for the repeats are retrieved, the ones of the single-copies suffer from extra-links combined with low, sometimes null, weights. On the one hand, in case of null weights, the circuits in Lathyrus_pubescens and Triosteum_pinnatifidum do not pass through contigs that must participate in the SCs. On the other hand, Podocarpus_totara possesses two objective-equivalent subpaths in an SC.

Surprisingly, our tool khloraascaf reversed some subparts of single-copies. This is due to extra-links, but they are specifically caused by the existence of very short IRs hidden in them (remember that the links correspond to paths between pairs of contigs in the cDBG). This behaviour is observed in *Carpodetus_serratus*, *Jasminum_tortuosum* and *Lophocereus_schottii*.

**Nucleotide sequences misassemblies** To avoid confusion, we always use nucl. seq. to denote “nucleotide sequence” to contrast with sequences of contigs/regions. In the following we analyse the instances presenting an unexpected behaviour for the Quast metrics, regarding if the sequences are found. Supplementary Table 18 permits verifying if the contigs that participate in the sequence have been correctly assembled. Except for *Lathyrus_pubescens* where one local misassembly is due to the missing contig in the sequence, all the (local) misassemblies in *Commiphora_foliacea*, *EucommiaUlmoides*, *Juniperus_scopulorum*, *Musa_ornata*, *Sciadopitys Verticalis*, *Taxus_baccata*, *Welwitschia_mirabilis* provided by Quast are found in the nucl. seq. of links.
<table><caption>**Table 16 – Sequence and Quast metrics for the initial synthetic data version.**
For each instance in the column **Instance**: **ILPs** provides the optimal (at most two) hierarchical problem successions; **Total** reports the number of eulerian circuits in the region graph (genome forms); SOC is the number of oriented contig sequences that equal to the reference oriented contig sequence; SOR is the number of oriented region sequences in bijection with the reference oriented region sequence; **%gnm** is the genome fraction of the best sequence produced by one of the genome forms; **#mis** are the number of misassemblies (left) and of local misassemblies (right).</caption><thead><tr><th rowspan="2">Instance</th><th rowspan="2">ILPs</th><th rowspan="2">Total</th><th colspan="2">Successful</th><th colspan="2">Quast</th></tr><tr><th>SOC</th><th>SOR</th><th>%gnm</th><th>#mis</th></tr></thead><tbody><tr><td rowspan="2">Abies_alba</td><td>dr-sc</td><td>1</td><td>0</td><td>0</td><td>99.88</td><td>3</td><td>0</td></tr><tr><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Acorus_americanus</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.99</td><td>0</td><td>0</td></tr><tr><td rowspan="2">Agathis_dammara</td><td>dr-ir-sc</td><td>10</td><td>0</td><td>2</td><td>80.51</td><td>1</td><td>2</td></tr><tr><td>ir-dr-sc</td><td>6</td><td>0</td><td>0</td><td>80.94</td><td>2</td><td>2</td></tr><tr><td>Azima_tetracantha</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.99</td><td>0</td><td>0</td></tr><tr><td>Begonia_pulchrifolia</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Carpodetus_serratus</td><td>ir-sc</td><td>2</td><td>0</td><td>2</td><td>100.00</td><td>2</td><td>0</td></tr><tr><td>Circaeaster_agrestis</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.99</td><td>0</td><td>0</td></tr><tr><td>Clematis_repens</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Commiphora_foliacea</td><td>ir-sc</td><td>4</td><td>1</td><td>4</td><td>98.82</td><td>0</td><td>5</td></tr><tr><td>Cucumis_hystrix</td><td>ir-sc</td><td>2</td><td>0</td><td>0</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>EucommiaUlmoides</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.99</td><td>0</td><td>2</td></tr><tr><td>Jasminum_tortuosum</td><td>ir-sc</td><td>2</td><td>0</td><td>2</td><td>99.99</td><td>2</td><td>0</td></tr><tr><td>Juniperus_scopulorum</td><td>sc</td><td>1</td><td>1</td><td>1</td><td>99.72</td><td>0</td><td>1</td></tr><tr><td>Lamprocapnos_spectabilis</td><td>dr-sc</td><td>1</td><td>0</td><td>0</td><td>35.86</td><td>3</td><td>0</td></tr><tr><td rowspan="2">Lathyrus_pubescens</td><td>dr-sc</td><td>1</td><td>0</td><td>0</td><td>98.47</td><td>2</td><td>2</td></tr><tr><td>ir-sc</td><td>2</td><td>0</td><td>2</td><td>98.50</td><td>0</td><td>2</td></tr><tr><td>Lophocereus_schottii</td><td>sc</td><td>1</td><td>0</td><td>1</td><td>99.88</td><td>1</td><td>0</td></tr><tr><td>Musa_ornata</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.97</td><td>0</td><td>2</td></tr><tr><td>Oenothera_glazioviana</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Pelargonium_nanum</td><td>ir-sc</td><td>2</td><td>0</td><td>2</td><td>96.44</td><td>4</td><td>0</td></tr><tr><td>Podocarpus_totara</td><td>sc</td><td>1</td><td>0</td><td>1</td><td>99.71</td><td>4</td><td>0</td></tr><tr><td>Porphyra_purpura</td><td>dr-sc</td><td>1</td><td>1</td><td>1</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Sagittaria_trifolia</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td rowspan="2">SciadopitysVerticillata</td><td>dr-sc</td><td>1</td><td>0</td><td>0</td><td>98.99</td><td>6</td><td>3</td></tr><tr><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.00</td><td>1</td><td>4</td></tr><tr><td>Sciaphila_densiflora</td><td>sc</td><td>1</td><td>1</td><td>1</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Selaginella_kraussiana</td><td>dr-sc</td><td>1</td><td>1</td><td>1</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Selaginella_vardei</td><td>dr-sc</td><td>1</td><td>1</td><td>1</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Taxus_baccata</td><td>sc</td><td>1</td><td>1</td><td>1</td><td>99.77</td><td>0</td><td>2</td></tr><tr><td>Triosteum_pinnatifidum</td><td>ir-dr-sc</td><td>2</td><td>0</td><td>2</td><td>99.28</td><td>0</td><td>2</td></tr><tr><td>Uvaria_macrophylla</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.90</td><td>0</td><td>0</td></tr><tr><td>Welwitschia_mirabilis</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>100.00</td><td>0</td><td>1</td></tr></tbody></table>

(the table continues on the next page)
<table><caption>Table 16, continued</caption><thead><tr><th rowspan="2">Instance</th><th rowspan="2">ILPs</th><th rowspan="2">Total</th><th colspan="2">Successful</th><th colspan="2">Quast</th></tr><tr><th>SOC</th><th>SOR</th><th>%gnm</th><th>#mis</th></tr></thead><tbody><tr><td>Wolffia_australiana</td><td>ir-sc</td><td>4</td><td>1</td><td>4</td><td>99.99</td><td>0 0</td></tr></tbody></table>

### 10.2.4 Modified version

In this section, we present a manually changed synthetic data version to succeed the scaffolding. The goal is to precisely evaluate the robustness of *khloraascaf*. We detail the modifications below<fnref n="9" />:

* **Agathis_dammara** The multiplicity of contig 1 is raised to be equal to 3.

* **Begonia_pulchrifolia** Contig 4 becomes the starter. So according to the multiplicity computation described in Section 10.2.1, the multiplicities of contigs 1 to 5 become 1, while this one of contig 0 raises to 2.

* **Carpodetus_serratus** Link $(10_r, 11_f)$ is deleted.

* **Jasminum_tortuosum** Link $(6_f, 4_f)$ is added.

* **Lamprocapnosspectabilis** Contig 8 becomes the starter. So the multiplicities of contigs 2, 4, 6, 10 and 11 increase by one, while the ones of contigs 0, 1 and 3 increase by two.

* **Lathyrus_pubescens** The weight of contig 12 raises to 0.01.

* **Lophocereus_schottii** Link $(10_f, 3_f)$ is deleted.

* **Pelargonium_nanum** The multiplicity of contig 2 raises from 3 to 4, without respecting the computation of the multiplicity described in Section 10.2.1.

* **Podocarpus_totara** Link $(6_f, 11_r)$ is deleted.

* **Triosteum_pinnatifidum** The weight of contig 7 raises to 0.01.

Table 17 provides all the metrics obtained by running *khloraascaf*. The instances are also solved very quickly (solver times < 3 sec., Table 20). It truthfully finds all the sequences for the modified synthetic data except for *Agathis_dammara* instance. Although all the repeats (both the direct and the inverted ones) have been

<footnote n="9">⁹For more details, please refer to https://khloraascaf-results.readthedocs.io/en/lat est/benchmark_3/synthetic_data_v2/.</footnote>
retrieved, one of the single-copy region has not been found. It can be explained by extra-links that create alternative paths with the same optimal value for *SCP*. Note that the oriented region sequence has been still retrieved.

All the (local) misassemblies provided by Quast for *Carpodetus_serratus*, *Lathyrus_pubescens*, *Pelargonium_nanum* and *Triosteum_pinnatifidum* just concern the nucl. seq. of links that is not a *khloraascaf* issue<fnref n="10" />.

■ **Table 17 – Sequence and Quast metrics for the modified synthetic data version.**

The caption is the same as in Table 16.

<table><thead><tr><th rowspan="2">Instance</th><th rowspan="2">ILPs</th><th rowspan="2">Total</th><th colspan="2">Successful</th><th colspan="2">Quast</th></tr><tr><th>SOC</th><th>SOR</th><th>%gnm</th><th>#mis</th></tr></thead><tbody><tr><td rowspan="2">Agathis_dammara</td><td>dr-ir-sc</td><td>10</td><td>0</td><td>2</td><td>99.97</td><td>3</td><td>2</td></tr><tr><td>ir-dr-sc</td><td>6</td><td>0</td><td>0</td><td>99.99</td><td>2</td><td>2</td></tr><tr><td>Begonia_pulchrifolia</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.99</td><td>0</td><td>0</td></tr><tr><td>Carpodetus_serratus</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>100.00</td><td>0</td><td>1</td></tr><tr><td>Jasminum_tortuosum</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.99</td><td>0</td><td>0</td></tr><tr><td>Lamprocapnos_spectabilis</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.63</td><td>0</td><td>3</td></tr><tr><td rowspan="2">Lathyrus_pubescens</td><td>dr-sc</td><td>1</td><td>0</td><td>0</td><td>99.01</td><td>2</td><td>4</td></tr><tr><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>99.73</td><td>0</td><td>1</td></tr><tr><td>Lophocereus_schottii</td><td>sc</td><td>1</td><td>1</td><td>1</td><td>100.00</td><td>0</td><td>0</td></tr><tr><td>Pelargonium_nanum</td><td>ir-sc</td><td>2</td><td>1</td><td>2</td><td>97.66</td><td>1</td><td>2</td></tr><tr><td>Podocarpus_totara</td><td>sc</td><td>1</td><td>1</td><td>1</td><td>99.75</td><td>0</td><td>0</td></tr><tr><td>Triosteum_pinnatifidum</td><td>ir-dr-sc</td><td>2</td><td>1</td><td>2</td><td>99.53</td><td>0</td><td>2</td></tr></tbody></table>

11 Conclusion

While the scaffolding problem is traditionally defined with distances data between the contigs, we show it is possible to avoid them in the case of the well-studied circular chloroplast genomes. Based on their specificities, we provide a new scaffolding formulation focused on revealing structural haplotypes.

Under the assumption that chloroplast genomes possess few repeats, we formalise their architectures as combinations of direct and inverted repeats, joined by single-copies, where the repeats are couples of identical (or reversed) nucleotide sequences. We tackle the chloroplast genome scaffolding as a discrete optimisation problem

<footnote n="10">¹⁰For more details on the results for the cleared version of the synthetic data, please refer to https://khloraascaf-results.readthedocs.io/en/latest/benchmark_3/v2_order_analysis/ and https://khloraascaf-results.readthedocs.io/en/latest/benchmark_3/v2_quota_analysis/.</footnote>
that yields three suboptimisation ones. We split the inherent multi-objective
problem into one optimisation problem per region type. As a consequence, it
is necessary to choose the order of subproblem resolutions as a function of the
results of previously solved problems. This is what has been addressed through the
hierarchical combination strategy. We model each subproblem with an ILP.

As our dedicated chloroplast scaffolding definition is a region-scaffolding-driven,
the region graph is a natural data structure to reveal distinct genome forms that
can coexist in a same cell. Indeed, particularly due to an IR flip-flop mechanism,
regions between the IRs can be reversed during the genome replication process.

Moreover, we prove the decision version of the Chloroplast Scaffolding Problem (*CHSP*) to be *NP-complete* in the general case, even though numerical results on perfect artificial data suggest there is a class of *MDCG* graphs where the problem is in *P*. Without surprise, the noisy artificial data benchmark confirms the theoretical complexity.

We have implemented our approach and the ILP formulations in a Python3 package, **khloraascaf**<fnref n="11" />, that we test on synthetic chloroplast contigs and links data. When the input data permit finding the solution, **khloraascaf** successfully retrieves the genome forms. Even if the decision problem is *NP-complete*, the small size of the input data enables to quickly solve optimally *CHSP*.

Our results show that the scaffolding-repeat problem formulations *DRP* and *IRP* seem to be sufficient to scaffold the repeats. This tends to validate our assumptions on the small number of repeats, and especially on the sufficiency of defining the repeats as pairs of equal (reversed) nucleotide sequences.

**12** Discussion and perspectives

While our scaffolding problem formulation seems to be sufficient to retrieve the
repeats, it seems it is not fully suitable for single-copies. If we have applied
the maximum-weighted circuit problem to scaffold the single-copies after having
scaffolded the repeats, it was only with the intention of staying in the context of
global optimisation. On the one hand, having weights on links may have been more
appropriate than just considering weights on the contigs: in some sense, that is the
purpose of distances. On the other hand, khloraascaf could initially scaffold the
repeats and then propose several solutions to link them, e.g. scored by the weights.

Concerning the tests on synthetic data, we should use a more traditional
assembly graph input: in fact, as revealed by comparing the khloraascaf results
with the reference genomes, the used link generation suffers from local, or worse,
global, misassemblies. As next step, we plan to inject khloraascaf into a state-of-

<footnote n="11">¹¹https://pypi.org/project/khloraascaf/</footnote>
the-art chloroplast genome pipeline, like GetOrganelle, and substitute what can
be identified as the scaffolding part by our method. Hence, we should be able to
relevantly compare khloraascaf approach with the state-of-the-art.

khloraascaf is sensible to the contig multiplicity computation. For now, a
contig multiplicity is obtained by normalising its coverage by this of the starter. A
better strategy may be to choose the smallest coverage for the normalisation as we
expect the multiplicities to be upper-bounds of the contigs use.

Another issue concerns the choice of the starter: while it depends on the result
of the mapping of matK gene map on the contigs, DRP, IRP and SCP problems
may be adapted for a set of candidate for the starter.

To generalise *CHSP* on non-equally (reversed) pair of regions for the repeats,
two combining ideas are proposed: on the one hand, we can add pairs of contigs to
the repeated fragment sets from the user-input. On the other hand, *CHSP* should
be able to handle the case when a single-copy region is in only one of the regions
of a repeat: for now, the contiguity constraint and the objective exclude this case.
The contiguity constraint can be adapted to accept only one contiguous region out
of the two.

From a user-case perspective, the region graph data structure can be used to
determine what genome forms are present in the read set, and in which proportion.
Indeed, as the region graph explicitly describes the junctions between the regions,
especially between the inverted repeats, one may extract the nucleotide sequences
of these junctions to answer the existence of the forms in the read set.
<figure><img src="image_91.png" /></figure>

# CONCLUSIONS AND PERSPECTIVES

Michel Petrucciani. (1994). Estate [Song]. On *Live*. Blue Note

<figure><img src="image_92.png" /></figure>

<figure><img src="image_93.png" /></figure>

## In this chapter

<table><tr><td>Fragment graph</td><td>134</td></tr><tr><td>Thesis contribution</td><td>134</td></tr><tr><td>Short-term future work</td><td>135</td></tr><tr><td>Long-term future work</td><td>135</td></tr><tr><td>Chloroplast genome scaffolding</td><td>136</td></tr><tr><td>Thesis contribution</td><td>136</td></tr><tr><td>Short-term future work</td><td>137</td></tr><tr><td>Long-term future work</td><td>138</td></tr></table>

Since the 1980s, the fragment assembly's aims and their formulations have evolved following the development of sequencing technologies and the data characteristics they produce. Two main stages compose the fragment assembly process. The read assembly stage aims to assemble the reads into contigs based on sequence overlaps. Then, the scaffolding stage evolves orienting and ordering the contigs to form scaffolds. Finally, the scaffolds would represent the chromosomes with the minimum number of gaps and the minimum number of nucleotide differences.

The literature describes numerous formulations of combinatorial optimisation problems for the scaffolding problem, and as many approaches for dealing with large, noisy datasets. In this thesis, we outline the subsampling strategies, the combinatorial optimisation problem formulations and the solving approaches described in the state-of-the-art methods. We analyse and compare the fragment graph structures of the literature based on an implementation design. We then propose a formulation for the specific scaffolding of chloroplast genome, which we model with an Integer Linear Programming (ILP) and solve in Python3. In subsequent sections, we combine the conclusions of Chapters III and IV to discuss
their contribution to this thesis.

<figure><img src="image_94.png" /></figure>

# Fragment graph

Graphs are mathematical structures that are useful to store links between two fragments. Two main stages of the fragment assembly employ graphs to handle the input data or to output the result. In the read assembly stage, especially in the OLC approach, the fragments correspond to the reads and the links are the overlaps between two reads. In the scaffolding stage, the fragments represent the contigs and links can come from paired-end read or long read alignments against the contigs.

The literature describes three graph structures: a digraph, a bigraph and an ungraph. While the first one highlights visually the double-strand sequencing, its weakness is that it requires two vertices for each read (Kececioglu, 1991). Myers (1995) is the first to employ a bigraph to store overlaps. The key idea is to aggregate the two orientations of a read into only one vertex, a link and its reverse into only one edge. Finally, the ungraph structure associates one vertex for each of the two extremities of a fragment to simplify the graph traversal (Huson et al., 2002). The ungraph handles two sets of edges. A fragment-edge connects the two vertices corresponding to the two extremities of a fragment (the tail and the head). Two vertices representing the extremities of a fragment are connected by a fragment-edge. The second one is the multiset of link-edges. A link-edge represents both a link and its reverse.

## Thesis contribution

Although some authors have provided conceptual descriptions of these graph structures, the latter have never been compared before from an implementation stat point of view. In this thesis, we have proposed an implementation design based on adjacency lists for each graph structure. At the opposite of a compressed sparse row or column, the adjacency lists are more suitable for adding or deleting vertices or edges in the graphs. We have also described transformation processes to pass from one implementation to another. We then have visualised the graph structures and their implementations in a map.

We have retained three implementations: *DGS* and *DGF* for the multidigraph and *BGU* for the multibigraph. After we have compared their memory consumption, we theoretically measure the cost function for each of their algorithm on elementary operations, such as iterating over the neighbours, adding or deleting a vertex or an edge. We come to the conclusion that if memory is the critical issue, then
the *BGU* implementation should be preferred. Otherwise, we recommend *DGS* and *DGF* for iterating over the neighbours and for deleting a vertex or an edge. *BGU* is preferable for adding vertices or edges. To conclude, *DGF* proves to be well-balanced and ideally tailored for dynamic graph operations. *DGF* is available in a Python3 package named `revsymg`<fnref n="1" /> and is easily installable via `PyPI`.

## Short-term future work

We plan to implement *DGS*, *DGF* and *BGU* in Rust or an equivalent compiled programming language to compare them experimentally and confirm the theoretically obtained costs. Similarly, as the selection of a graph structure impacts further fragment assembly modelling, we aim to suggest integer linear models for path and vertex or edge coverage searches in each graph, analyse them, and compare them in theory and practice.

## Long-term future work

This thesis is the first to propose fragment graph implementations for comparison. The graphs are expressed as adjacency lists. It would be useful to propose implementations based on compressed matrices and compare them mutually and with adjacency lists.

However, these structures are suitable for RAM storage. If the size of the instances exceeds the RAM storage capacity, the entire or partial graph should remain on disk. An efficient disk storage strategy may benefit from using the disk or, even better, CPU caches.

The GFA<fnref n="2" /> and PAF<fnref n="3" /> file formats, which store assembly graphs and matches between pairs of oriented fragments, could serve as a source of inspiration for disk storage methods. They are reminiscent of bidirected or undirected graphs. One strategy might be to keep the whole graph in this form on disk and convert a part of interest into a *DGS* or *DGF* in RAM.

This strategy may rely on neighbourhood analyses or the specific topology of fragment graphs. It should therefore be related to the research on graph algorithms, partitioning and matrix computation.

Although this work illustrates the use of fragment graphs for the assembly process, other areas involving nucleotide sequences may benefit from this work. For instance, genome graphs only represent variation between genomes of individuals of the same species, but the presence of complementary reverse regions implies

<footnote n="1">¹https://pypi.org/project/revsymg/</footnote>

<footnote n="2">²https://github.com/GFA-spec/GFA-spec</footnote>

<footnote n="3">³https://github.com/lh3/miniasm/blob/master/PAF.md</footnote>
double-stranded behaviour. In this context, one could examine the adaptation of fragment graphs.

<figure><img src="image_95.png" /></figure>

# Chloroplast genome scaffolding

While the scaffolding problem is traditionally defined with distances data between the contigs, we show it is possible to avoid them in the case of the well-studied circular chloroplast genomes. Based on their specificities, we provide a new scaffolding formulation focused on revealing structural haplotypes.

## Thesis contribution

Under the assumption that chloroplast genomes possess few repeats, we formalise their architectures as combinations of direct and inverted repeats, joined by single-copies, where the repeats are couples of identical (or reversed) nucleotide sequences. We tackle the chloroplast genome scaffolding as a discrete optimisation problem that yields three suboptimisation ones. We split the inherent multi-objective problem into one optimisation problem per region type. As a consequence, it is necessary to choose the order of subproblem resolutions as a function of the results of previously solved problems. This is what has been addressed through the hierarchical combination strategy. We model each subproblem with an ILP.

As our dedicated chloroplast scaffolding definition is a region-scaffolding-driven, the region graph is a natural data structure to reveal the distinct genome forms that can coexist in a same cell. Indeed, particularly due to an IR flip-flop mechanism, regions between the IRs can be reversed during the genome replication process.

Moreover, we prove the Chloroplast Scaffolding Problem (CHSP) to be *NP*-complete in the general case, even though numerical results on perfect artificial data suggest there is a class of MDCG graphs where the problem is in *P*. Without surprise, the noisy artificial data benchmark confirms the theoretical complexity.

We have implemented our approach and the ILP formulations in a Python3 package, khloraascaf<fnref n="4" />, that we test on synthetic chloroplast contigs and links data. When the input data permit finding the solution, khloraascaf successfully retrieves the genome forms. Even if the decision problem is *NP*-complete, the small size of the input data enables to quickly solve optimally CHSP.

Our results show that the scaffolding-repeat problem formulations DRP and IRP seem to be sufficient to scaffold the repeats. This tends to validate our assumptions on the small number of repeats, and especially on the sufficiency of defining the repeats as pairs of equal (reversed) nucleotide sequences.

<footnote n="4">⁴https://pypi.org/project/khloraascaf/</footnote>
## Short-term future work

While our scaffolding problem formulation seems to be sufficient to retrieve the repeats, it seems it is not fully suitable for single-copies. If we have applied the maximum-weighted circuit problem to scaffold the single-copies after having scaffolded the repeats, it was only with the intention of staying in the context of global optimisation. On the one hand, having weights on links may have been more appropriate than just considering weights on the contigs: in some sense, that is the purpose of distances. On the other hand, *khloraascaf* could initially scaffold the repeats and then propose several solutions to link them, e.g. scored by the existence-weights.

Concerning the tests on synthetic data, we should use a more traditional assembly graph input: in fact, as revealed by comparing the *khloraascaf* results with the reference genomes, the used link generation suffers from local, or worse, global, misassemblies. As next step, we plan to inject *khloraascaf* into a state-of-the-art chloroplast genome pipeline, like *GetOrganelle*, and substitute what can be identified as the scaffolding part by our method. Hence, we should be able to relevantly compare *khloraascaf* approach with the state-of-the-art.

*khloraascaf* is sensible to the contig multiplicity computation. For now, a contig multiplicity is obtained by normalising its coverage by this of the starter. A better strategy may be to choose the smallest coverage for the normalisation as we expect the multiplicities to be upper-bounds of the contigs use.

Another issue concerns the choice of the starter: while it depends on the result of the mapping of *matK* gene map on the contigs, *DRP*, *IRP* and *SCP* problems may be adapted for a set of candidate for the starter.

To generalise *CHSP* on non-equally (reversed) pair of regions for the repeats, two combining ideas are proposed: on the one hand, we can add pairs of contigs to the repeated fragment sets from the user-input. On the other hand, *CHSP* should be able to handle the case when a single-copy region is in only one of the regions of a repeat: for now, the contiguity constraint and the objective exclude this case. The contiguity constraint can be adapted to accept only one contiguous region out of the two.

From a user-case perspective, the region graph data structure can be used to determine what genome forms are present in the read set, and in which proportion. Indeed, as the region graph explicitly describes the junctions between the regions, especially between the inverted repeats, one may extract the nucleotide sequences of these junctions to answer the existence of the forms in the read set.
Long-term future work

There are numerous generic methods for genome assembly and scaffolding in the literature. When applying generic methods to the genomes of specific organisms, the application often results in data subsampling strategies and hyper-parameterisation. We aimed to develop a scaffolding method that is aware of the specific characteristics of the data at the core of the problem formulation.

The challenge lies in the mathematical formalisation of biological knowledge. It is a matter of clearly postulating the knowledge behind the specific method. Postulating provides a solid basis for developing an approach to solve the problem, and facilitate discussion.

As we succeeded in formulating the scaffolding for chloroplast genomes, we would be able to propose a specific formulation for the other organelles in a plant cell. Mitochondria, for instance, are likely to have the same genomic structures as the chloroplasts. The next step could be to combine the dedicated formulations. If they each represent one objective, future work will involve multi-objective studies.

We can also adapt our scaffolding formulation for the repeats to identify the regions resulting from an event of repeat degeneration. In this case, we could apply our methods to genomes of different organisms, such as bacterial, plant or animal genomes, without the need to assemble them.
BIBLIOGRAPHY

Adey, A., Kitzman, J. O., Burton, J. N., Daza, R., Kumar, A., Christiansen, L., Ronaghi, M., Amini, S., Gunderson, K. L., Steemers, F. J., and Shendure, J. (2014). In vitro, long-range sequence information for de novo genome assembly via transposase contiguity. *Genome Research*, 24(12):2041–2049.

Aganezov, S., Avdeyev, P., Alexeev, N., Rong, Y., and Alekseyev, M. A. (2022). Orienting Ordered Scaffolds: Complexity and Algorithms. *SN Computer Science*, 3(4):308.

Andonov, R., Djidjev, H., François, S., and Lavenier, D. (2019). Complete assembly of circular and chloroplast genomes based on global optimization. *Journal of Bioinformatics and Computational Biology*, 17(3):1950014.

Ankenbrand, M. J., Pfaff, S., Terhoeven, N., Qureischi, M., Gündel, M., Weiß, C. L., Hackl, T., and Förster, F. (2018). chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data. *Journal of Open Source Software*, 3(21):464.

Bakker, F. T., Lei, D., Yu, J., Mohammadin, S., Wei, Z., van de Kerke, S., Gravendeel, B., Nieuwenhuis, M., Staats, M., Alquezar-Planas, D. E., and Holmer, R. (2016). Herbarium genomics: Plastome sequence assembly from a range of herbarium specimens using an Iterative Organelle Genome Assembly pipeline. *Biological Journal of the Linnean Society*, 117(1):33–43.

Bankevich, A., Nurk, S., Antipov, D., Gurevich, A. A., Dvorkin, M., Kulikov, A. S., Lesin, V. M., Nikolenko, S. I., Pham, S., Prjibelski, A. D., Pyshkin, A. V., Sirotkin, A. V., Vyahhi, N., Tesler, G., Alekseyev, M. A., and Pevzner, P. A. (2012). SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. *Journal of Computational Biology: A Journal of Computational Molecular Cell Biology*, 19(5):455–477.

Batzoglou, S., Jaffe, D. B., Stanley, K., Butler, J., Gnerre, S., Mauceli, E., Berger, B., Mesirov, J. P., and Lander, E. S. (2002). ARACHNE: A whole-genome shotgun assembler. *Genome Research*, 12(1):177–189.

Bendich, A. J. (1987). Why do chloroplasts and mitochondria contain so many copies of their genome? *BioEssays*, 6(6):279–282.

Bendich, A. J. (2004). Circular chloroplast chromosomes: The grand illusion. *The Plant Cell*, 16(7):1661–1666.
Bock, R. and Knoop, V., editors (2012). *Genomics of Chloroplasts and Mitochondria*, volume 35 of *Advances in Photosynthesis and Respiration*. Springer Netherlands, Dordrecht.

Bodily, P. M., Fujimoto, M. S., Snell, Q., Ventura, D., and Clement, M. J. (2016). ScaffoldScaffolder: Solving contig orientation via bidirected to directed graph reduction. *Bioinformatics*, 32(1):17–24.

Boetzer, M., Henkel, C. V., Jansen, H. J., Butler, D., and Pirovano, W. (2011). Scaffolding pre-assembled contigs using SSPACE. *Bioinformatics*, 27(4):578–579.

Boetzer, M. and Pirovano, W. (2012). Toward almost closed genomes with GapFiller. *Genome Biology*, 13(6):R56.

Briot, N., Chateau, A., Coletta, R., de Givry, S., Leleux, P., and Schiex, T. (2014). An integer linear programming approach for genome scaffolding. In *WCB: Workshop on Constraint-Based Methods for Bioinformatics*, 10th Workshop on Constraint-Based Methods for Bioinformatics (WCB), 2014, page 16 p., Lyon, France.

Burton, J. N., Adey, A., Patwardhan, R. P., Qiu, R., Kitzman, J. O., and Shendure, J. (2013). Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions. *Nature Biotechnology*, 31(12):1119–1125.

Chateau, A. and Giroudeau, R. (2014). Complexity and Polynomial-Time Approximation Algorithms around the Scaffolding Problem. In Dediu, A.-H., Martín-Vide, C., and Truthe, B., editors, *Algorithms for Computational Biology*, Lecture Notes in Computer Science, pages 47–58, Cham. Springer International Publishing.

Chateau, A. and Giroudeau, R. (2015). A complexity and approximation framework for the maximization scaffolding problem. *Theoretical Computer Science*, 595:92–106.

Cheng, H., Concepcion, G. T., Feng, X., Zhang, H., and Li, H. (2021). Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. *Nature Methods*, 18(2):170–175.

Chikhi, R. and Rizk, G. (2012). Space-Efficient and Exact de Bruijn Graph Representation Based on a Bloom Filter. In Raphael, B. and Tang, J., editors, *Algorithms in Bioinformatics*, Lecture Notes in Computer Science, pages 236–248, Berlin, Heidelberg. Springer.
Chin, C.-S., Peluso, P., Sedlazeck, F. J., Nattestad, M., Concepcion, G. T., Clum, A., Dunn, C., O'Malley, R., Figueroa-Balderas, R., Morales-Cruz, A., Cramer, G. R., Delledonne, M., Luo, C., Ecker, J. R., Cantu, D., Rank, D. R., and Schatz, M. C. (2016). Phased diploid genome assembly with single-molecule real-time sequencing. *Nature Methods*, 13(12):1050–1054.

Clark, S. C., Egan, R., Frazier, P. I., and Wang, Z. (2013). ALE: A generic assembly likelihood evaluation framework for assessing the accuracy of genome and metagenome assemblies. *Bioinformatics*, 29(4):435–443.

Coissac, E., Hollingsworth, P. M., Lavergne, S., and Taberlet, P. (2016). From barcodes to genomes: Extending the concept of DNA barcoding. *Molecular Ecology*, 25(7):1423–1428.

Coombe, L., Zhang, J., Vandervalk, B. P., Chu, J., Jackman, S. D., Birol, I., and Warren, R. L. (2018). ARKS: Chromosome-scale scaffolding of human genome drafts with linked read kmers. *BMC Bioinformatics*, 19(1):234.

Daniell, H., Lin, C.-S., Yu, M., and Chang, W.-J. (2016). Chloroplast genomes: Diversity, evolution, and applications in genetic engineering. *Genome Biology*, 17(1):134.

Davot, T., Chateau, A., Fossé, R., Giroudeau, R., and Weller, M. (2022). On a greedy approach for genome scaffolding. *Algorithms for Molecular Biology*, 17(1):16.

Dayarian, A., Michael, T. P., and Sengupta, A. M. (2010). SOPRA: Scaffolding algorithm for paired reads via statistical optimization. *BMC Bioinformatics*, 11(1):345.

Delahaye, C. and Nicolas, J. (2021). Sequencing DNA with nanopores: Troubles and biases. *PLOS ONE*, 16(10):e0257521.

Deng, X.-W., Wing, R. A., and Gruissem, W. (1989). The chloroplast genome exists in multimeric forms. *Proceedings of the National Academy of Sciences*, 86(11):4156–4160.

Dierckxsens, N., Mardulyn, P., and Smits, G. (2017). NOVOPlasty: De novo assembly of organelle genomes from whole genome data. *Nucleic Acids Research*, 45(4):e18.

Donmez, N. and Brudno, M. (2013). SCARPA: Scaffolding reads with practical algorithms. *Bioinformatics*, 29(4):428–434.
Edmonds, J. and Johnson, E. L. (1970). Matching: A well-solved class of integer linear programs. In *In: Combinatorial Structures and Their Applications (Gordon and Breach)*, pages 89–92.

Ekim, B., Berger, B., and Chikhi, R. (2021). Minimizer-space de Bruijn graphs: Whole-genome assembly of long reads in minutes on a personal computer. *Cell Systems*, 12(10):958–968.e6.

François, S., Andonov, R., Lavenier, D., and Djidjev, H. (2018). Global Optimization for Scaffolding and Completing Genome Assemblies. *Electronic Notes in Discrete Mathematics*, 64:185–194.

Franklin, R. E. and Gosling, R. G. (1953). The structure of sodium thymonucleate fibres. I. The influence of water content. *Acta Crystallographica*, 6(8):673–677.

Freudenthal, J. A., Pfaff, S., Terhoeven, N., Korte, A., Ankenbrand, M. J., and Förster, F. (2020). A systematic comparison of chloroplast genome assembly tools. *Genome Biology*, 21:254.

Gao, S., Nagarajan, N., and Sung, W.-K. (2011). Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences. In Bafna, V. and Sahinalp, S. C., editors, *Research in Computational Molecular Biology*, Lecture Notes in Computer Science, pages 437–451, Berlin, Heidelberg. Springer.

Ghurye, J., Rhie, A., Walenz, B. P., Schmitt, A., Selvaraj, S., Pop, M., Phillippy, A. M., and Koren, S. (2019). Integrating Hi-C links with assembly graphs for chromosome-scale assembly. *PLOS Computational Biology*, 15(8):e1007273.

Gingeras, T., Milazzo, J., Sciaky, D., and Roberts, R. (1979). Computer programs for the assembly of DNA sequences. *Nucleic Acids Research*, 7(2):529–543.

Goldfarb, D. and Idnani, A. (1983). A numerically stable dual method for solving strictly convex quadratic programs. *Mathematical Programming*, 27(1):1–33.

Gould, S. B. (2012). Algae's complex origins. *Nature*, 492(7427):46–48.

Gritsenko, A. A., Nijkamp, J. F., Reinders, M. J., and de Ridder, D. (2012). GRASS: A generic algorithm for scaffolding next-generation sequencing assemblies. *Bioinformatics*, 28(11):1429–1437.

Gurevich, A., Saveliev, V., Vyahhi, N., and Tesler, G. (2013). QUAST: Quality assessment tool for genome assemblies. *Bioinformatics*, 29(8):1072–1075.
Gusfield, D. (2019). The RNA-Folding Problem. In Gusfield, D., editor, *Integer Linear Programming in Computational and Systems Biology: An Entry-Level Text and Course*, pages 105–121. Cambridge University Press, Cambridge.

Hernandez, D., François, P., Farinelli, L., Østerås, M., and Schrenzel, J. (2008). De novo bacterial genome sequencing: Millions of very short reads assembled on a desktop computer. *Genome Research*, 18(5):802–809.

Hierholzer, C. and Wiener, C. (1873). Ueber die Möglichkeit, einen Linienzug ohne Wiederholung und ohne Unterbrechung zu umfahren. *Mathematische Annalen*, 6(1):30–32.

Hon, T., Mars, K., Young, G., Tsai, Y.-C., Karalius, J. W., Landolin, J. M., Maurer, N., Kudrna, D., Hardigan, M. A., Steiner, C. C., Knapp, S. J., Ware, D., Shapiro, B., Peluso, P., and Rank, D. R. (2020). Highly accurate long-read HiFi sequencing data for five complex genomes. *Scientific Data*, 7(1):399.

Huson, D. H., Reinert, K., and Myers, E. W. (2002). The greedy path-merging algorithm for contig scaffolding. *Journal of the ACM*, 49(5):603–615.

Idury, R. M. and Waterman, M. S. (1995). A New Algorithm for DNA Sequence Assembly. *Journal of Computational Biology*, 2(2):291–306.

Jain, C. (2023). Coverage-preserving sparsification of overlap graphs for long-read assembly. *Bioinformatics*, 39(3):btad124.

Jin, J.-J., Yu, W.-B., Yang, J.-B., Song, Y., dePamphilis, C. W., Yi, T.-S., and Li, D.-Z. (2020). GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes. *Genome Biology*, 21(1):241.

Kamath, G. M., Shomorony, I., Xia, F., Courtade, T. A., and Tse, D. N. (2017). HINGE: Long-read assembly achieves optimal repeat resolution. *Genome Research*, 27(5):747–756.

Karp, R. M. (1972). Reducibility among Combinatorial Problems. In Miller, R. E., Thatcher, J. W., and Bohlinger, J. D., editors, *Complexity of Computer Computations: Proceedings of a Symposium on the Complexity of Computer Computations, Held March 20–22, 1972, at the IBM Thomas J. Watson Research Center, Yorktown Heights, New York, and Sponsored by the Office of Naval Research, Mathematics Program, IBM World Trade Corporation, and the IBM Research Mathematical Sciences Department*, The IBM Research Symposia Series, pages 85–103. Springer US, Boston, MA.
Kececioglu, J. D. (1991). *Exact and Approximation Algorithms for DNA Sequence Reconstruction*. PhD thesis, The University of Arizona.

Kececioglu, J. D. and Myers, E. W. (1995). Combinatorial algorithms for DNA sequence assembly. *Algorithmica*, 13(1):7.

Kidwell, M. G. (2002). Transposable elements and the evolution of genome size in eukaryotes. *Genetica*, 115(1):49-63.

Kim, K.-J. and Lee, H.-L. (2005). Widespread Occurrence of Small Inversions in the Chloroplast Genomes of Land Plants. *Molecules and Cells*, 19(1):104-113.

Kobayashi, T., Takahara, M., Miyagishima, S.-y., Kuroiwa, H., Sasaki, N., Ohta, N., Matsuzaki, M., and Kuroiwa, T. (2002). Detection and Localization of a Chloroplast--encoded HU-Like Protein That Organizes Chloroplast Nucleoids. *The Plant Cell*, 14(7):1579-1589.

Kolmogorov, M., Yuan, J., Lin, Y., and Pevzner, P. A. (2019). Assembly of long, error-prone reads using repeat graphs. *Nature Biotechnology*, 37(5):540-546.

Korbel, J. O., Urban, A. E., Affournit, J. P., Godwin, B., Grubert, F., Simons, J. F., Kim, P. M., Palejev, D., Carriero, N. J., Du, L., Taillon, B. E., Chen, Z., Tanzer, A., Saunders, A. C. E., Chi, J., Yang, F., Carter, N. P., Hurles, M. E., Weissman, S. M., Harkins, T. T., Gerstein, M. B., Egholm, M., and Snyder, M. (2007). Paired-End Mapping Reveals Extensive Structural Variation in the Human Genome. *Science*, 318(5849):420-426.

Koren, S., Treangen, T. J., and Pop, M. (2011). Bambus 2: Scaffolding metagenomes. *Bioinformatics*, 27(21):2964-2971.

Koren, S., Walenz, B. P., Berlin, K., Miller, J. R., Bergman, N. H., and Phillippy, A. M. (2017). Canu: Scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. *Genome Research*, 27(5):722-736.

Kumar, R. A., Oldenburg, D. J., and Bendich, A. J. (2014). Changes in DNA damage, molecular integrity, and copy number for plastid DNA and mitochondrial DNA during maize development. *Journal of Experimental Botany*, 65(22):6425-6439.

Li, H. (2016). Minimap and miniasm: Fast mapping and de novo assembly for noisy long sequences. *Bioinformatics*, 32(14):2103-2110.

Li, M. (1990). Towards a DNA sequencing theory (learning a string). In *Proceedings [1990] 31st Annual Symposium on Foundations of Computer Science*, pages 125-134 vol.1.
Li, R., Zhu, H., Ruan, J., Qian, W., Fang, X., Shi, Z., Li, Y., Li, S., Shan, G., Kristiansen, K., Li, S., Yang, H., Wang, J., and Wang, J. (2010). De novo assembly of human genomes with massively parallel short read sequencing. *Genome Research*, 20(2):265-272.

Lokshtanov, D., Saurabh, S., and Sikdar, S. (2009). Simpler Parameterized Algorithm for OCT. In Fiala, J., Kratochvíl, J., and Miller, M., editors, *Combinatorial Algorithms*, Lecture Notes in Computer Science, pages 380–384, Berlin, Heidelberg. Springer.

Long, L., Li, Y., Wang, S., Liu, Z., Wang, J., and Yang, M. (2023). Complete chloroplast genomes and comparative analysis of Ligustrum species. *Scientific Reports*, 13(1):212.

Luo, J., Lyu, M., Chen, R., Zhang, X., Luo, H., and Yan, C. (2019). SLR: A scaffolding algorithm based on long reads and contig classification. *BMC Bioinformatics*, 20(1):1–11.

Luo, J., Wang, J., Zhang, Z., Li, M., and Wu, F.-X. (2017). BOSS: A novel scaffolding algorithm based on an optimized scaffold graph. *Bioinformatics*, 33(2):169–176.

Luo, R., Liu, B., Xie, Y., Li, Z., Huang, W., Yuan, J., He, G., Chen, Y., Pan, Q., Liu, Y., Tang, J., Wu, G., Zhang, H., Shi, Y., Liu, Y., Yu, C., Wang, B., Lu, Y., Han, C., Cheung, D. W., Yiu, S.-M., Peng, S., Xiaoqian, Z., Liu, G., Liao, X., Li, Y., Yang, H., Wang, J., Lam, T.-W., and Wang, J. (2012). SOAPdenovo2: An empirically improved memory-efficient short-read de novo assembler. *GigaScience*, 1(1):18.

Mandric, I. and Zelikovsky, A. (2015). ScaffMatch: Scaffolding algorithm based on maximum weight matching. *Bioinformatics*, 31(16):2632-2638.

Marçais, G. and Kingsford, C. (2011). A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. *Bioinformatics*, 27(6):764-770.

McKain, M. and afinit (2017). Mrmckain/Fast-Plast: Fast-Plast v.1.2.6. Zenodo.

Medvedev, P. and Brudno, M. (2009). Maximum Likelihood Genome Assembly. *Journal of Computational Biology*, 16(8):1101–1116.

Medvedev, P., Georgiou, K., Myers, G., and Brudno, M. (2007). Computability of Models for Sequence Assembly. In Giancarlo, R. and Hannenhalli, S., editors, *Algorithms in Bioinformatics*, Lecture Notes in Computer Science, pages 289–301, Berlin, Heidelberg. Springer.
Mehdi, K., Gibrat, J.-F., and Elloumi, M. (2017). Generations of Sequencing Technologies: From First to Next Generation. *Electromagnetic Biology and Medicine*, 9(3):8 p.

Miller, C. E., Tucker, A. W., and Zemlin, R. A. (1960). Integer Programming Formulation of Traveling Salesman Problems. *Journal of the ACM*, 7(4):326-329.

Miller, J. R., Delcher, A. L., Koren, S., Venter, E., Walenz, B. P., Brownley, A., Johnson, J., Li, K., Mobarry, C., and Sutton, G. (2008). Aggressive assembly of pyrosequencing reads with mates. *Bioinformatics* (Oxford, England), 24(24):2818-2824.

Myers, E. W. (1995). Toward simplifying and accurately formulating fragment assembly. *Journal of Computational Biology: A Journal of Computational Molecular Cell Biology*, 2(2):275-290.

Myers, E. W. (2005). The fragment assembly string graph. *Bioinformatics*, 21(suppl_2):ii79-ii85.

Myers, E. W., Sutton, G. G., Delcher, A. L., Dew, I. M., Fasulo, D. P., Flanigan, M. J., Kravitz, S. A., Mobarry, C. M., Reinert, K. H., Remington, K. A., Anson, E. L., Bolanos, R. A., Chou, H. H., Jordan, C. M., Halpern, A. L., Lonardi, S., Beasley, E. M., Brandon, R. C., Chen, L., Dunn, P. J., Lai, Z., Liang, Y., Nusskern, D. R., Zhan, M., Zhang, Q., Zheng, X., Rubin, G. M., Adams, M. D., and Venter, J. C. (2000). A whole-genome assembly of Drosophila. *Science* (New York, N.Y.), 287(5461):2196-2204.

O'Connor, M., Peifer, M., and Bender, W. (1989). Construction of Large DNA Segments in Escherichia coli. *Science*, 244(4910):1307-1312.

Ohyama, K., Fukuzawa, H., Kohchi, T., Shirai, H., Sano, T., Sano, S., Umesono, K., Shiki, Y., Takeuchi, M., Chang, Z., Aota, S.-i., Inokuchi, H., and Ozeki, H. (1986). Chloroplast gene organization deduced from complete sequence of liverwort Marchantia polymorpha chloroplast DNA. *Nature*, 322(6079):572-574.

Palmer, J. D. (1983). Chloroplast DNA exists in two orientations. *Nature*, 301(5895):92-93.

Palmer, J. D. (1985). Comparative Organization of Chloroplast Genomes. *Annual Review of Genetics*, 19(1):325-354.

Peltola, H., Söderlund, H., and Ukkonen, E. (1984). SEQAID: A DNA sequence assembling program based on a mathematical model. *Nucleic Acids Research*, 12(1 Pt 1):307-321.
Pervez, M. T., ul Hasnain, M. J., Abbas, S. H., Moustafa, M. F., Aslam, N., and Shah, S. S. M. (2022). A Comprehensive Review of Performance of Next-Generation Sequencing Platforms. *BioMed Research International*, 2022:e3457806.

Pevzner, P. A. (1989). L-Tuple DNA Sequencing: Computer Analysis. *Journal of Biomolecular Structure and Dynamics*, 7(1):63-73.

Pevzner, P. A., Tang, H., and Waterman, M. S. (2001). An Eulerian path approach to DNA fragment assembly. *Proceedings of the National Academy of Sciences of the United States of America*, 98(17):9748-9753.

Poloczek, M. and Szegedy, M. (2012). Randomized Greedy Algorithms for the Maximum Matching Problem with New Analysis. In *2012 IEEE 53rd Annual Symposium on Foundations of Computer Science*, pages 708-717.

Pop, M., Kosack, D. S., and Salzberg, S. L. (2004). Hierarchical Scaffolding With Bambus. *Genome Research*, 14(1):149-159.

Qin, M., Wu, S., Li, A., Zhao, F., Feng, H., Ding, L., and Ruan, J. (2019). LRScaf: Improving draft genomes using long noisy reads. *BMC Genomics*, 20(1):955.

Räihä, K.-J. and Ukkonen, E. (1981). The Shortest Common Supersequence Problem over Binary Alphabet is NP-Complete. *Theor. Comput. Sci.*

Rice, E. S. and Green, R. E. (2019). New Approaches for Genome Assembly and Scaffolding. *Annual Review of Animal Biosciences*, 7:17-40.

Roy, R. S., Chen, K. C., Sengupta, A. M., and Schliep, A. (2012). SLIQ: Simple Linear Inequalities for Efficient Contig Scaffolding. *Journal of Computational Biology*, 19(10):1162-1175.

Ruan, J. and Li, H. (2020). Fast and accurate long-read assembly with wtdbg2. *Nature Methods*, 17(2):155-158.

Sahlin, K., Vezzi, F., Nystedt, B., Lundeberg, J., and Arvestad, L. (2014). BESST - Efficient scaffolding of large fragmented assemblies. *BMC Bioinformatics*, 15(1):281.

Salmela, L., Mäkinen, V., Välimäki, N., Ylinen, J., and Ukkonen, E. (2011). Fast scaffolding with small independent mixed integer programs. *Bioinformatics*, 27(23):3259-3265.
Sancho, R., Cantalapiedra, C. P., López-Alvarez, D., Gordon, S. P., Vogel, J. P., Catalán, P., and Contreras-Moreira, B. (2018). Comparative plastome genomics and phylogenomics of Brachypodium: Flowering time signatures, introgression and recombination in recently diverged ecotypes. *New Phytologist*, 218(4):1631–1644.

Sanger, F., Coulson, A. R., Hong, G. F., Hill, D. F., and Petersen, G. B. (1982). Nucleotide sequence of bacteriophage λ DNA. *Journal of Molecular Biology*, 162(4):729–773.

Sanger, F., Nicklen, S., and Coulson, A. R. (1977). DNA sequencing with chain-terminating inhibitors. *Proceedings of the National Academy of Sciences*, 74(12):5463–5467.

Schrijver, A. (2003). *Combinatorial Optimization: Polyhedra and Efficiency*. Springer Science & Business Media.

Seyer, P., Kowallik, K. V., and Herrmann, R. G. (1981). A physical map of Nicotiana tabacum plastid DNA including the location of structural genes for ribosomal RNAs and the large subunit of ribulose bisphosphate carboxylase/oxygenase. *Current Genetics*, 3(3):189–204.

Shafin, K., Pesout, T., Lorig-Roach, R., Haukness, M., Olsen, H. E., Bosworth, C., Armstrong, J., Tigyi, K., Maurer, N., Koren, S., Sedlazeck, F. J., Marschall, T., Mayes, S., Costa, V., Zook, J. M., Liu, K. J., Kilburn, D., Sorensen, M., Munson, K. M., Vollger, M. R., Monlong, J., Garrison, E., Eichler, E. E., Salama, S., Haussler, D., Green, R. E., Akeson, M., Phillippy, A., Miga, K. H., Carnevali, P., Jain, M., and Paten, B. (2020). Nanopore sequencing and the Shasta toolkit enable efficient de novo assembly of eleven human genomes. *Nature Biotechnology*, 38(9):1044–1053.

Shapiro, M. B. (1967). An Algorithm for Reconstructing Protein and RNA Sequences. *Journal of the ACM*, 14(4):720–731.

Shinozaki, K., Ohme, M., Tanaka, M., Wakasugi, T., Hayashida, N., Matsubayashi, T., Zaita, N., Chunwongse, J., Obokata, J., Yamaguchi-Shinozaki, K., Ohto, C., Torazawa, K., Meng, B., Sugita, M., Deno, H., Kamogashira, T., Yamada, K., Kusuda, J., Takaiwa, F., Kato, A., Tohdoh, N., Shimada, H., and Sugiura, M. (1986). The complete nucleotide sequence of the tobacco chloroplast genome: Its gene organization and expression. *The EMBO Journal*, 5(9):2043–2049.

Sommer, D. D., Delcher, A. L., Salzberg, S. L., and Pop, M. (2007). Minimus: A fast, lightweight genome assembler. *BMC bioinformatics*, 8:64.
Staden, R. (1979). A strategy of DNA sequencing employing computer programs. *Nucleic Acids Research*, 6(7):2601-2610.

Staden, R. (1980). A new computer method for the storage and manipulation of DNA gel reading data. *Nucleic Acids Research*, 8(16):3673-3694.

Sun, J., Wang, Y., Liu, Y., Xu, C., Yuan, Q., Guo, L., and Huang, L. (2020). Evolutionary and phylogenetic aspects of the chloroplast genome of Chaenomeles species. *Scientific Reports*, 10(1):11466.

Thode, V. A., Oliveira, C. T., Loeuille, B., Siniscalchi, C. M., and Pirani, J. R. (2021). Comparative analyses of Mikania (Asteraceae: Eupatorieae) plastomes and impact of data partitioning and inference methods on phylogenetic relationships. *Scientific Reports*, 11(1):13267.

Trevors, J. T. (1996). Genome size in bacteria. *Antonie van Leeuwenhoek*, 69(4):293-303.

Tsai, C.-H. and Strauss, S. H. (1989). Dispersed repetitive sequences in the chloroplast genome of Douglas-fir. *Current Genetics*, 16(3):211-218.

Turmel, M., Otis, C., and Lemieux, C. (2017). Divergent copies of the large inverted repeat in the chloroplast genomes of ulvophycean green algae. *Scientific Reports*, 7(1):994.

Venter, J. C., Smith, H. O., and Hood, L. (1996). A new strategy for genome sequencing. *Nature*, 381(6581):364-366.

Wang, W. and Lanfear, R. (2019). Long-Reads Reveal That the Chloroplast Genome Exists in Two Distinct Versions in Most Plants. *Genome Biology and Evolution*, 11(12):3372-3381.

Watson, J. D. and Crick, F. H. C. (1953). The Structure of Dna. *Cold Spring Harbor Symposia on Quantitative Biology*, 18:123-131.

Weller, M., Chateau, A., and Giroudeau, R. (2015). Exact approaches for scaffolding. *BMC Bioinformatics*, 16(14):S2.

Xiao-Ming, Z., Junrui, W., Li, F., Sha, L., Hongbo, P., Lan, Q., Jing, L., Yan, S., Weihua, Q., Lifang, Z., Yunlian, C., and Qingwen, Y. (2017). Inferring the evolutionary mechanism of the chloroplast genome size by comparing whole-chloroplast genome sequences in seed plants. *Scientific Reports*, 7(1):1555.

Zerbino, D. R. and Birney, E. (2008). Velvet: Algorithms for de novo short read assembly using de Bruijn graphs. *Genome Research*, 18(5):821-829.
<figure><img src="image_96.png" /></figure>
# APPENDIX

## In this chapter

<table><tr><td>1</td><td>Repeated fragment set functions</td><td>A1</td></tr><tr><td>2</td><td>Reduction of the repeated fragment sets</td><td>A2</td></tr><tr><td>2.1</td><td>Repeated fragment set reductions.</td><td>A2</td></tr><tr><td>2.2</td><td>Pairs of repeated fragment set reductions.</td><td>A5</td></tr><tr><td>2.3</td><td>Adjacent repeated fragment set reductions.</td><td>A8</td></tr><tr><td>3</td><td>Metrics</td><td>A17</td></tr><tr><td>3.1</td><td>Quast metrics.</td><td>A17</td></tr><tr><td>4</td><td>Supplementary results</td><td>A18</td></tr><tr><td>4.1</td><td>v1 scaffolding benchmark.</td><td>A18</td></tr><tr><td>4.2</td><td>v2 scaffolding benchmark.</td><td>A20</td></tr></table>

## 1 Repeated fragment set functions

dirfrag: $V \to DirF$

$$v \mapsto \begin{cases} (v, v') & \text{if } occ_v | 2 \\ & \text{where } occ_{v'} = occ_v + 1 \\ (v', v) & \text{else } occ_v \nmid 2 \\ & \text{where } occ_{v'} = occ_v - 1 \\ or_v & = or_{v'} \end{cases}$$

invfrag: $V \to InvF$

$$v \mapsto \begin{cases} (v, v') & \text{if } or_v = 0 \\ & \text{where } or_{v'} = 1 \land occ_{v'} = occ_v + 1 \\ (v', v) & \text{else } or_v = 1 \\ & \text{where } or_{v'} = 0 \land occ_{v'} = occ_v - 1 \end{cases}$$
$diradj: E \to E$

$(u, v) \mapsto (u', v')$

$ctg_{u'} = ctg_{u}$

$or_{u'} = or_{u}$

$occ_{u'} = \begin{cases} occ_{u} + 1 & \text{if } occ_{u} | 2 \\ occ_{u} - 1 & \text{otherwise} \end{cases}$

$ctg_{v'} = ctg_{v}$

$or_{v'} = or_{v}$

$occ_{v'} = \begin{cases} occ_{v} + 1 & \text{if } occ_{v} | 2 \\ occ_{v} - 1 & \text{otherwise} \end{cases}$

$invadj: E \to E$

$(u, v) \mapsto (v', u')$

$ctg_{v'} = ctg_{v}$

$or_{v'} = 1 - or_{v}$

$occ_{v'} = 2 \lfloor \frac{occ_{v}}{2} \rfloor + (1 - or_{v})$

$ctg_{u'} = ctg_{u}$

$or_{u'} = 1 - or_{u}$

$occ_{u'} = 2 \lfloor \frac{occ_{u}}{2} \rfloor + (1 - or_{u})$

# 2 Reduction of the repeated fragment sets

In this section, we give prove the minimality and the sufficiency of the repeated fragment sets (*RepF*, *PRepF*, *ARepF*) definitions.

## 2.1 Repeated fragment set reductions

Here we consider the reduction operations for *DirF* and *InvF*, and conclude on their minimality. Only two vertices $i, j \in V$ with the same identifier, i.e. $ctg_i = ctg_j$ can form a repeated fragment. The relative orientations between *i* and *j* is determined according to the type of the repeat the fragment is associated to. The occurrences of the two vertices in a repeated fragment must differ.
Two combinatorial reductions are necessary — commutative and pairing reductions:

i. commutative reduction means that $\forall (i, j) \in RepF, (j, i) \notin RepF$;

ii. pairing reduction means that it is not necessary to pair all the occurrences cross all the occurrences, and we can group by consecutive occurrences without any intersection.

Reminder of the *DirF* definition:

$$DirF = \bigcup_{c \in \mathcal{R}} \left\{ \begin{array}{l} (i, j) \in V^2 \text{ s.t.} \\ \text{ctg}_i = \text{ctg}_j = c \\ \wedge \text{or}_i = \text{or}_j \in \{f, r\} \\ \wedge \text{occ}_i = \text{occ}_j - 1 = 2k \\ 0 \le k < \left\lfloor \frac{\text{mult}(c)}{2} \right\rfloor \end{array} \right\}$$

**► Proposition 2.1: Commutative reduction for DirF**

<figure><img src="image_97.png" /></figure>

$$\forall (i, j) \in DirF, (j, i) \notin DirF$$

**▷ Proof**

$$\text{Let } (i, j) \in \text{DirF. } \text{occ}_j = \text{occ}_i + 1 \text{ so } (j, i) \notin \text{DirF}$$

<figure><img src="image_98.png" /></figure>

**► Proposition 2.2: Pairing reduction for DirF**

<figure><img src="image_99.png" /></figure>

The direct fragment set contains a sufficient number of direct fragments to form all the solutions.

**▷ Proof**

The idea is to find the minimum occurrence difference between two direct fragments such that they can both be chosen to form DRs. Let $(i, j) \in DirF$ be a direct fragment and suppose that it is chosen to participate into a DR.

Let $k, l \in V$ | $ctg_k = ctg_l = ctg_i = ctg_j \wedge or_k = or_l = or_i = or_j$ be a candidate for the next direct fragment:

- the occurrences of $k$ and $l$ must differ from the ones of $i$ and $j$, as a vertex can occupy at most one position, so $occ_j < occ_k$;
— according to Proposition 2.1, $occ_k < occ_l$ and since occurrences are integers, $occ_k = occ_l - 1$.

So $occ_i = occ_j - 1 \le occ_k - 2 = occ_l - 3 \implies occ_i + 3 \le occ_l$. The smallest possible values for $occ_k$ and $occ_l$ equal $occ_i + 2$ and $occ_j + 2$.

⊲

### ▶ Proposition 2.3

Direct fragment set $DirF$ is a minimal set for $DRP$.

### ▷ Proof

Firstly, we must keep the direct fragment for both the forward ($or_i = or_j = 0$) and the reverse orientations ($or_i = or_j = 1$) as the scaffolding problem aims to order and orientate the contigs, so by definition we cannot determine which orientation will participate in the order. Secondly, making step bigger than 2k prevents to use all the occurrences of a contig, so it contradicts the contig multiplicity definition.

⊲

Reminder of the $InvF$ definition:

$$InvF = \bigcup_{c \in \mathbb{R}} \left\{ \begin{array}{l} (i, j) \in V^2 \text{ s.t.} \\ \text{ctg}_i = \text{ctg}_j = c \\ \wedge or_i = f \wedge or_j = r \\ \wedge occ_i = occ_j - 1 = 2k \\ 0 \le k < \left\lfloor \frac{\text{mult}(c)}{2} \right\rfloor \end{array} \right\}$$

### ▶ Proposition 2.4: Commutative reduction for $InvF$

<figure><img src="image_100.png" /></figure>

$$\forall (i, j) \in InvF, (j, i) \notin InvF$$

### ▷ Proof

Let $(i, j) \in InvF$. $occ_j = occ_i + 1$ so $(j, i) \notin InvF$

⊲
► **Proposition 2.5: Pairing reduction for $InvF$**

The inverted fragment set contains a sufficient number of inverted fragments to form all the solutions.

▷ Proof

The idea is to find the minimum occurrence difference between two inverted fragments such that they can both be chosen to form IRs. Let $(i, j) \in InvF$ be an inverted fragment and suppose that it is chosen to participate into an IR. Let $k, l \in V$ $|$ $ctg_k = ctg_l = ctg_i = ctg_j \wedge or_k \neq or_l$ be a candidate for the next direct fragment:

- to respect the property given by Proposition 2.4, $or_k = 1 - or_l = 0$

- the occurrences of $k$ and $l$ must differ from the ones of $i$ and $j$ when their orientations are matching, as a vertex can occupy at most one position. Furthermore, for a given occurrence and a given identifier, the two orientations are mutually exclusive, so $occ_j < occ_k$;

- according to Proposition 2.4, $occ_k < occ_l$ and since occurrences are integers, $occ_k = occ_l - 1$.

So $occ_i = occ_j - 1 \leq occ_k - 2 = occ_l - 3 \implies occ_i + 3 \leq occ_l$. The smallest possible values for $occ_k$ and $occ_l$ equal $occ_i + 2$ and $occ_j + 2$.

<figure><img src="image_101.png" /></figure>

► **Proposition 2.6**

Inverted fragment set $InvF$ is a minimal set for $IRP$.

▷ Proof

Step bigger than 2k prevents to use all the occurrences of a contig, so it contradicts the contig multiplicity definition.

<figure><img src="image_102.png" /></figure>

## 2.2 Pairs of repeated fragment set reductions

Here we consider the reduction operations for $PDirF$ and $PInvF$, and conclude on their minimality. Just the *commutative reduction* is necessary.
Reminder of pairs of direct fragments set definition:

$$ PDirF = \left\{ \begin{array}{l} ((i, j), (k, l)) \in DirF^2 \text{ s.t.} \\ ctg_j < ctg_k \\ \lor \\ ctg_j = ctg_k \land occ_j < occ_k \end{array} \right\} $$

### ▶ Proposition 2.7: Commutative reduction for PDirF

$$ \forall ((i, j), (k, l)) \in PDirF, ((k, l), (i, j)) \notin PDirF $$

### ▷ Proof

Let $((i, j), (k, l)) \in PDirF$:

- if $ctg_j < ctg_k$ thus $ctg_l > ctg_i$ so $((k, l), (i, j)) \notin PDirF$

- else $ctg_j = ctg_k \land occ_j < occ_k$ thus $occ_l > occ_i$ so $((k, l), (i, j)) \notin PDirF$

△

### ▶ Proposition 2.8

$$ \text{Pair of direct fragment set } PDirF \text{ is a minimal set for } DRP. $$

### ▷ Proof

By reductio ad absurdum: if PDirF can be minimised to a PDirF' set, then there exists a pair of direct fragments $(p, q) \in PDirF | (p, q) \notin PDirF'$.

Let $p = (i, j)$ and $q = (k, l)$, and let $q_2 = (m, n) | (p, q_2) \in PDirF' \land (q, q_2) \in PDirF'$. i, j, k, l, m, n vertices can simultaneously participate in the solution.

Let build a counter-example: suppose that we have the order i k l j m n. This implies the direct fragments p and q fall into a forbidden case for DRP, but this is not detected by Constraints C14 and C15 because $(p, q) \notin PDirF'$. According to Constraint C16, the direct fragments p and $q_2$ can both match, and direct fragments q and $q_2$ too. So $m_p = m_q = m_{q_2} = 1$: absurd because p and q are nested so $m_p$ and $m_q$ should be equal to 0.

△
Reminder of pairs of inverted fragments set definition:

$$ PInvF = \left\{ \begin{array}{l} ((i, j), (k, l)) \in InvF^2 \text{ s.t.} \\ \text{ctg}_j < \text{ctg}_k \\ \lor \\ \text{ctg}_j = \text{ctg}_k \land \text{occ}_j < \text{occ}_k \end{array} \right\} $$

► **Proposition 2.9: Commutative reduction for PInvF**

$$ \forall ((i, j), (k, l)) \in PInvF, ((k, l), (i, j)) \notin PInvF $$

▷ Proof

Let $((i, j), (k, l)) \in PInvF$:

- if $ctg_j < ctg_k$ thus $ctg_l > ctg_i$ so $((k, l), (i, j)) \notin PInvF$

- else $ctg_j = ctg_k \land occ_j < occ_k$ thus $occ_l > occ_i$ so $((k, l), (i, j)) \notin PInvF$

▷

► **Proposition 2.10**

$$ \text{Pair of inverted fragment set } PInvF \text{ is a minimal set for } \mathcal{IRP}. $$

▷ Proof

By reductio ad absurdum: if $PInvF$ can be minimised to a $PInvF'$ set, then there exists a pair of inverted fragments $(p, q) \in PInvF | (p, q) \notin PInvF'$.

Let be $p = (i, j)$ and $q = (k, l)$, and let $q_2 = (m, n) | (p, q_2) \in PInvF' \land (q, q_2) \in PInvF'$. $i, j, k, l, m, n$ vertices can simultaneously participate in the solution.

Let build a counter-example: suppose that we have the order $ikjlmn$. This implies the inverted fragments $p$ and $q$ fall into a forbidden case for $\mathcal{DRP}$, but this is not detected by Constraints C14 and C15 because $(p, q) \notin PInvF'$. According to Constraint C16, the inverted fragments $p$ and $q_2$ can both match, and inverted fragments $q$ and $q_2$ too. So $m_p = m_q = m_{q_2} = 1$: absurd because $p$ and $q$ intersect so $m_p$ and $m_q$ should be equal to 0.

▷
## 2.3 Adjacent repeated fragment set reductions

Here we consider the reduction operations for $ADirF$ and $AInvF$, and conclude on their minimality.

Reminder of adjacent direct fragments set definition:

$$ADirF = \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ \quad ctg_u \neq ctg_v \\ \quad \wedge occ_u = 2k \\ \quad 0 \le k < \left\lfloor \frac{mult_u}{2} \right\rfloor \\ \quad \wedge occ_v = 2k' \\ \quad 0 \le k' < \left\lfloor \frac{mult_v}{2} \right\rfloor \end{array} \right\} \cup \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ \quad ctg_u = ctg_v \\ \quad \wedge (or_u = 0 \lor or_v = 0) \\ \quad \wedge occ_u = 2k \\ \quad \wedge occ_v = 2k' \\ \quad 0 \le k < k' < \left\lfloor \frac{mult_u}{2} \right\rfloor \end{array} \right\}$$

This set is the minimum exhaustive set of canonical edges which represent adjacent direct fragments. Two combinatorial reductions are necessary — *canonical reduction for different identifiers* (Propositions 2.11 and 2.12) and *same identifiers occurrences* reductions (Propositions 2.13 and 2.14).

### ► Proposition 2.11: Canonical reductions for $ADirF$ for different identifiers

$ADirF$ set contains only two of the eight edges that exist between two adjacent direct fragments when the vertex' identifiers are different, i.e.
$$ \forall (i,j) \in \text{DirF}, \forall (k,l) \in \text{DirF} \mid \text{ctg}_i \neq \text{ctg}_k: $$

$$ \begin{array}{lcl} (i,k) & \in & \text{ADirF} \\ \wedge(\bar{k}, \bar{i}) & \in & \text{ADirF} \\ \wedge(i,l) & \notin & \text{ADirF} \\ (i,k) \in E & \iff & \wedge(\bar{l}, \bar{i}) \notin \text{ADirF} \\ \iff (i,l) \in E & \iff & \wedge(j,k) \notin \text{ADirF} \\ \iff (j,k) \in E & \iff & \wedge(\bar{k}, \bar{j}) \notin \text{ADirF} \\ \iff (j,l) \in E & \iff & \wedge(j,l) \notin \text{ADirF} \\ \wedge(\bar{l}, \bar{j}) & \notin & \text{ADirF} \\ (k,i) & \in & \text{ADirF} \\ \wedge(\bar{i}, \bar{k}) & \in & \text{ADirF} \\ \wedge(k,j) & \notin & \text{ADirF} \\ (k,i) \in E & \iff & \wedge(\bar{j}, \bar{k}) \notin \text{ADirF} \\ \iff (k,j) \in E & \iff & \wedge(l,i) \notin \text{ADirF} \\ \iff (l,i) \in E & \iff & \wedge(\bar{i}, \bar{l}) \notin \text{ADirF} \\ \iff (l,j) \in E & \iff & \wedge(l,j) \notin \text{ADirF} \\ \wedge(\bar{j}, \bar{l}) & \notin & \text{ADirF} \end{array} $$

### ▷ Proof

The reciprocal are trivial as $\text{ADirF} \subset E$.

All the edges using only $i, \bar{i}, k$ or $\bar{k}$ belong to $\text{ADirF}$ as all the vertex' occurrences are even. While all the edges using $j, \bar{j}, l$ or $\bar{l}$ cannot belong to $\text{ADirF}$ as one of the vertex' occurrence is odd.

▲

### ► Proposition 2.12: Minimum $\text{ADirF}$ set for different identifiers

When the vertex' identifiers are different, it is not possible to reduce the number of canonical edges based on occurrences comparisons.

### ▷ Proof

In order to reduce $\text{ADirF}$ when the vertex' identifiers are different, we have to be more restrictive on the occurrences constraints. To reduce the loops on $k$ and $k'$, we have to find an order between the occurrences of $u$ and $v$.
However, for each order relation between the occurrences ($occ_u \le occ_v$ and $occ_u \ge occ_v$), we can find two associated counter-examples. For the sake of clarity, we will find a counter-example for the case $occ_u \le occ_v$. The reasoning is analogous for the second case.

Let a, b and c be three contigs in C such that all their identifiers are different, $mult_a = 2$, $mult_b = mult(c) = 1$, and $(a_f, b_f) \in \mathcal{L}$, $(b_f, a_f) \in \mathcal{L}$ and $(a_f, c_f) \in \mathcal{L}$.

Let $fa_1 = (i, j)$, $fa_2 = (i', j')$ be the two direct fragments associated to a and $fb = (k, l)$, $fc = (m, n)$ the two respectively associated to b and c. The idea behind the proof is to know if the adjacent direct fragments $fa_1$ $fb$ $fa_2$ $fc$ can participate into a solution path $Path$. In such a case, $Path$ is sub-defined by edges $i \to k \to i' \to m$ (and $j \to l \to j' \to n$). Remark that $k$ can permute with $l$, and $m$ with $n$.

According to ADirF set definition, $(i, k) \in ADirF$, $(k, i') \in ADirF$ and $(i', m) \in ADirF$. If we constrain the occurrences by the relation $occ_u \le occ_v$, it implies that:

$$occ_i \le occ_k (= occ_l - 1)$$

$$occ_l - 1 \le occ_{i'} - 1$$

$$occ_{i'} \le occ_m$$

As $fa_1 \ne fa_2$, $occ_i \ne occ_{i'} - 1$:

if $occ_i < occ_{i'} - 1$ then when we focus on the first edges path (in the "best" case n permutes with m) $occ_k(= 0) \le occ_{i'}(= 2) \le occ_n(= 1)$, so it is absurd;

else $occ_i > occ_{i'} - 1$ then when we focus on the first edges path, you can permute $i'$ and $i$ and you fall into the same absurd implication as above.

Δ

► **Proposition 2.13: Occurences reduction for equal identifiers in ADirF**

There exists an order between the occurrences of the vertices of an edge between two direct fragments that does not reduce the number of feasible and distinct solutions when the vertex' identifiers are equal.
▷ Proof

Let $u_1, u_2, \dots, u_n, n = 2\lfloor\frac{mult_u}{2}\rfloor$ be same identifiers vertices participating in a solution path for $\mathcal{IRP}$, in this order (i.e. $\forall k \in \llbracket 1, n \rrbracket$, $u_k$ is before $u_{k+1}$ in the solution path).

It can be provenª that: regardless the equation set of $k-1$ order relation $R_{k,k+1} \in \{<, >\}$ between two vertices $ctg_{uk}$ and $ctg_{uk+1}$ (i.e. $ctg_{uk}R_{k,k+1}ctg_{uk+1}, \forall k \in \llbracket 1, n \rrbracket$), there exists a permutation between all distinct occurrences of $u$ such that all the relations $R_{k,k+1}$ are satisfied. Thus, let defined the following equation set (used in ADirF):

$$R_{k,k+1} := \begin{cases} < & \text{if } or_{uk} = 0 \lor or_{uk+1} = 0 \\ > & \text{else } or_{uk} = or_{uk+1} = 1 \end{cases}$$

Have we got enough edges to build a path with $\frac{n}{2}$ adjacent direct fragments? On the one hand, the maximum number of edges we would need is: $2(\lfloor\frac{mult_u}{2}\rfloor - 1)$ (the 2 is caused by direct fragments). On the other hand, ADirF provides:

$$2 \sum_{k'=1}^{\lfloor\frac{mult_u}{2}\rfloor-1} 1 = 2\left(\lfloor\frac{mult_u}{2}\rfloor - 1\right)$$

edges, (the 2 is caused by diradj function).

▷

► Proposition 2.14: Canonical reductions for ADirF for equal identifiers

ADirF set contains only one of the eight edges that exist between two adjacent direct fragments when the vertex' identifiers are equal, i.e.

<footnote>ªWith Lucas Robidou we are currently writing a paper and an associated algorithm</footnote>
$$ \forall (i,j) \in DirF, \forall (k,l) \in DirF \mid ctg_i = ctg_k \land occ_i < occ_k: $$

$$ \begin{array}{l} (i,k) \in ADirF \\ \land(\bar{k}, \bar{i}) \notin ADirF \\ \land(i,l) \notin ADirF \\ (i,k) \in E \\ \iff (i,l) \in E \iff \land(\bar{l}, \bar{i}) \notin ADirF \\ \iff (j,k) \in E \iff \land(j,k) \notin ADirF \\ \iff (j,l) \in E \iff \land(\bar{k}, \bar{j}) \notin ADirF \\ \phantom{\iff (j,k) \in E \iff \land(j,k) \notin ADirF} \land(j,l) \notin ADirF \\ \phantom{\iff (j,k) \in E \iff \land(j,k) \notin ADirF} \land(\bar{l}, \bar{j}) \notin ADirF \\ (k,i) \notin ADirF \\ \phantom{\iff (k,i) \notin ADirF} \land(\bar{i}, \bar{k}) \in ADirF \\ (k,i) \in E \\ \iff (k,j) \in E \iff \land(\bar{j}, \bar{k}) \notin ADirF \\ \iff (l,i) \in E \iff \land(l,i) \notin ADirF \\ \iff (l,j) \in E \iff \land(\bar{i}, \bar{l}) \notin ADirF \\ \phantom{\iff (l,i) \in E \iff \land(l,i) \notin ADirF} \land(l,j) \notin ADirF \\ \phantom{\iff (l,i) \in E \iff \land(l,j) \notin ADirF} \land(\bar{j}, \bar{l}) \notin ADirF \end{array} $$

### ▷ Proof

The reciprocals are trivial as $ADirF \subset E$.

Only $(i, k) \in ADirF$ and $(\bar{i}, \bar{k}) \in ADirF$ as the occurrences of $i, \bar{i}, k$ or $\bar{k}$ are even, and $occ_i = occ_{\bar{i}} < occ_k = occ_{\bar{k}}$.

For the other cases, either they contradict the even-occurrences constraint or the constraint on the occurrences order.

◁
Reminder of adjacent inverted fragments set definition:

$$ AInvF = \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ \quad ctg_u < ctg_v \\ \quad \wedge occ_u = 2k + or_u \\ \quad 0 \le k < \left\lfloor \frac{mult_u}{2} \right\rfloor \\ \quad \wedge occ_v = 2k' + or_v \\ \quad 0 \le k' < \left\lfloor \frac{mult_v}{2} \right\rfloor \end{array} \right\} \cup \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ \quad ctg_u = ctg_v \\ \quad \wedge (or_u = 0 \lor or_v = 0) \\ \quad \wedge occ_u - or_u = 2k \\ \quad \wedge occ_v - or_v = 2k' \\ \quad 0 \le k < k' < \left\lfloor \frac{mult_u}{2} \right\rfloor \end{array} \right\} $$

This set is the minimum exhaustive set of canonical edges which represent adjacent inverted fragments. Two combinatorial reductions are necessary — *canonical reduction for different identifiers* (Propositions 2.15 and 2.16) and *same identifiers occurrences* reductions (Propositions 2.17 and 2.18).

► **Proposition 2.15: Canonical reductions for AInvF for different identifiers**

$AInvF$ set contains only one of the four edges that exist between two adjacent inverted fragments when the vertex' identifiers are different, i.e.

$$ \forall (i, j) \in InvF, \forall (k, l) \in InvF \mid ctg_i < ctg_k: $$

$$ \begin{array}{l} (i, k) \in AInvF \\ (i, k) \in E \iff \begin{array}{l} \wedge(\bar{k}, \bar{i}) \notin AInvF \\ \wedge(l, j) \notin AInvF \\ \wedge(\bar{j}, \bar{l}) \notin AInvF \end{array} \\ (k, i) \notin AInvF \\ (k, i) \in E \iff \begin{array}{l} \wedge(\bar{i}, \bar{k}) \notin AInvF \\ \wedge(j, l) \in AInvF \\ \wedge(\bar{l}, \bar{j}) \in AInvF \end{array} \\ (i, l) \in AInvF \\ (i, l) \in E \iff \begin{array}{l} \wedge(\bar{l}, \bar{i}) \notin AInvF \\ \wedge(\bar{j}, \bar{k}) \notin AInvF \end{array} \end{array} $$
$$ (j,k) \in AInvF \\ (j,k) \in E \iff \begin{aligned}[t] & \wedge(\bar{k}, \bar{j}) \notin AInvF \\ & \wedge(l, i) \notin AInvF \\ & \wedge(\bar{i}, \bar{l}) \notin AInvF \end{aligned} $$

### ▷ Proof

The reciprocals are trivial as $AInvF \subset E$.

Let $(i, j) \in InvF$, $(k, l) \in InvF$ be two adjacent inverted fragments such that $ctg_i < ctg_k$. We will focus on demonstrating the case $(i, k) \in E$. The others follow the same logics.

- $(i, k) \in AInvF$ as $\exists n \in [0, \lfloor \frac{mult_i}{2} \rfloor] \lfloor \text{ such that } occ_i = 2n + or_i = 2n, \text{ and } \exists n' \in [0, \lfloor \frac{mult_k}{2} \rfloor] \lfloor \text{ such that } occ_k = 2n' + or_k = 2n' + 1; \}$

- $(\bar{k}, \bar{i}) \notin AInvF$ as $ctg_{\bar{k}} > ctg_{\bar{i}}$;

- $(l, j) \notin AInvF$ as $ctg_l > ctg_j$;

- $(\bar{j}, \bar{l}) \notin AInvF$ as $or_{\bar{j}} = 0$ and $occ_{\bar{j}} \nmid 2$ so $\exists n \in [0, \lfloor \frac{mult_j}{2} \rfloor] \lfloor \text{ such that } occ_{\bar{j}} = 2n. \}$

### △

### ▶ Proposition 2.16: Minimum $AInvF$ set for different identifiers

When the vertex' identifiers are different, it is not possible to reduce the number of canonical edges based on occurrences comparisons.

### ▷ Proof

In order to reduce $AInvF$ when the vertex' identifiers are different, we have to be more restrictive on the occurrences constraints. To reduce the loops on $k$ and $k'$, we have to find an order between the occurrences of $u$ and $v$.

However, for each order relation between the occurrences ($occ_u \le occ_v$ and $occ_u \ge occ_v$), we can find two associated counter-examples. For the sake of clarity, we will find a counter-example for the case $occ_u \le occ_v$. The reasoning is analogous for the second case.

Let $c$ and $d$ be two contigs in $C$ such that $ctg_c < ctg_d$, $mult(c) = 4$ and
$mult_d = 2$, and $(c_f, d_f) \in \mathcal{L}$ and $(d_f, c_f) \in \mathcal{L}$. Let $p_1 = (i, j)$, $p_2 = (i', j')$ be the two inverted fragments associated to $c$ and $q = (k, l)$ the one associated to $d$. The idea behind the proof is to know if the adjacent inverted fragments $p_1$ and $p_2$ can participate into a solution path $Path$. In such a case, $Path$ is sub-defined by edges $i \to k \to i'$ and $j' \to l \to j$.

According to AInvF set definition, $(i, k) \in AInvF$ and $(j', l) \in AInvF$. If we constrain the occurrences by the relation $occ_u \le occ_v$, thus on the one hand, $occ_i \le occ_k(= occ_l - 1)$, and on the other hand $occ_l - 1 \ge occ_{j'} - 1$. As $p_1 \ne p_2$, $occ_i \ne occ_{j'} - 1$:

if $occ_i < occ_{j'} - 1$ then $occ_k(= 0) > occ_i(= 0)$ so it is absurd;

else $occ_i > occ_{j'} - 1$ then $occ_i(= 2) \le occ_k(= 0)$ so it is absurd.

Δ

► **Proposition 2.17: Occurences reduction for equal identifiers in AInvF**

There exists an order between the occurrences of the vertices of an edge between two inverted fragments that does not reduce the number of feasible and distinct solutions when the vertex' identifiers are equal.

▷ Proof

Let $u_1, u_2, \dots, u_n, n = 2\lfloor \frac{mult_u}{2} \rfloor$ be same identifiers vertices participating in a solution path for $\mathcal{IRP}$, in this order (i.e. $\forall k \in [1, n]$, $u_k$ is before $u_{k+1}$ in the solution path).

It can be provenª that: regardless the equation set of $k - 1$ order relation $R_{k,k+1} \in \{<, >\}$ between two vertices $ctg_{uk}$ and $ctg_{uk+1}$ (i.e. $ctg_{uk}R_{k,k+1}ctg_{uk+1}, \forall k \in [1, n]$), there exists a permutation between all distinct occurrences of $u$ such that all the relations $R_{k,k+1}$ are satisfied. Thus, let define the following equation set (used in AInvF):

$$R_{k,k+1} := \begin{cases} < & \text{if } or_{uk} = 0 \lor or_{uk+1} = 0 \\ > & \text{else } or_{uk} = or_{uk+1} = 1 \end{cases}$$

Have we got enough edges to build a path with $\frac{n}{2}$ adjacent inverted fragments? On the one hand, the maximum number of edges we would need is: $2(\lfloor \frac{mult_u}{2} \rfloor - 1)$ (the 2 is caused by inverted fragments). On the other
hand, AInvF provides:

$$2 \sum_{k'=1}^{\lfloor \frac{mult_u}{2} \rfloor - 1} 1 = 2 \left( \lfloor \frac{mult_u}{2} \rfloor - 1 \right)$$

edges, (the 2 is caused by invadj function).

△

► **Proposition 2.18: Canonical reductions for AInvF for equal identifiers**

AInvF set contains only one of the four edges that exist between two adjacent inverted fragments when the vertex' identifiers are different, i.e.

$$\forall (i, j) \in InvF, \forall (k, l) \in InvF \mid ctg_i = ctg_k \land occ_i < occ_k$$

$$\begin{align*}
(i, k) &\in AInvF \\
(i, k) \in E &\iff \begin{array}{l} \land(\bar{k}, \bar{i}) \notin AInvF \\ \land(l, j) \notin AInvF \\ \land(\bar{j}, \bar{l}) \notin AInvF \end{array} \\
(i, l) &\in AInvF \\
(i, l) \in E &\iff \begin{array}{l} \land(k, j) \notin AInvF \\ \land(\bar{l}, \bar{i}) \notin AInvF \\ \land(\bar{j}, \bar{k}) \notin AInvF \end{array} \\
(j, k) &\in AInvF \\
(j, k) \in E &\iff \begin{array}{l} \land(\bar{k}, \bar{j}) \notin AInvF \\ \land(l, i) \notin AInvF \\ \land(\bar{i}, \bar{l}) \notin AInvF \end{array}
\end{align*}$$

Nota bene (k, i) ∈ E case is not here because (i, k) does not change the solution as their sequences are equal. This non-redundancy advantage is allowed by occurrence reduction in Proposition 2.17.

▷ Proof

The reciprocals are trivial as AInvF ⊂ E.

<footnote>ªWith Lucas Robidou we are currently writing a paper and an associated algorithm</footnote>

<footnote>Let (i, j) ∈ InvF, (k, l) ∈ InvF be two adjacent inverted fragments such</footnote>
that $ctg_i = ctg_k$ and $occ_i < occ_k$. We will focus on demonstrating the case $(i, k) \in E$. The others follow the same logics.

- $(i, k) \in AInvF$ as $\exists n \in [0, \lfloor \frac{mult_i}{2} \rfloor] \lfloor such that occ_i = 2n + or_i = 2n, \text{ and} \\ \exists n' \in [0, \lfloor \frac{mult_k}{2} \rfloor] \lfloor such that occ_k = 2n' + or_k = 2n' + 1; \}$

- $(\bar{k}, \bar{i}) \notin AInvF$ as $occ_{\bar{k}} > occ_{\bar{i}}$;

- $(l, j) \notin AInvF$ as $occ_l > occ_j$;

- $(\bar{j}, \bar{l}) \notin AInvF$ as $or_{\bar{j}} = 0$ and $occ_{\bar{j}} \neq 2$ so $\exists n \in [0, \lfloor \frac{mult_{\bar{j}}}{2} \rfloor] \lfloor such that \\ occ_{\bar{j}} = 2n.$

Δ

# 3 Metrics

## 3.1 Quast metrics

Quast metric descriptions can be found at https://quast.sourceforge.net/docs/manual.html#sec3.1. Here we adapt the description of the ones we use in this paper.

\# misassemblies is the number of positions in the contigs (breakpoints) that satisfy one of the following criteria:

- the left flanking sequence aligns over 1 kbp away from the right flanking sequence on the reference;

- flanking sequences overlap on more than 1 kbp;

- flanking sequences align to different strands.

\# local misassemblies is the number of positions in the contigs (breakpoints) that satisfy the following conditions:

- the gap or overlap between left and right flanking sequences is less than 1 kbp, and larger than 200 bp (the maximum indel length);

- the left and right flanking sequences both are on the same strand.

Genome fraction (%) is the percentage of aligned bases in the reference genome. A base in the reference genome is aligned if there is at least one contig with at
least one alignment to this base. Contigs from repetitive regions may map to multiple places, and thus may be counted multiple times.

# 4 Supplementary results

## 4.1 v1 scaffolding benchmark

<table><caption>Table 18 – Benchmark 3 v1 contig Quast</caption><thead><tr><th>Instance</th><th>ILPs</th><th>%gnm</th><th>#mis</th></tr></thead><tbody><tr><td rowspan="2">Abies_alba</td><td>dr-sc</td><td>99.32</td><td>0</td></tr><tr><td>ir-sc</td><td>99.32</td><td>0</td></tr><tr><td>Acorus_americanus</td><td>ir-sc</td><td>99.52</td><td>0</td></tr><tr><td rowspan="2">Agathis_dammara</td><td>dr-ir-sc</td><td>80.19</td><td>0</td></tr><tr><td>ir-dr-sc</td><td>80.19</td><td>0</td></tr><tr><td>Azima_tetracantha</td><td>ir-sc</td><td>99.93</td><td>0</td></tr><tr><td>Begonia_pulchrifolia</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Carpodetus_serratus</td><td>ir-sc</td><td>98.61</td><td>0</td></tr><tr><td>Circaeaster_agrestis</td><td>ir-sc</td><td>99.56</td><td>0</td></tr><tr><td>Clematis_repens</td><td>ir-sc</td><td>99.93</td><td>0</td></tr><tr><td>Commiphora_foliacea</td><td>ir-sc</td><td>98.34</td><td>0</td></tr><tr><td>Cucumis_hystrix</td><td>ir-sc</td><td>99.49</td><td>0</td></tr><tr><td>EucommiaUlmoides</td><td>ir-sc</td><td>99.01</td><td>0</td></tr><tr><td>Jasminum_tortuosum</td><td>ir-sc</td><td>99.77</td><td>0</td></tr><tr><td>Juniperus_scopulorum</td><td>sc</td><td>98.82</td><td>0</td></tr><tr><td>Lamprocapnos_spectabilis</td><td>dr-sc</td><td>35.11</td><td>0</td></tr><tr><td rowspan="2">Lathyrus_pubescens</td><td>dr-sc</td><td>98.20</td><td>0</td></tr><tr><td>ir-sc</td><td>98.20</td><td>0</td></tr><tr><td>Lophocereus_schottii</td><td>sc</td><td>98.63</td><td>0</td></tr><tr><td>Musa_ornata</td><td>ir-sc</td><td>99.08</td><td>0</td></tr><tr><td>Oenothera_glazioviana</td><td>ir-sc</td><td>98.50</td><td>0</td></tr><tr><td>Pelargonium_nanum</td><td>ir-sc</td><td>96.04</td><td>0</td></tr><tr><td>Podocarpus_totara</td><td>sc</td><td>98.76</td><td>0</td></tr><tr><td>Porphyra_purpura</td><td>dr-sc</td><td>99.95</td><td>0</td></tr><tr><td>Sagittaria_trifolia</td><td>ir-sc</td><td>98.55</td><td>0</td></tr><tr><td rowspan="2">SciadopitysVerticillata</td><td>dr-sc</td><td>96.24</td><td>0</td></tr><tr><td>ir-sc</td><td>96.24</td><td>0</td></tr><tr><td>Sciaphila_densiflora</td><td>sc</td><td>99.87</td><td>0</td></tr><tr><td>Selaginella_kraussiana</td><td>dr-sc</td><td>99.84</td><td>0</td></tr><tr><td>Selaginella_vardei</td><td>dr-sc</td><td>99.98</td><td>0</td></tr><tr><td>Taxus_baccata</td><td>sc</td><td>98.07</td><td>0</td></tr><tr><td>Triosteum_pinnatifidum</td><td>ir-dr-sc</td><td>98.30</td><td>0</td></tr><tr><td>Uvaria_macrophylla</td><td>ir-sc</td><td>98.35</td><td>0</td></tr></tbody></table>

(the table continues on the next page)
<table><caption>Table 18, continued</caption><thead><tr><th>Instance</th><th>ILPs</th><th>%gnm</th><th>#mis</th></tr></thead><tbody><tr><td>Welwitschia_mirabilis</td><td>ir-sc</td><td>99.17</td><td>0 0</td></tr><tr><td>Wolffia_australiana</td><td>ir-sc</td><td>99.83</td><td>0 0</td></tr></tbody></table>

<table><caption>Table 19 – Benchmark 3 v1 ILP stats</caption><thead><tr><th>Instance</th><th>|C|</th><th>|L|</th><th>|V|</th><th>|E|</th><th>Time</th></tr></thead><tbody><tr><td>Abies_alba</td><td>9</td><td>28</td><td>20</td><td>36</td><td>0.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.02</td></tr><tr><td>Acorus_americanus</td><td>5</td><td>16</td><td>16</td><td>40</td><td>0.04</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.05</td></tr><tr><td>Agathis_dammara</td><td>20</td><td>54</td><td>50</td><td>96</td><td>0.20</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.64</td></tr><tr><td>Azima_tetracantha</td><td>7</td><td>20</td><td>16</td><td>28</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Begonia_pulchrifolia</td><td>6</td><td>14</td><td>12</td><td>14</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Carpodetus_serratus</td><td>13</td><td>44</td><td>36</td><td>76</td><td>0.14</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.21</td></tr><tr><td>Circaeaster_agrestis</td><td>13</td><td>36</td><td>44</td><td>112</td><td>0.37</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>4.06</td></tr><tr><td>Clematis_repens</td><td>6</td><td>16</td><td>18</td><td>38</td><td>0.04</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.04</td></tr><tr><td>Commiphora_foliacea</td><td>16</td><td>38</td><td>38</td><td>58</td><td>0.04</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.04</td></tr><tr><td>Cucumis_hystrix</td><td>9</td><td>22</td><td>22</td><td>40</td><td>0.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.03</td></tr><tr><td>Eucommia-ulmoides</td><td>13</td><td>32</td><td>32</td><td>52</td><td>0.03</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.05</td></tr><tr><td>Jasminum_tortuosum</td><td>10</td><td>26</td><td>30</td><td>68</td><td>0.07</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.13</td></tr><tr><td>Juniperus_scopulorum</td><td>10</td><td>28</td><td>20</td><td>28</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Lamprocapnos_spectabilis</td><td>12</td><td>36</td><td>26</td><td>46</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Lathyrus_pubescens</td><td>17</td><td>44</td><td>36</td><td>52</td><td>0.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.02</td></tr><tr><td>Lophocereus_schottii</td><td>12</td><td>36</td><td>24</td><td>36</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Musa_ornata</td><td>15</td><td>32</td><td>46</td><td>82</td><td>0.24</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.87</td></tr><tr><td>Oenothera_glazioviana</td><td>11</td><td>30</td><td>30</td><td>56</td><td>0.07</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.08</td></tr><tr><td>Pelargonium_nanum</td><td>11</td><td>32</td><td>34</td><td>108</td><td>0.17</td></tr></tbody></table>

(the table continues on the next page)
<table><caption>Table 19, continued</caption><thead><tr><th>Instance</th><th>|C|</th><th>|L|</th><th>|V|</th><th>|E|</th><th>Time</th></tr></thead><tbody><tr><td>Podocarpus_totara</td><td>13</td><td>42</td><td>26</td><td>42</td><td>0.43</td></tr><tr><td>Porphyra_purpura</td><td>3</td><td>8</td><td>8</td><td>16</td><td>0.00<br>0.03</td></tr><tr><td>Sagittaria_trifolia</td><td>15</td><td>32</td><td>34</td><td>46</td><td>0.00<br>0.00</td></tr><tr><td>SciadopitysVerticillata</td><td>25</td><td>70</td><td>52</td><td>78</td><td>0.02<br>0.02<br>0.04</td></tr><tr><td>Sciaphila_densiflora</td><td>4</td><td>8</td><td>8</td><td>8</td><td>0.00<br>0.00</td></tr><tr><td>Selaginella_kraussiana</td><td>5</td><td>12</td><td>14</td><td>26</td><td>0.01<br>0.01</td></tr><tr><td>Selaginella_vardei</td><td>3</td><td>8</td><td>8</td><td>16</td><td>0.00<br>0.00</td></tr><tr><td>Taxus_baccata</td><td>17</td><td>42</td><td>34</td><td>42</td><td>0.00<br>0.00</td></tr><tr><td>Triosteum_pinnatifidum</td><td>13</td><td>40</td><td>32</td><td>66</td><td>0.07<br>0.15</td></tr><tr><td>UvariaMacrophylla</td><td>13</td><td>30</td><td>42</td><td>86</td><td>0.24<br>0.84</td></tr><tr><td>Welwitschia_mirabilis</td><td>13</td><td>34</td><td>30</td><td>48</td><td>0.02<br>0.03</td></tr><tr><td>Wolffia_australiana</td><td>10</td><td>28</td><td>24</td><td>44</td><td>0.02<br>0.02</td></tr></tbody></table>

## 4.2 v2 scaffolding benchmark

<table><caption>Table 20 – Benchmark 3 v2 ILP stats</caption><thead><tr><th>Instance</th><th>|C|</th><th>|L|</th><th>|V|</th><th>|E|</th><th>Time</th></tr></thead><tbody><tr><td>Agathis_dammara</td><td>20</td><td>54</td><td>52</td><td>108</td><td>0.22<br>0.83</td></tr><tr><td>Begonia_pulchrifolia</td><td>6</td><td>14</td><td>14</td><td>22</td><td>0.00<br>0.00</td></tr><tr><td>Carpodetus_serratus</td><td>13</td><td>42</td><td>36</td><td>74</td><td>0.13<br>0.20</td></tr><tr><td>Jasminum_tortuosum</td><td>10</td><td>28</td><td>30</td><td>70</td><td>0.08<br>0.13</td></tr><tr><td>Lamprocapnos_spectabilis</td><td>12</td><td>36</td><td>48</td><td>180</td><td>0.68<br>2.67</td></tr></tbody></table>

(the table continues on the next page)
<table><caption>Table 20, continued</caption><thead><tr><th>Instance</th><th>|C|</th><th>|L|</th><th>|V|</th><th>|E|</th><th>Time</th></tr></thead><tbody><tr><td>Lathyrus_pubescens</td><td>17</td><td>44</td><td>36</td><td>52</td><td>0.02</td></tr><tr><td>Lophocereus_schottii</td><td>12</td><td>34</td><td>24</td><td>34</td><td>0.00</td></tr><tr><td>Pelargonium_nanum</td><td>11</td><td>32</td><td>36</td><td>128</td><td>0.33</td></tr><tr><td>Podocarpus_totara</td><td>13</td><td>40</td><td>26</td><td>40</td><td>0.84</td></tr><tr><td>Triosteum_pinnatifidum</td><td>13</td><td>40</td><td>32</td><td>66</td><td>0.00</td></tr><tr><td>Triosteum_pinnatifidum</td><td>13</td><td>40</td><td>32</td><td>66</td><td>0.07</td></tr><tr><td>Triosteum_pinnatifidum</td><td>13</td><td>40</td><td>32</td><td>66</td><td>0.18</td></tr></tbody></table>
<figure><img src="image_103.png" /></figure>
<figure><img src="image_104.png" /></figure>

# ACRONYMS

A | B | C | D | F | I | L | M | N | O | P | S | T

## A

**ALE** Assembly Likelihood Estimation 55

## B

**BAC** Bacterial Artificial Chromosome 9, 11, 17

**BFS** Breadth-First Search 20, 50

**BOG** Best Overlap Graph 17, 18

## C

**CSC** Compressed Sparse Column 58

**CSR** Compressed Sparse Row 58

## D

**DBG** De Bruijn Graph 18–21, 38, 53

**DNA** Deoxyribonucleic acid 2–4, 6–9, 11–14, 22, 24, 51, 54

**DR** Direct Repeat 6, 7, 88, 91–93, 111

## F

**FIFO** First In First Out 111

## I

**ILP** Integer Linear Programming 48, 89, 90, 131, 133, 136

**IR** Inverted Repeat 6–8, 88, 91–93, 111, 116, 117, 136

## L

**LIFO** Last In First Out 111, 117
**LSC** Long Single-Copy 6, 7, 88, 93

**M**

**MILP** Mixed-Integer Linear Programming 48, 50

**MIQP** Mixed-Integer Quadratic Programming 51, 54

**MWA2M** Maximum-Weight Acyclic 2-Matching 48, 49

**N**

**NGS** Next Generation Sequencing technologies 10, 18

**O**

**OLC** Overlap-Layout-Consensus 17–22, 24, 43, 86, 134

**ONT** Oxford Nanopore Technology 12, 15

**P**

**PCR** Polymerase Chain Reaction 10, 11

**S**

**S&E** Seed-and-Extend 49, 52, 54

**SBH** Sequencing-by-Hybridisation 18

**SBS** Sequencing By Synthesis 9–11

**SC** Single-Copy 91, 92

**SCS** Shortest Common Superstring 15–17, 20, 48

**SMRT** Single Molecule Run Time technologies 11

**SSC** Short Single-Copy 6, 7, 88

**T**

**TGS** Third Generation Sequencing technologies 11, 17
<figure><img src="image_105.png" /></figure>

# SYMBOLS

Fragment | Graph | Common sets

## Fragment

C Set of contigs 43–45, 48–50, 90–92, 97, 98, 102, 117, A21

f Forward orientation 25, 28, 29, 31, 35, 37, 59, 64, 65, 68, 73, 75, 77, 79–82, see
also forward & {f, r}

$\mathcal{F}_f$ Set of forward fragments 25, 59–61, 63–65, 69

$\mathcal{F}_r$ Set of reverse fragments 25, 59, 61, 63, 65, 69

$\mathcal{F}$ Set of fragments 25–28, 30, 31, 34, 36, 58, 59, 61, 63, 65–67, 70, 82, 84

$\mathcal{L}_f$ (Multi)set of forward links 59, see also $\mathcal{L}$ & $\mathcal{L}_r$

$\mathcal{L}'$ Underlying set of the multiset $\mathcal{L}$ 26–28, 31, 35, 37, see also $\mathcal{L}$ & $\mathcal{F}$

$\mathcal{L}_r$ (Multi)set of reverse links 59, see also $\mathcal{L}$ & $\mathcal{L}_f$

$\mathcal{L}$ (Multi)set of links 26–28, 30–32, 34, 36, 43, 44, 48–50, 58, 59, 61, 63, 65–68, 70, 82, 90–92, 97, 101, 117, 119, 124, 125, A19, G3, see also $\mathcal{L}'$ & $\mathcal{F}$

{f, r} Set of orientations 25, 26, 31, 44, 64, 66, 73–75, 91, 97, 98, 110, 115, 120, A3

$r$ Reverse orientation 25, 28, 29, 31, 35, 37, 64, 65, 68, 77, 80, see also reverse &
{f, r}

$\mathcal{R}_{aw}$ Set of raw reads 15

$\vdash$ Reverse operation 25–27, 29, 32, 35, 44, 45, 55, 59, 61, 63, 65, 68, 69, 72, 73, 76, 78, 79, 81, 91, 92, 101–104, 106, 108–110, 115, 116, 122, 127, A17, see also
reverse

$\Sigma_F$ Fragment label set 25, 59, 60, 62, 68

$\Sigma_{\mathcal{F}_{un}}$ Unoriented fragment label set 25, 59, 64, 66

$\Sigma_L$ Link label set 26

$\Sigma_{\mathcal{L}_{can}}$ Canonical link label set 59, 60, 62, 64, 66, 68
$\mathcal{F}_{un}$ Set of unoriented fragments 25, 31

**Graph**

$G$ Generic graph see also undirected graph, directed graph, bidirected graph, vertex & edge

$DG$ Directed fragment multigraph, $DG = (V, E, \Phi)$ 27–29, 38, 39, 48, 55, 59, 97, see also directed graph

$BG$ Bidirected fragment multigraph, $BG = (V, E, \Phi, \text{attre})$ 30–33, 38, 39, 44, 49, 50, 66, see also bidirected graph

$UG$ Undirected fragment multigraph, $UG = (V, E_F, E_L, \Phi_L)$ 33–36, 38, 39, 44, 46–51, 68, see also undirected graph

$V$ Vertex set 28–32, 34, 35, 46, 59, 60, 62, 64–66, 68, 69, 72, 76, 78–80, 97, G4, see also vertex

$V_f$ Forward vertex set 64, see also vertex & forward

$E$ Edge set 28–32, 46, 59–66, 97, G4, see also edge

$E_F$ Fragment-edge set 34–36, 68, 69, G4, see also edge, fragment & UG

$E_L$ Link-edge set 34–36, 68, G4, see also edge, link & UG

$\Phi$ Incidence function 28–32, 59, 61, 63, 65, 66, 114–116, G5, see also edge

$\Phi_L$ Incidence function for link-edges 34–36, 68, G4, see also edge, link & UG

$N$ Neighbour set 59, 66, 72, 74–77, 80–82, see also edge

$N^-$ Predecessor set 60, 62, 64, 68, 72, 73, 75–77, 79–82, see also edge, $N \& N^+$

$N^+$ Successor set 60, 62, 64, 68, 72–82, see also edge, $N \& N^-$

$N^F$ Fragment extremity neighbour set see also edge, $N \& UGA$

$N^L$ Link neighbour set 68, 69, see also edge, $N \& UGA$

$\Sigma_V$ Vertex label set 60, 62, 64, 66, 68, 70, 73–75, see also vertex

$\Sigma_{V_f}$ Forward vertex label set 64, see also vertex

$\Sigma_E$ Edge label set 60, 62, 64, 66, 68, 73–75, see also edge
*DGA* All oriented fragments multidigraph 59, 61, 63, see also *DG* & directed multi-graph

*DGS* Oriented fragments' successors multidigraph 61, 63, 68, 70, 72, 84, 86, 134, 135, see also *DG* & directed multigraph

*DGF* Forward fragments multidigraph 63, 65, 66, 68, 72, 84, 86, 134, 135, see also *DG* & directed multigraph

*BGU* Unoriented fragments multibigraph 66, 68, 72, 84, 86, 134, 135, see also *BG* & bidirected multigraph

*UGA* All fragments' extremities undirected graph multiungraph 68, 70, see also *UG* & undirected multigraph

*MDCG* Multiplied doubled contig graph, $MDCG = (V, E, vwex)$ 97, 98, 102, 103, 113, 116, 131, 136, see also *DG*, directed graph, contig & link

*RegGraph* Region graph, *RegGraph* = (*Vreg*, *Ereg*, *Φ*) 114, 116, see also *DG* & directed multigraph

*Vreg* Set of oriented region 114–116, *G5*, see also *RegGraph* & *V*

*Ereg* Multiset of links between two oriented regions 114–116, *G5*, see also *RegGraph* & *E*

**Common sets**

*{0, 1}* Set of binaries 104, 107

*N* The set of natural (positive) integer 26, 29, 32, 35, 59, 90, 97, 110, 117–119

*R* The set of real numbers 90, 97, 103, 104, 106, 109

*$\Sigma_{nuc}$* Nucleotide alphabet set *$\Sigma_{nuc}$* = {*A*, *C*, *G*, *T*} 2, 25, 90
<figure><img src="image_106.png" /></figure>

# COMPUTATIONAL TERMS

Complexity | Database | File formats | Programs and packages

## Complexity

*NP* Nondeterministic polynomial time 117

*NP*-complete Nondeterministic polynomial-time complete 20, 22, 45, 47, 48, 90, 119, 121–123, 131, 136

*NP*-hard Nondeterministic polynomial-time hardness 16, 18, 46, 49, 119, 121

*P* Contains all decision problems that can be solved by a deterministic Turing machine using a polynomial amount of computation time 20, 48, 131, 136

## Database

CpGDB Chloroplast Genome Database 123

NCBI National Center for Biotechnology Information 123

PyPI Python Package Index vii, viii, 86, 135

## File formats

GFA Graphical Fragment Assembly (GFA) Format 135

PAF Pairwise mApping Format 135

## Programs and packages

ARACHNE A whole-genome shotgun assembler 17

ARKS Alignment-free linked read genome scaffolding methodology 39

Bambus General-purpose scaffolder based on paired-reads 30, 39, 45

Bambus 2 Metagenome scaffolder 30, 39, 45

BESST Scaffolder based on paired-end data 33, 39

BOSS Scaffolder based on paired-end or mate-pair read set 30, 39
CABOG Celera Assembler with the Best Overlap Graph 18

Canu Canu is a fork of the Celera Assembler, designed for high-noise single-molecule sequencing 17, 18

CBC Open-source mixed integer linear solver 121

Celera A whole genome assembler originally developed at Celera Genomics for the assembly of the human genome 17

chloroExtractor Pipeline for DNA extraction of chloroplast DNA from whole genome plant data 52

Chloroplast assembly protocol A set of scripts for the assembly of chloroplast genomes out of whole-genome sequencing reads 52

Edena De novo assembler for Illumina reads 30

FALCON Diploid aware genome assembler designed for Pacific Biosciences long read data 28

Fast-Plast Chloroplast short-read assembly pipeline 52

Flye De novo assembler for single-molecule sequencing reads, such as those produced by PacBio and Oxford Nanopore Technologies 20

fragScaff Scaffolder based on linked-read data 33, 39

GapFiller Paired-read based gap filler 54

GAT Scaffolder based on paired-reads and global optimisation method 28, 39, 52

GetOrganelle Organelle genome assembly toolkit x, 52, 89, 132, 137

GRASS Generic short-read scaffolder based on paired-reads or on a reference genome 30, 39

Gurobi Linear, integer and non-linear solver 121, 122

hifiasm Haplotype-resolved de novo assembler initially designed for PacBio HiFi reads 18, 28

HINGE OLC long read assembler based on an idea called *hinging* 28

IOGA Iterative organellar genome assembly 52

Jellyfish *k*-mer counting tool 52
khloraascaf Chloroplast contig scaffolder aware of structural haplotypes based on
several integer linear programs viii, 121–123, 126, 127, 129–132, 136, 137

LACHESIS Chromosome-scale scaffolder based on chromatin interactions 30, 39

LRScaf Hybrid scaffolder based on long reads 33, 39

mdBG Minimizer-space de Bruijn graphs for whole-genome assembly 20

Minia Short-read assembler based on a de Bruijn graph 124

Minimus Assembly pipelines designed specifically for small data-sets, such as the set
of reads covering a specific gene 30

MIP Paired-read scaffolder based on mixed integer programming 30, 39

NOVOPlasty *De novo* assembler and heteroplasmy/variance caller for short circular genomes 52

Opera Paired-end read scaffolder based on the graph bandwidth problem 30, 39

ORG.Asm Organellar read assembler 52

Quast Quality assessment tool for genome assemblies 126, 127, 130, A17

revsymg Python implementation of DNA fragment graph vii, 86, 135

SALSA2 Hi-C scaffolder 33, 39

ScaffMatch Short read scaffolding tool based on Maximum-Weight Matching 33, 39

ScaffoldScaffolder Scaffolder based on a bidirected to directed graph reduction 30

Scaftools Integer linear programming approach for genome scaffolding 33, 39

SCARPA Scaffolder which combines fixed-parameter tractable and bounded algorithms
with linear programming 33, 39

Shasta *De novo* assembler for long reads, optimized for Oxford Nanopore (ONT)
reads 28

SLIQ Simple linear inequalities based mate-pair reads filtering and scaffolding 39

SLR Hybrid scaffolder based on long reads 30, 39

SOAPdenovo Short-read assembly method that can build a *de novo* draft assembly for
the human-sized genomes 20, G10
SOAPdenovo2 Successor of SOAPdenovo 53

SOPRA Scaffolding algorithm for paired reads via statistical optimisation 30, 39

SPAdes St. Petersburg genome assembler — is an assembly toolkit containing various assembly pipelines 20, 54

SSPACE Paired-read scaffolder 30, 39, 54

Velvet Short read *de novo* assembler using de Bruijn graphs 20, 54

wtdbg2 *De novo* sequence assembler for long noisy reads produced by PacBio or Oxford Nanopore Technologies 20
<figure><img src="image_107.png" /></figure>

# GLOSSARY

A | B | C | D | E | F | G | K | L | M | N | O | P | R | S | U | V | W

## A

**alignment** Nucleotide comparison between at least two sequences

**assembly graph** One of the output of genome assembly method *see also* read assembly

## B

**bidirected graph** Graph that contains three edge types 30, 45, 86, 134

**bidirected multigraph** A graph where the same bioriented edge between the same two vertices can occur several times 30, 31, 86, 134, *see also* multigraph & bidirected graph

## C

**chloroplast** Organelle in the plants' cells specialised in the photosynthesis process 2, 4–8, 22, 51, 90, 93, 133

**chromosome** A chromosome is a DNA molecule that contains the genetic information for an organism 4, 13, 21, 24, 133, *see also* genome & DNA

**contig** Result of the assembly of the reads, nucleotide sequence longer than the reads 21, 22, 37, 38, 46–48, 50, 54, 86, 90, 130, 133, 134, 136, *see also* read assembly & scaffolding

## D

**de Bruijn graph** Graph structure where the vertices are word of length *k*, and edges are overlaps of length *k* − 1. Can be vertex-centric or edge-centric. 18–21, 24, 53, *see also* *k*-mer

**directed graph** Graph where the edges are directed 16, 17, 27, 86, 89, 90, 97, 134

**directed multigraph** A graph where the same oriented edge between the same two vertices can occur several times 28, 29, 59, 86, 97, 114, 134, *see also* multigraph & directed graph

## E
**edge** Component of a graph that connects two vertices 17, 22, 28–31, 33–35, 46, 48, 49, 58–67, 69, 70, 72, 76–82, 86, 134, 135, see also graph & vertex

F

**forward** Original nucleotide sequence orientation 14, 25, 33, 34, 59, 62–65, 73, 75, 76, 90, see also reverse

**fragment** Generic nucleotide sequence 2, 8–10, 12–14, 22, 24–28, 30, 31, 33–35, 58, 59, 61–63, 65–70, 73–77, 81, 82, 85, 86, 133, 134

**fragment assembly** Generic term for the DNA fragment assembly, can also denote the whole process of assembling the reads to obtain the whole DNA sequences 8, 13, 15, 16, 22, 37, 51, 133, see also read assembly & scaffolding

**fragment-edge** Component of a graph that connects two vertices 33–35, 68–70, 86, 134, see also graph & vertex

**fragment-neighbour** Vertex connected to a given one forming a fragment-edge 68, see also neighbour, fragment-edge & UG

G

**genome** Entire set of DNA instructions found in a cell 2–10, 13–17, 22, 90, 133, see also DNA

**graph** Object composed of vertices connected by edges 22, 28, 58, 70, 72, 85, 86, 133–135, see also vertex & edge

K

**k-mer** Word of length *k* 18–21, 52–54, 124, see also read assembly, SBH,DBG & de Bruijn graph

L

**link** Ordered pair of oriented fragment 22, 28, 30, 33, 34, 38, 58–70, 73, 74, 86, 91, 134, see also fragment

**link-edge** Component of a graph that connects two vertices 33–35, 46, 48, 49, 68, 70, 86, 134, see also graph & vertex

**link-neighbour** Vertex connected to a given one forming a link-edge 68–70, see also neighbour, link-edge & UG
M

**mitochondrion** Organelle found in the cells of most eukaryotes specialised in the
aerobic respiration 4, 5, 8

**multigraph** A graph where there can be several edges between the same two vertices,
*see also graph*

N

**neighbour** Vertex connected to a given one 58, 61, 63, 65-68, 70, 72, 78-80, 86, 134

**nucleotide** Basic building block of nucleic acids (RNA and DNA) 2, 3, 6, 8-12, 16,
38, *see also* DNA

O

**organelle** Specialised subunit, usually within a cell, that has a specific function 4, 5, 8

**overlap** Suffix-prefix alignment type between two sequences 9, 13, 14, 16-18, 48, 86,
133, 134, *see also* alignment & link

**overlap graph** Graph structure that stores overlaps between reads 17, 22, *see also*
overlap

P

**path** A sequence of vertices such that two consecutive vertices are connected by an
edge in the graph 17, 29, 32, 33, 35, 36, *see also* walk, graph, vertex & edge

**predecessor** In-vertex of one edge of one given vertex 59-62, 64, 65, 68, 72-75, *see*
*also successor*

R

**read** DNA fragment from one DNA's strand, output by a sequencer 9-18, 20, 21, 24,
37, 38, 51-53, 58, 86, 133, 134, *see also sequencer*

**read assembly** Fragment assembly stage to obtain longer nucleotide sequence (contigs) from the reads 22, 28, 33, 38, 47, 48, 53, 133, *see also* fragment assembly & contig

**reverse** Reverse-complement nucleotide sequence orientation 2, 3, 14, 25-31, 33, 34,
59-64, 66, 67, 69, 73, 75, 78-81, 86, 90, 134, *see also forward*
S

**scaffold** Sequence of oriented contigs separated by gap **21, 37, 133**, see *also* **scaffolding**

**scaffolding** Orienting and ordering the contigs **22, 28, 33, 37, 38, 86, 133, 134**, see *also* **fragment assembly, read assembly, contig & scaffold**

**sequencer** Sequencing technology machine **8**

**sequencing** Method to generate nucleotide fragments **2, 8–11, 13, 15, 24, 51**

**strand** The DNA molecule is made up of two strands, each of which has a complementary sequence to the other **2, 3, 8–10, 12**

**string graph** Graph structure where vertices are genomic regions, and edges overlaps between them **17**, see *also* **overlap graph & overlap**

**structural haplotype** Copies of the genome with different structures **7, 22, 51, 55, 136**

**successor** Out-vertex of one edge of one given vertex **59–65, 68, 72, 74, 75**, see *also* **predecessor**

U

**undirected graph** Graph where the edges are undirected **33, 35, 45, 86, 134**

**undirected multigraph** A unoriented graph where there can be several edges between the same two vertices **31, 34, 66, 68**, see *also* **multigraph & directed graph**

**unitig** Special case of a contig, it represents a non-ambiguous assembled sequence **18, 19, 21, 28, 48, 54**, see *also* **contig, read assembly & scaffolding**

V

**vertex** Component of a graph **17, 22, 28–31, 33, 34, 36, 46, 49, 51, 58–70, 72–82, 86, 134, 135**, see *also* **graph & edge**

W

**walk** A sequence of vertices such that two consecutive vertices are connected by an edge in the graph **35**, see *also* **path, graph, vertex & edge**
<figure><img src="image_108.png" /></figure>
**Titre :** Assemblage de fragments ADN : structures de graphes et échafaudage de génomes de chloroplastes

**Mot clés :** Assemblage de génomes, Programmation linéaire en nombres entiers, Répétitions génomiques, Haplotypes structuraux

**Résumé :** L'obtention de la séquence nucléotidique d'une molécule ADN nécessite sa fragmentation par des technologies de séquençage et l'assemblage des fragments. Ces fragments sont appelés lectures. Elles souffrent d'erreurs de séquençage et sont considérées sous deux orientations : celle de leur brin ADN d'origine ou l'inverse-complémentaire pour l'autre brin. L'assemblage se base sur des chevauchements deux à deux entre des lectures orientées, et est composé de trois phases : l'assemblage des lectures pour obtenir des contigs (des séquences plus longues que les lectures), l'échafaudage des contigs,

pour obtenir des échafaudages (des ordres de contigs orientés), et la complétion des échafaudages (trouver les séquences de nucléotides séparant les contigs orientés dans les échafaudages).

Dans ce manuscrit, nous comparons des structures de graphes représentant des relations de successions entre des séquences ADN orientées, utiles à différentes phases de l'assemblage. Puis, nous nous penchons sur le problème de l'échafaudage dédié aux génomes de chloroplastes en proposant une nouvelle formulation, une résolution exacte et une implémentation.

**Title:** DNA fragment assembly: graph structures and chloroplast genome scaffolding

**Keywords:** Genome assembly, Integer linear programming, Genome repeats, Structural haplotypes

**Abstract:** To obtain the nucleotide sequence of a DNA molecule, the molecule is fragmented using a sequencing technology and the fragments are assembled. These fragments are called reads. They are subject to sequencing errors and must be considered in two orientations: that of their original DNA strand, or the reverse-complementary for the other strand. Assembly is based on pairwise overlaps between oriented reads and consists in three phases: assembling the reads to obtain contigs (sequences longer than the reads), scaffolding the contigs to obtain scaffolds (orders of oriented contigs), and completing the scaffolds (finding the nucleotide sequences separating the oriented contigs in the scaffolds).

In this manuscript, we compare graph structures representing succession relations between oriented DNA sequences, useful at different phases of assembly. Then, we address the scaffolding problem dedicated to chloroplast genomes by proposing a new formulation, an exact resolution and an implementation.