# Global exact optimisations for chloroplast structural haplotype scaffolding

Victor Epain, Rumen Andonov

► To cite this version:

Victor Epain, Rumen Andonov. Global exact optimisations for chloroplast structural haplotype scaffolding. Algorithms for Molecular Biology, 2024, 19 (5), pp.1-35. 10.1186/s13015-023-00243-1. hal-04134429v3

HAL Id: hal-04134429

https://inria.hal.science/hal-04134429v3

Submitted on 5 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Global exact optimisations for chloroplast structural
haplotype scaffolding

Victor Epain¹ and Rumen Andonov¹*

¹*GenScale, Centre Inria de l'Université de Rennes, IRISA, 263 Avenue Général Leclerc, Rennes, 35700, France.

*Corresponding author(s). E-mail(s): rumen.andonov@irisa.fr;

**Abstract**

**Background** Scaffolding is an intermediate stage of fragment assembly. It consists in orienting and ordering the contigs obtained by the assembly of the sequencing reads. In the general case, the problem has been largely studied with the use of distances data between the contigs. Here we focus on a dedicated scaffolding for the chloroplast genomes. As these genomes are small, circular and with few specific repeats, numerous approaches have been proposed to assemble them. However, their specificities have not been sufficiently exploited.

**Results** We give a new formulation for the scaffolding in the case of chloroplast genomes as a discrete optimisation problem, that we prove the decision version to be *NP*-Complete. We take advantage of the knowledge of chloroplast genomes and succeed in expressing the relationships between a few specific genomic repeats in mathematical constraints. Our approach is independent of the distances and adopts a genomic regions view, with the priority on scaffolding the repeats first. In this way, we encode the structural haplotype issue in order to retrieve several genome forms that coexist in the same chloroplast cell. To solve exactly the optimisation problem, we develop an integer linear program that we implement in Python3 package khloraascaf. We test it on synthetic data to investigate its performance behaviour and its robustness against several chosen difficulties.

**Conclusions** We succeed to model biological knowledge on genomic structures to scaffold chloroplast genomes. Our results suggest that modelling genomic regions is sufficient for scaffolding repeats and is suitable for finding several solutions corresponding to several genome forms.

**Keywords:** Genome assembly, Inverted repeats, Integer linear programming, NP-Complete problem

# 1 Background

DNA molecule is a support of living mechanism information found in all the living organisms. Its sequence can be seen as a word on the alphabet $\Sigma_{nuc} = \{A, C, G, T\}$. Genetic, genomic and epi-genetic analysis lead to combinatorial problems that need to be solved by computing approaches and thus require to pass from the DNA molecule

to a word representation in a computer. This pro-
cess is described as sequencing the molecule using
a sequencing technology.

*The scaffolding in the process of fragment assembly*

Current technologies are not yet able to read an
entire DNA molecule. They output a huge amount
of small overlapping erroneous DNA sequences,
named *reads*. Moreover, due to the double-strand structure of the DNA, the reads come from either one strand or its complement in reverse order, and the sequencing technologies cannot assure that two reads were sequenced from the same strand. Thus, each read must be considered in two orientations: the one given by the technologies (defined as the *forward* orientation), and its reverse-complement (defined as the *reverse* orientation) (e.g. *ATGCCA* and *TGGCAT* are each other's reverse-complement).

From the reads, genome assembly methods aim to find the longest true DNA sequences. Note that it is sufficient to only find one strand out of two, as the other is obtained by the reverse-complement transformation of the first one. Assembly methods are often split into three major stages: (i) *assembly* of reads based on their overlaps to obtain longer sequences (*contigs*); (ii) *scaffolding*, that aims to obtain an order of oriented contigs, potentially separated by nucleotide distances (*scaffolds*); (iii) *gap-filling*, that aims to complete the assembly by filling the gaps between the contigs in the scaffolds. Here we focus on the scaffolding step uniquely.

The vast majority of proposed scaffolding formulations are based on distances data between the reads (paired-end or mate-pair data) that are adapted for the contigs to obtain scaffolds. Read distances can either be counted to represent a confidence linking information between two contigs (Chateau and Giroudeau, 2015; Mandric and Zelikovsky, 2015), either be considered precisely for contig nucleotide positioning (Andonov et al., 2019), or combined both approaches as in (Huson et al., 2002; Wang et al., 2002; Salmela et al., 2011).

### Chloroplast genome architecture and structural haplotypes

In this paper, we address the scaffolding problem for the particular class of chloroplast genomes. 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 loosed 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 long and short single-copies (*LSC, 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 1 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 IR on either side of the reversed subsequence.

### Chloroplast scaffolding approaches

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.

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
<figure><img src="image_2.png" /><figcaption>**Figure 1** The most studied chloroplast genome form is circular and very often quadripartite. For each of the three figures, coloured arrows represent nucleotide sequences. *LSC* and *SSC* stand for long and short single copies (purple and red), respectively. They correspond to regions (subsequences) that are not repeated in the genome. On the opposite, *IR* and *DR* stand for inverted and direct repeat (green and blue), respectively. (a) This architecture is the most common one and is defined as a quadripartite architecture. The two green *IR* arrows face each other and illustrates that one is the reverse-complement sequence of the other; (b) the two blue *DR* arrows are in the same direction that illustrates both have the same nucleotide sequence; (c) the two types of repeat can simultaneously exist in the chloroplast genome, and DRs are shorter than IRs.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2** During the DNA replication of the chloroplast genome, one of the region between the two inverted repeats can be reversed (c.f. the red region *SSC*). This provokes the existence of several forms of the genome in the same chloroplast (*heteroplasmy*).</figcaption></figure>

of a repeat have identical (or reversed) nucleotide sequence; (iii) structural haplotypes can be seen as permutations in a sequence of oriented contigs.

## Organisation of the paper

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 directed graph (digraph) where we apply several Integer Linear Programming (ILP) strategies (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

<table><caption>**Table 1** Toy example of input data. Left: set of contigs *C*. Right: set of links *L*. For the sake of space, for no one of the links in the table, its reverse is given, although it belongs to *L*.</caption><thead><tr><th>contig</th><th>mult</th><th>wex</th><th>contig</th><th>orient</th><th>contig</th><th>orient</th></tr></thead><tbody><tr><td>a</td><td>1</td><td>0.70</td><td>a</td><td>f</td><td>c</td><td>r</td></tr><tr><td>b</td><td>2</td><td>0.83</td><td>a</td><td>r</td><td>c</td><td>r</td></tr><tr><td>c</td><td>2</td><td>0.17</td><td>b</td><td>r</td><td>c</td><td>f</td></tr><tr><td>d</td><td>1</td><td>0.43</td><td>b</td><td>f</td><td>d</td><td>f</td></tr><tr><td></td><td></td><td></td><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*: $\mathcal{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 1 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. 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 1 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 1 provides the properties the sequence of oriented contigs must respect.

**Definition 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 \text{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 directed repeats, the inverted repeats and the single-copies.

**Definition 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} = (\bar{c}_{n-1}, \dots, \bar{c}_1, \bar{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 1.

**Definition 3 (Direct repeat – DR).** A DR is a couple of regions ($dr_i, dr_j$) where $dr_i = dr_j$.

**Definition 4 (Inverted repeat – IR).** An IR is a couple of regions ($ir_k, ir_l$) where $ir_l = \bar{ir}_k$.

**Definition 5 (Repeat).** A repeat is the generic term to denote either DR or IR. The length replen($R$) of a repeat $R = (r_i, r_j)$ equals the lengths of $r_i$ and $r_j$ ($replen(R) = |r_i| + |r_j|$).

**Definition 6 (Single-copy – SC).** A SC is a region that is not part of a repeat.

**Definition 7 (Region weight).** The weight *rwex*($r$) of a region $r$ is defined as $\text{rwex}(r) = \sum_{c \in r} \text{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 8 (Sequence of oriented regions).** Consider a sequence $SOR = (r_0, r_1, \dots, r_{n-1})$:

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

$$ - \forall c \in \mathcal{C}, \sum_{r \in 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

<footnote n="1">¹The notation $x[i]$ denotes the $i^{th}$ element in $x$. Indices start from 0.</footnote>
are longer than 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 3 illustrates the impact of degenerations on quadripartite structures. Indeed, in Figure 3b, 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.

**Definition 9 (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 10); (ii) IRP for the inverted repeats (Definition 11) and (iii) SCP for the single-copies (Definition 12). We tackle CHSP in a hierarchical succession of DRP, IRP and SCP (Definition 14). Any intermediate problem must preserve the regions found by its predecessors (Definition 13).

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 10 (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$;

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

**Definition 11 (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 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$;

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

Figure 4 provides examples of valid oriented contig positioning for each common genome structure (Figure 1). Although Figure 8 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 12 (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.
<figure><img src="image_4.png" /><figcaption>**Figure 3** 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.</figcaption></figure>

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

**Definition 13** (Chloroplast scaffolding problem succession). *DRP*, *IRP* and *SCP* (Definitions 10 to 12) 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 1 and 4. 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 14** (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.

# 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
<figure><img src="image_5.png" /><figcaption>**Figure 4** Chloroplast repeat scaffolding. 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_0, a_1, b_0, b_1, c_0, c_1$ and $d_0, d_1$ 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_2$) since the IR contains more contigs than the DR.</figcaption></figure>

of the scaffolding problems from Definitions 10 to 12.

**Definition 15** (Multiplied Doubled Contig Graph – MDCG). Given a set of contigs $C$, their multiplicities and the link set $\mathcal{L}$, the multiplied doubled contig graph MDCG = ($V, E, vwex$) is defined such that:

$$V = \left\{ \begin{array}{l} v_{f,0}, \dots, v_{f,n-1}, \\ v_{r,0}, \dots, v_{r,n-1} \end{array} \middle| \begin{array}{l} c \in C \\ n = \text{mult}(c) \end{array} \right\}$$

is the set of all the forward and reverse occurrences of all the contigs ($|V| = 2 \sum_{c \in C} \text{mult}(c)$). The vertices are associated with four functions:

contig: $V \to 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) = vwex(\text{contig}(v))$.

$$E = \left\{ (u, v) \in V^2 \text{ s.t. } \begin{array}{l} \text{contig}(u) \text{ vor}(u) \\ \text{contig}(v) \text{ vor}(v) \end{array} \in \mathcal{L} \right\}$$

is the set of multiplied links ($|E| = \sum_{(c,d) \in \mathcal{L}} \text{mult}(c)\text{mult}(d)$).

Figure 5 illustrates the MDCG representing the example data given in Table 1.

<figure><img src="image_6.png" /><figcaption>**Figure 5** 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 = \text{contig}(c_{r,1})$, in its reverse orientation ($\text{vor}(c_{r,1}) = r$), and in its second occurrence ($\text{vocc}(c_{r,1}) = 1$). The colours are the same as the ones in Figure 1a 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.</figcaption></figure>

## 4.2 Repeated fragment sets

A repeat is a couple containing two identical (or reverse for IR) sequences of oriented contigs (Definitions 3 and 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 16.** 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 5. We then precise the set of repeated fragments for each repeat type. Denote by $\mathcal{R} = \{c \in \mathcal{C} | \text{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 $\mathcal{C}$, i.e. $\forall c, d \in \mathcal{C}, c \neq d \iff c < d \lor c > d$.

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

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

**Definition 18.** An inverted fragment $(u, v) \in V^2$ is a repeated fragment, such that the orientations of $u$ and $v$ differ.

$$ \mathrm{InvF} = \bigcup_{c \in \mathcal{R}} \left\{ \begin{array}{l} (i, j) \in V^2 \text{ s.t.} \\ ctg_i = 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\} $$

Figure 6 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 \to DirF$ and *invfrag*: $V' \subset V \to InvF$ (abstracted with the *repfrag* function, c.f. Section A for their definitions).

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

**Definition 19** (Set of pairs of direct fragments).

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

**Definition 20** (Set of pairs of inverted fragments).

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

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

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

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

**Definition 22.** An adjacent direct fragment is an edge between two direct fragments. Let 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{\text{mult}_u}{2} \right\rfloor \\ \wedge occ_v = 2k' \\ 0 \le k' < \left\lfloor \frac{\text{mult}_v}{2} \right\rfloor \end{array} \right\} $$
<figure><img src="image_7.png" /><figcaption>(a) Direct fragments and their pairs</figcaption></figure>

<caption>(b) Inverted fragments and their pairs</caption>

<figure><img src="image_8.png" /><figcaption>**Figure 7** Adjacent repeated fragment sets examples. 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})$.</figcaption></figure>

$$ \bigcup \left\{ \begin{array}{l} (u, v) \in E \text{ s.t.} \\ \quad ctg_u = ctg_v \\ \quad \wedge (or_u = f \lor or_v = f) \\ \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\} $$

**Definition 23.** 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.} \\
ctg_u < ctg_v \\
\quad \land occ_u = 2k + or_u \\
\quad 0 \le k < \left\lfloor \frac{mult_u}{2} \right\rfloor \\
\quad \land 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.} \\
ctg_u = ctg_v \\
\quad \land (or_u = f \lor or_v = f) \\
\quad \land occ_u - or_u = 2k \\
\quad \land occ_v - or_v = 2k' \\
\quad 0 \le k < k' < \left\lfloor \frac{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 7. In addition, we add two
functions to retrieve the adjacent repeated frag-
ments from the edges in Θ(1): diradj: E → E
and invadj: E → E (abstracted with the repadj
function, c.f. Section A for their definitions).

5 Integer Linear
Programming (ILP)
formulation

Modelling DRP, IRP and SCP from Defini-
tions 10 to 12 requires finding a valid circuit in
MDCG.

**Definition 24** (Valid circuit in MDCG). Given a graph MDCG = (V, E) and a starting vertex s, where ctgs is the starting contig, ors = f and occs = 0. A circuit cp in MDCG is valid if:

- it starts and ends with s;

- ∀ v ∈ cp, v̄ /∈ cp, where ctg_v = ctg_v̄, or_v = 1 - or_v̄ and occ_v = occ_v̄;

- consecutive vertices u and v in cp are connec-
ted by an edge (u, v) ∈ 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} 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 cir-
cuit 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 M x_e \quad \forall e \in E \quad (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_{v\bar{s}} = 0 \qquad \forall v \in N_{\bar{s}}^{-} \quad (C4)
$$

$$
x_{\bar{s}v} = 0 \qquad \forall v \in N_{\bar{s}}^{+} \quad (C5)
$$

$$
\forall v \in V \setminus \{s, \bar{s}\} : \\
i_v + i_{\bar{v}} \le 1 \tag{C6}
$$
$$ \sum_{u \in N_v^-} x_{uv} \leq i_v \leq \sum_{w \in N_v^+} x_{vw} \quad (C7) $$

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

$$ \begin{aligned} x_e &\in \{0, 1\} & \forall e \in E \\ i_v &\in [0, 1] & \forall v \in V \setminus \{s, \bar{s}\} \\ f_e &\in \mathbb{R}_{\ge 0} & \forall e \in E \end{aligned} $$

## 5.2 Repeated regions constraints

The following constraints are general to define ILPs for *DRP* and *IRP*. Definitions 10 and 11 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 10 and 11, and given *PDirF* and *PInvF* (Definitions 19 and 20), denote by *ForbidDR* and *ForbidIR* the sets of forbidden quartet vertices for the DRs and IRs, respectively:

$$ \begin{aligned} \text{ForbidDR} &= \left\{ \begin{array}{l} (i, k, l, j) \\ (k, i, j, l) \end{array} \middle| \begin{array}{l} i, j \in p, i \neq j \\ \wedge k, l \in q, k \neq l \\ \forall (p, q) \in \text{PDirF} \end{array} \right\} \\ \text{ForbidIR} &= \left\{ \begin{array}{l} (i, k, j, l) \\ (k, i, l, j) \end{array} \middle| \begin{array}{l} i, j \in p, i \neq j \\ \wedge k, l \in q, k \neq l \\ \forall (p, q) \in \text{PInvF} \end{array} \right\} \end{aligned} $$

Figure 8 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 respectively associated with *ForbidDR* and

*ForbidIR* sets, such that:

$$ \text{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 \text{PDirF} \end{array} \right\} $$

$$ \text{AlphaIR} = \left\{ \begin{array}{l} (i, k), (i, l), (j, k), (j, l) \\ \text{s.t. } ((i, j), (k, l)) \in \text{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 2 gives the correspondence of the sets depending on the problem to solve.

<table><caption>**Table 2** ILP sets and functions corresponding table.</caption><thead><tr><th>ILP</th><th>RepF</th><th>PRepF</th><th>ARepF</th></tr></thead><tbody><tr><td>DRP</td><td>DirF</td><td>PDirF</td><td>ADirF</td></tr><tr><td>IRP</td><td>InvF</td><td>PInvF</td><td>AlnvF</td></tr></tbody></table><table><thead><tr><th>ILP</th><th>Forbid</th><th>Alpha</th><th>repfrag</th><th>repadj</th></tr></thead><tbody><tr><td>DRP</td><td>ForbidDR</td><td>AlphaDR</td><td>dirfrag</td><td>diradj</td></tr><tr><td>IRP</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 existing flow to define the position $pos(v)$ of a vertex $v \in V$, i.e. $pos(v) = \sum_{w \in N_v^+} f_w$.
<figure><img src="image_9.png" /><figcaption>**Figure 8** Non-exhaustive illustrations for authorised and forbidden order cases for two repeated fragments ((i, j), (k, l)) ∈ PRepF. (a) and (b) Authorised and forbidden orders for PDirF; (c) and (d) Authorised and forbidden orders for PInvF.</figcaption></figure>

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 \text{RepF}: \qquad m_{ij} \le i_i \qquad (C9) $$

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

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

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

$$ \text{pos}(u) + \text{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}: $$

$$ 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 ADirF^* \cup AInvF^* : \\ 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:

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

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

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

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

$$ m_q \le m_p \quad \forall (p, q) \in \text{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 an 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 \text{Repeats}} \text{replen}(r) - |\text{Repeats}| $$

$$ \begin{aligned} & = \sum_{p \in \text{RepF}} 2m_p - \left( \sum_{p \in \text{RepF}} m_p - \sum_{e \in \text{ARepF}} \text{isadj}_e \right) \\ & = \sum_{p \in \text{RepF}} m_p + \sum_{e \in \text{ARepF}} \text{isadj}_e \end{aligned} $$

where *Repeats* is the set of repeats.

### DRP and IRP models

$$ \max \sum_{p \in \text{RepF}} m_p + \sum_{e \in \text{ARepF}} \text{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}\}} \text{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 14, two problem combinations are opposed. The combinations are kept depending on the value of the problems' objective functions.

**Definition 25** (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 25 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 3.

<table><caption>**Table 3** 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></tr></thead><tbody><tr><th>0</th><td>0</td><td>-</td><td>sc</td></tr><tr><th>&gt; 0</th><td>0</td><td>-</td><td>dr-sc</td></tr><tr><th>&gt; 0</th><td>&gt; 0</td><td>-</td><td>dr-ir-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 = \begin{pmatrix} \text{rid}_0 & \text{ror}_0 \\ \vdots & \vdots \\ \text{rid}_{2m+n-1} & \text{ror}_{2m+n-1} \end{pmatrix} \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 9 illustrates the regions extracted by Algorithm 4 in the toy example.

<figure><img src="image_10.png" /><figcaption>**Figure 9** Extracting the genome architecture in MDCG from a CHSP solution. 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 4 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.</figcaption></figure>

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 1 aims to determine the 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 2 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 3). At the end, Algorithm 4 builds $COR$ and $SOR$ from an ILP solution.

To build the repeated regions, we use a First In First Out (FIFO) queue for the DRs, and a Last In First Out (LIFO) queue for the IRs. Each queue $rep_queue$, is associated with three methods: $rep_queuePUT(x)$ append $x$ to the FIFO/LIFO;
rep_queue.IS_EMPTY returns true if the queue is empty;

rep_queue.PEEK returns the first/last value in the
FIFO/LIFO;

rep_queue.POP 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 ∈ X?” for an object x in a set X is in Θ(1).
Algorithm 1 is in Θ(1). Algorithm 2 is in O(|SCs|),
where |SCs| is the number of contigs composing
the single-copy region that contains the starting
vertex. Algorithm 3 is in Θ(1), and so Algorithm 4
is in O(|V| + |E|).

**Algorithm 1** 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$) ∈ DirF* **then**

3: return DR

4: if invfrag($v$) ∈ InvF* **then**

5: return IR

6: return SC

8 Multiple genome forms

The sequence of oriented regions SOR repres-
ents one chloroplast genome form. Recall that
especially with the LSC-IR-SSC-IR architecture
(Figure 1a), 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 as-
sembly graph: the region graph. The discovery
of multiple genome forms is associated with the
search of Eulerian circuits in this graph. Figure 10
illustrates the procedure for the toy example.

**Definition 26 (Region graph RegGraph).** Given the $(m+n)$-uplet COR of forward regions and the sequence SOR of oriented regions, RegGraph = $(Vreg, Ereg, \Phi)$ denotes a directed multigraph

**Algorithm 2** Get the initial vertex of the single-copy 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** INITIAL VERTEX( )

2: v ← s

3: u ← u | x_{us}^* = 1

4: reg_code ← GETRegion_code(u)

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

6: v ← u

7: u ← u | x_{uv}^* = 1

8: reg_code ← GETRegion_code(u)

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

10: return s

11: return v

**Algorithm 3** 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 ≠ region_code **then**

3: return False

4: if region_code = DR **then**

5: return $x_{diradj(u,v)}^*$ = 1

6: if region_code = IR **then**

7: return $x_{invadj(u,v)}^*$ = 1

8: return False

(multidigraph) named the region graph, such that:

$$
Vreg = \left\{ \begin{array}{l} v_{0,f}, \dots, v_{m+n-1,f}, \\ v_{0,r}, \dots, v_{m+n-1,r} \end{array} \middle| \left|COR\right| = m + n \right\}
$$

is the set of oriented regions ($|Vreg| = 2|COR|$).
For each vertex, bijective function vreg: Vreg $\leftrightarrow$
$\{r \in COR\} \times \{f, r\}$ provides the oriented region it
represents.
**Algorithm 4 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** ILPSolutionToRegions() 
2:     **prev_region_code** ← SC 
3:     SOR ← [0f] ▷ 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** ← INITIALVertex() 
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** ← GETRegionCode(v) 
12:        **if** region_code = SC **then** 
13:            **if** prev_region_code ≠ region_code **then** 
14:                **region_index** ← |COR| 
15:                SOR.**APPEND**(region_index_f) 
16:                COR.**APPEND**([]) 
17:                COR[region_index].**APPEND**(ctg_v, or_v) 
18:            **else** 
19:                **repf** ← dirfrag(v) if region_code = DR else invfrag(v) 
20:                **if** repf ∉ repf_rep **then** 
21:                    **if** ¬REPEAT_IS_CONTINGUOUS(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**(dirfrag(v) if region_code = DR else invfrag(v)) 
28:                        **repf-rep[repf]** ← region_index 
29:                    **else** 
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
<figure><img src="image_11.png" /><figcaption>**Figure 10** The toy example's solution given in Figure 9 results in the graph visualised in Figure 5 will give a region graph that shows an LSC-IR-SSC-IR architecture as the one given in Figure 1a. 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 2. (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.</figcaption></figure>

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 Ereg$ denotes its reverse where $\Phi(e) = (u, v), \Phi(\bar{e}) = (\bar{v}, \bar{u}) \mid (Ereg = 2|SOR|)$.

$\Phi$: $Ereg \to \{(u, v) \mid 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 \mid \Phi(e) = (vreg^{-1}(r_i), vreg^{-1}(r_j))$.

Figure 10 illustrates the region graph for the toy example's solution given in Figure 9. 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 27** (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 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 1).

*Proof.* Let RegGraph = (Vreg, Ereg, Φ) be a region graph. Denote by $euc = (v_0, v_1, \dots, v_{2m+n-1})$ an Eulerian circuit in RegGraph. For each two consecutive oriented regions $v_i, v_j$ in $euc$, there exists an edge $e \in Ereg$ such that $\Phi(e) = (v_i, v_j)$. According to Definition 26, $vreg(v_i) = r_i$ and $vreg(v_j) = r_j$ are also consecutive in the oriented region sequence that has originally built RegGraph, otherwise in its reverse. Thus, according to Algorithm 4, there is an edge $(u, v) \in E$ in MDCG such that $(ctg_u, or_u)$ and $(ctg_v, or_v)$ equal $r_i[|r_i| - 1]$ and $r_j[0]$, respectively. □

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 10 gives the resulting region graph obtained from the input data given in Table 1.

We accept all the Eulerian circuits, although they may contradict the repeated region interval constraints given in Definitions 10 and 11. For example, the oriented region sequence $(0_f, 1_f, 2_f, 3_f, 1_f, 3_r)$ respects the definitions, where $(1_f, 1_f)$ is a DR and $(3_f, 3_r)$ is an IR. One of the Eulerian circuit produces the oriented region sequence $(0_f, 1_f, 2_f, 3_f, 1_r, 3_r)$. Region $1_f$ now evolves in a new IR $(1_f, 1_r)$. The order between the regions of the two IRs contradicts Definition 11.

## 9 NP-completeness

To prove the $\mathcal{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 28** (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 $\sum_{ir \in IR} replen ir) \ge k$ and $|IR| \le m'$, where IR is the solution set of inverted repeats?
**Proposition 2 (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 accumulative repeat length equals at least *k* and the number of repeats equals at most *m*′.

*Proof.* Algorithm 5 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 accumulative repeat length equals at least *k* and the number of repeats equals at most *m*′.

A LIFO *ir\lifo* is associated with four meth-
ods:

*ir\lifoPUT(x)* append *x* to the LIFO;

*ir\lifo.IS_EMPTY* returns true if the LIFO is empty;

*ir\lifo.PEEK* returns the last value in the LIFO;

*ir\lifo.POP* deletes the last value in the LIFO and returns it.

Algorithm 6 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 5
and Algorithm 6 are linear according the size of
SOR and COR. Remember we assume that the
belonging test “*is* x ∈ X?” for an object x in a
set X is in Θ(1).
□

**Proposition 3.** *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 \mathcal{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

**Algorithm 5** 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* ∈ N.

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

1: **function** IS_SOR_VALID() 
2:     **IR** ← `{}` ▷ Set of inverted repeats 
3:     **count_reg** ← HASHTABLE() ▷ Keep the orientations of a region in a hash table 
4:     ▷ Identify the inverted repeats ▷ 
5:     **for** (*rid*, *ror*) ∈ SOR **do** 
6:         **if** *rid* ∉ **count_reg** **then** 
7:             **count_reg**[*rid*] ← [*ror*] 
8:             **else** 
9:                 **if** |**count_reg**[*rid*] = 2| **then** 
10:                 **return** False 
11:             **if** *ror* ≠ **count_reg**[*rid*][0] **then** 
12:                 **IR**.ADD(*rid*) ▷ Repeats are pairs of regions 
13:             **count_reg**.APPEND(*ror*) 
14:         ▷ Verify parameters *k* and *m*′ 
15:     **if** $\sum_{ir \in IR} |COR[ir]| < k$ or $|IR| > m'$ **then** 
16:         **return** False 
17:     ▷ Verify the order between IRs ▷ 
18:     *ir\lifo* ← LIFO() 
19:     **for** (*rid*, *ror*) ∈ SOR **do** 
20:         **if** *rid* ∈ **IR** **then** 
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() 
□

vertices between *s* and *t* in the longest path). We
shall build an instance transform function *tf* such
that $I \in \mathcal{LPSTP} \iff tf(I) \in \mathcal{DIRP}$. Function
*tf* transforms graph *G* to graph *G'* = (*V'*, *E*'),
vertex *s* ∈ *V* to vertices *s*<sub>*f*</sub>, *s*'<sub>*f*</sub> ∈ *V*', *t* ∈ *V* to
*t*<sub>*f*</sub>, *t*<sub>*r*</sub> ∈ *V*', *k* to *k*' = 2*k*, and fix parameter *m*'<sub>*i*</sub> = 1.
Figure 11 illustrates the transformation.

All four subgraphs $G'_{or,i}$ in Figure 11b can be
seen as copies of $G_{V\setminus\{s,t\}} = (V \setminus \{s,t\}, E_{V\setminus\{s,t\}})$
(the subgraph induced by the vertex set $V \setminus \{s,t\}$)
in Figure 11a. For each $or \in \{f,r\}$, for each $i \in$
$\{0,1\}$, $G'_{or,i} = (V'_{or,i}, E'_{or,i})$ such that:
**Algorithm 6 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: &nbsp;&nbsp;&nbsp;&nbsp;**if** ((`prev_c`, `prev_or`), (`c`, `or`)) ∉ *L* **then**

11: &nbsp;&nbsp;&nbsp;&nbsp;**return** False

12: &nbsp;&nbsp;&nbsp;&nbsp;**if** *c* ∉ `contig_count` **then**

13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`contig_count`[*c*] ← 1

14: &nbsp;&nbsp;&nbsp;&nbsp;**else**

15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`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: &nbsp;&nbsp;&nbsp;&nbsp;**if** `contig_count`[*c*] > `mult`(*c*) **then**

20: &nbsp;&nbsp;&nbsp;&nbsp;**return** False

— there is a bijective function `vtrans_or,i`: $V \setminus \{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 \setminus \{s,t\}} \hookrightarrow E'_{or,i}$ where $\forall (u,v) \in E_{V \setminus \{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:

<figure><img src="image_12.png" /><figcaption>(a) Digraph *G* for *LPSTP*</figcaption></figure>

<caption>(b) Digraph $G'$ for *DIRP*</caption>

**Figure 11** 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 \setminus \{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.

$$— \forall (s, w) \in E:$$

$$\begin{aligned} \text{esttrans}(s, w) = & \{(s_f, vtrans_{f,i}(w)) \in V'^2 \mid i \in \{0,1\}\} \\ & \cup \{(vtrans_{r,i}(w), s'_f) \in V'^2 \mid i \in \{0,1\}\} \end{aligned}$$

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

$$\begin{aligned} \text{esttrans}(u, t) = & \bigcup_{i \in \{0,1\}} \left\{ \begin{pmatrix} vtrans_{f,i}(u) & t_f \\ t_f & vtrans_{r,i}(u) \end{pmatrix} \in V'^2 \right\} \\ & \cup \left\{ \begin{pmatrix} vtrans_{f,i}(u) & t_r \\ t_r & vtrans_{r,i}(u) \end{pmatrix} \in V'^2 \right\} \end{aligned}$$

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 11b, where *InvF* = $\{(i,j) \in V_{f,0} \times V_{r,1} \mid 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 $\mathcal{LPSTP}$.

To conclude, as there is a linear time com-
plexity transform function tf such that $I \in$
$\mathcal{LPSTP} \iff tf(I) \in \mathcal{DIRP}$, $\mathcal{DIRP}$ is
$\mathcal{NP}$-Hard. $\square$

From Propositions 2 and 3 we conclude that
$\mathcal{DIRP}$ is $\mathcal{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 fol-
lowing 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.

## 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 re-
gions: $SC1 - IR - SC2 - \overline{IR}$. In the following,
we run what corresponds to *IRP* computation
in khloraascaf on two types of growing gener-
ated 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 to-
pology as the one illustrated in Figure 5. For
instance, testing perfect artificial data acts as
a control for further tests. Table 4 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 relax-
ation time seems to fit an exponential distribution
as well as for the B&B time, as shown in Figure 12,
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 multipli-
city 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 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.

<footnote n="2">²https://khloraa-scaffolding.readthedocs.io/en/latest</footnote>

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

<footnote n="4">⁴To reproduce the results, please refer to https://khloraascaf-results.readthedocs.io/en/latest/benchmark_5/.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12** Solver running time distributions for perfect artificial data. Points are measured times, the red curves correspond to the best $ae^{bx}$ function applied on IR length axis.</figcaption></figure>

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

### 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, 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

<footnote n="5">⁵For more details, please refer to https://khlaaascaf-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><caption>**Table 5** Gurobi solver metrics on noisy artificial growing data. The column descriptions are the same as the one in Table 4. 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.</caption><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>19 930</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>16 821</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>50 246</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>27 822</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>129 174</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>42 292</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>92 009</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>84 639</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>210 798</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>228 198</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>691 619</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>86 592</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>287 146</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>91 936</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>269 958</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>116 315</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>351 501</td></tr></tbody></table>

(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}|$$

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, \lfloor cov_c / cov_s - 0.1 \rfloor)$. 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 Section C.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 6 provides all the metrics defined above for the 31 instances. The instances are solved very quickly (solver times < 4.5 sec., Table D2). 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.

#### 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 Cumix_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 D1 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.

<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>
<table><caption>**Table 6** 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>Eucommia-ulmoides</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">Sciadopitys CENTICILLATA</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><tr><td>Wolffia_australiana</td><td>ir-sc</td><td>4</td><td>1</td><td>4</td><td>99.99</td><td>0</td><td>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.

**Lamprocapnos_spectabilis** 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.

<footnote n="9">⁹For more details, please refer to https://khloraascaf-results.readthedocs.io/en/latest/benchmark_3/synthetic_data_v2/.</footnote>
Table 7 provides all the metrics obtained by running khloraascaf. The instances are also solved very quickly (solver times < 3 sec., Table D3). 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 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" />.

## 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 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

<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_quast_analysis/.</footnote>

<footnote n="11">¹¹https://pypi.org/project/khloraascaf/</footnote>
<table><caption>**Table 7** Sequence and Quast metrics for the modified synthetic data version. The caption is the same as in Table 6.</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">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>

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-eqally (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.

**Supplementary information.** Supplementary definitions and proofs can be found in Section A and Section B. Metrics for the numerical results can be found in Section C.1, and supplementary results are in Section D.

**Availability of data and materials.** The Python3 package khloraascaf can be installed with PyPI at https://pypi.org/project/khloraascaf/ and the numerical results can be reproduced thanks to the instructions given in https://khloraascaf-results.readthedocs.io/en/latest/benchmark.3/.

## Declarations

**Competing interests.** The authors declare that they have no competing interest.

**Authors' contributions.** Both authors conceived the method, designed the ILP models and the algorithms, gave the proof of NP-completeness, wrote the paper. V.E. implemented the Python3 khloraascaf package and tested it on all the data.

**Acknowledgments.** We are grateful to Dominique Lavenier who provided the data for the synthetic data. Many thanks to Jacques Nicolas for his suggestions about formalism aspects. We are thankful to Abdelkader Ainouche for numerous discussions on chloroplast genomics. Sven Schrinner and Gunnar Klau are acknowledged for
the valuable discussions about the problem's complexity. We are particularly grateful to the reviewers for their careful reading, suggestions and comments, which we found very pertinent.

**Funding.** Not applicable.

# Appendix A Repeated fragment set functions

$$ dirfrag: \{v \in p, \forall p \in DirF\} \rightarrow 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 \in p, \forall p \in InvF\} \rightarrow InvF $$

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

$$ diradj: E \rightarrow 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 \rightarrow 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_v \\ occ_{u'} = 2 \lfloor \frac{occ_v}{2} \rfloor + (1 - or_v) $$

# Appendix B 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.

## B.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.} \\ ctg_i = ctg_j = c \\ \wedge or_i = or_j \in \{f, r\} \\ \wedge occ_i = occ_j - 1 = 2k \\ 0 \le k < \lfloor \frac{\text{mult}(c)}{2} \rfloor \end{array} \right\} $$

**Proposition 4** (Commutative reduction for DirF).

$$ \forall (i, j) \in DirF, (j, i) \notin DirF $$
Proof. Let $(i, j) \in DirF$. $occ_j = occ_i + 1$ so $(j, i) \notin DirF$ $\square$

**Proposition 5** (Pairing reduction for DirF). 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 \land 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 4, $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$. $\square$

**Proposition 6.** 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. $\square$

Reminder of the InvF definition:

$$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 \land or_i = f \land or_j = r \\ \quad \land occ_i = occ_j - 1 = 2k \\ \quad 0 \le k < \left\lfloor \frac{\text{mult}(c)}{2} \right\rfloor \end{array} \right\}$$

**Proposition 7** (Commutative reduction for InvF).

$$\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$ $\square$

**Proposition 8** (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 \land or_k \ne or_l$ be a candidate for the next direct fragment:

- to respect the property given by Proposition 7, $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 7, $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$. $\square$

**Proposition 9.** 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. $\square$

## B.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.} \\ \quad ctg_j < ctg_k \\ \quad \lor \\ \quad ctg_j = ctg_k \land occ_j < occ_k \end{array} \right\}$$

**Proposition 10** (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 11.** Pair of direct fragment set $PDirF$ 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 \mid (p, q) \notin PDirF'$.

Let $p = (i, j)$ and $q = (k, l)$, and let $q_2 = (m, n) \mid (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 $iklmn$. 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.} \\ ctg_j < ctg_k \\ \lor \\ ctg_j = ctg_k \land occ_j < occ_k \end{array} \right\} $$

**Proposition 12** (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 13.** Pair of inverted fragment set $PInvF$ is a minimal set for $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 \mid (p, q) \notin PInvF'$.

Let be $p = (i, j)$ and $q = (k, l)$, and let $q_2 = (m, n) \mid (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 $iklmn$. This implies the inverted 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 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. □

## B.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.} \\ ctg_u \neq ctg_v \\ \land occ_u = 2k \\ \qquad 0 \le k < \left\lfloor \frac{mult_u}{2} \right\rfloor \\ \land occ_v = 2k' \\ \qquad 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 \\ \land (or_u = f \lor or_v = f) \\ \land occ_u = 2k \\ \land occ_v = 2k' \\ \qquad 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 14 and 15) and same identifiers occurrences reductions (Propositions 16 and 17).
**Proposition 14** (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}{l} (i,k) \in A\text{DirF} \\ \wedge(\bar{k}, \bar{i}) \in A\text{DirF} \\ (i,k) \in E \\ \wedge(i,l) \notin A\text{DirF} \\ \iff (i,l) \in E \\ \wedge(l, \bar{i}) \notin A\text{DirF} \\ \iff (j,k) \in E \\ \wedge(j,k) \notin A\text{DirF} \\ \iff (j,l) \in E \\ \wedge(\bar{k}, \bar{j}) \notin A\text{DirF} \\ \wedge(j,l) \notin A\text{DirF} \\ \wedge(l, \bar{j}) \notin A\text{DirF} \\ (k,i) \in A\text{DirF} \\ \wedge(\bar{i}, \bar{k}) \in A\text{DirF} \\ (k,i) \in E \\ \wedge(k,j) \notin A\text{DirF} \\ \iff (k,j) \in E \\ \wedge(\bar{j}, \bar{k}) \notin A\text{DirF} \\ \iff (l,i) \in E \\ \wedge(l,i) \notin A\text{DirF} \\ \iff (l,j) \in E \\ \wedge(\bar{i}, \bar{l}) \notin A\text{DirF} \\ \wedge(l,j) \notin A\text{DirF} \\ \wedge(\bar{j}, \bar{l}) \notin A\text{DirF} \end{array} $$

*Proof.* The reciprocal are trivial as $ADirF \subset E$.

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

**Proposition 15** (Minimum 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 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 $\mathcal{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 A\text{DirF}$, $(k,i') \in A\text{DirF}$ and $(i',m) \in A\text{DirF}$. If we constrain the occurrences by the relation $occ_u \le occ_v$, it implies that:

$$ \begin{array}{l} \text{occ}_i \le \text{occ}_k(= \text{occ}_l - 1) \\ \text{occ}_l - 1 \le \text{occ}_{i'} - 1 \\ \text{occ}_{i'} \le \text{occ}_m \end{array} $$

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. $\square$

**Proposition 16** (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<fnref n="12" /> 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

<footnote n="12">¹²With Lucas Robidou we are currently writing a paper and an associated algorithm</footnote>
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 17** (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.

$$
\forall (i,j) \in \mathrm{DirF}, \forall (k,l) \in \mathrm{DirF} \mid \mathrm{ctg}_i = \mathrm{ctg}_k \wedge \mathrm{occ}_i < \mathrm{occ}_k:
$$

$$
\begin{array}{l @{\qquad} l @{\qquad} l}
(i,k) \in ADirF & & \\
\wedge(\bar{k}, \bar{i}) \notin ADirF & & \\
(i,k) \in E & & \\
\Leftrightarrow (i,l) \in E & & \wedge(\bar{l}, \bar{i}) \notin ADirF \\
\Leftrightarrow (j,k) \in E & \Leftrightarrow & \wedge(j,k) \notin ADirF \\
\Leftrightarrow (j,l) \in E & & \wedge(\bar{k}, \bar{j}) \notin ADirF \\
& & \wedge(j,l) \notin ADirF \\
& & \wedge(\bar{l}, \bar{j}) \notin ADirF \\
& & \\
(k,i) \notin ADirF & & \\
\wedge(\bar{i}, \bar{k}) \in ADirF & & \\
(k,i) \in E & & \\
\Leftrightarrow (k,j) \in E & \Leftrightarrow & \wedge(\bar{j}, \bar{k}) \notin ADirF \\
\Leftrightarrow (l,i) \in E & & \wedge(l,i) \notin ADirF \\
\Leftrightarrow (l,j) \in E & & \wedge(\bar{i}, \bar{l}) \notin ADirF \\
& & \wedge(l,j) \notin ADirF \\
& & \wedge(\bar{j}, \bar{l}) \notin ADirF
\end{array}
$$

Proof. The reciprocals are trivial as ADirF ⊂ 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:

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

This set is the minimum exhaustive set of ca-
nonical edges which represent adjacent inverted
fragments. Two combinatorial reductions are ne-
cessary — *canonical reduction for different identi-
fiers* (Propositions 18 and 19) and *same identifiers
occurrences reductions* (Propositions 20 and 21).

**Proposition 18** (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 \mathrm{InvF}, \forall (k,l) \in \mathrm{InvF} \mid \mathrm{ctg}_i < \mathrm{ctg}_k:
$$

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

*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$ such that $occ_i = 2n + or_i = 2n$, 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 $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 $\not\exists n \in [0, \lfloor \frac{mult_{\bar{j}}}{2} \rfloor] \lfloor$ such that $occ_{\bar{j}} = 2n$.

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

**Proposition 19** (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.

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

**Proposition 20** (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<fnref n="13" /> 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). $\square$

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

**Proposition 21** (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.

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

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

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

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

**Nota bene** $(k, i) \in 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 20.

*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$ 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 n \rfloor$ such that $occ_i = 2n + or_i = 2n$, and $\exists n' \in [0, \lfloor \frac{mult_k}{2} \rfloor] \lfloor n' \rfloor$ 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}} \nmid 2$ so $\exists n \in [0, \lfloor \frac{mult_j}{2} \rfloor] \lfloor n \rfloor$ such that $occ_{\bar{j}} = 2n$.

□

# Appendix C Metrics

## C.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

# Appendix D Supplementary results

## References

Andonov R, Djidjev H, François S, et al (2019) Complete assembly of circular and chloroplast genomes based on global optimization. Journal of Bioinformatics and Computational Biology 17(3):1950014. https://doi.org/10.1142/S0219720019500148

Ankenbrand MJ, Pfaff S, Terhoeven N, et al (2018) chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data. Journal of Open Source Software 3(21):464. https://doi.org/10.21105/joss.00464

Bendich AJ (2004) Circular chloroplast chromosomes: The grand illusion. The Plant Cell 16(7):1661–1666. https://doi.org/10.1105/tpc.160771

Bock R, Knoop V (eds) (2012) Genomics of Chloroplasts and Mitochondria, Advances in Photosynthesis and Respiration, vol 35.
<table><caption>**Table D1** 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>Abies_alba</td><td>dr-sc</td><td>99.32</td><td>0</td></tr><tr><td></td><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>Agathis_dammara</td><td>dr-ir-sc</td><td>80.19</td><td>0</td></tr><tr><td></td><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>ir-sc</td><td>98.61</td><td>0</td></tr><tr><td>Carpodetus_serratus</td><td>ir-sc</td><td>99.56</td><td>0</td></tr><tr><td>Circaeaster_agrestis</td><td>ir-sc</td><td>99.93</td><td>0</td></tr><tr><td>Clematis_repens</td><td>ir-sc</td><td>98.34</td><td>0</td></tr><tr><td>Commiphora_foliacea</td><td>ir-sc</td><td>99.49</td><td>0</td></tr><tr><td>Cucumis_hystrix</td><td>ir-sc</td><td>99.77</td><td>0</td></tr><tr><td>Eucommia-ulmoides</td><td>ir-sc</td><td>99.01</td><td>0</td></tr><tr><td>Jasminum_tortuosum</td><td>ir-sc</td><td>98.82</td><td>0</td></tr><tr><td>Juniperus_scopulorum</td><td>sc</td><td>35.11</td><td>0</td></tr><tr><td>Lamprocapnos_spectabilis</td><td>dr-sc</td><td>98.20</td><td>0</td></tr><tr><td>Lathyrus_pubescens</td><td>dr-sc</td><td>98.20</td><td>0</td></tr><tr><td>Lophocereus_schottii</td><td>ir-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>SciadopitysVerticillata</td><td>dr-sc</td><td>96.24</td><td>0</td></tr><tr><td></td><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><tr><td>Welwitschia_mirabilis</td><td>ir-sc</td><td>99.17</td><td>0</td></tr><tr><td>Wolffia_australiana</td><td>ir-sc</td><td>99.83</td><td>0</td></tr></tbody></table>

Springer Netherlands, Dordrecht, https://doi.org/10.1007/978-94-007-2920-9

Chateau A, Giroudeau R (2015) A complexity and approximation framework for the maximization scaffolding problem. Theoretical Computer Science 595:92–106. https://doi.org/10.1016/j.tcs.2015.06.023

Chikhi R, Rizk G (2012) Space-Efficient and Extract de Bruijn Graph Representation Based on a Bloom Filter. In: Raphael B, Tang J (eds) Algorithms in Bioinformatics. Springer, Berlin, Heidelberg, Lecture Notes in Computer Science, pp 236–248, https://doi.org/10.1007/978-3-642-33122-0_19

Coissac E, Hollingsworth PM, Lavergne S, et al (2016) From barcodes to genomes: Extending the concept of DNA barcoding. Molecular Ecology 25(7):1423–1428. https://doi.org/10.1111/mec.13549

<table><caption>**Table D2** 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>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>Carpodetus_serratus</td><td>13</td><td>44</td><td>36</td><td>76</td><td>0.14</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>Commiphora_foliacea</td><td>16</td><td>38</td><td>38</td><td>58</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>Jasminum_tortuosum</td><td>10</td><td>26</td><td>30</td><td>68</td><td>0.05</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>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>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><tr><td></td><td></td><td></td><td></td><td></td><td>0.43</td></tr><tr><td>Podocarpus_totara</td><td>13</td><td>42</td><td>26</td><td>42</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.03</td></tr><tr><td>Porphyra_purpura</td><td>3</td><td>8</td><td>8</td><td>16</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Sagittaria_trifolia</td><td>15</td><td>32</td><td>34</td><td>46</td><td>0.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.02</td></tr><tr><td>SciadopitysVerticillata</td><td>25</td><td>70</td><td>52</td><td>78</td><td>0.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.04</td></tr><tr><td>Sciaphila_densiflora</td><td>4</td><td>8</td><td>8</td><td>8</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Selaginella_kraussiana</td><td>5</td><td>12</td><td>14</td><td>26</td><td>0.01</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.01</td></tr><tr><td>Selaginella_vardei</td><td>3</td><td>8</td><td>8</td><td>16</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.00</td></tr><tr><td>Taxus_baccata</td><td>17</td><td>42</td><td>34</td><td>42</td><td>0.00</td></tr><tr><td></td><td></td><td></td><td></td><td></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></td><td></td><td></td><td></td><td></td><td>0.15</td></tr><tr><td>Uvaria_macrophylla</td><td>13</td><td>30</td><td>42</td><td>86</td><td>0.24</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.84</td></tr><tr><td>Welwitschia_mirabilis</td><td>13</td><td>34</td><td>30</td><td>48</td><td>0.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.03</td></tr><tr><td>Wolffia_australiana</td><td>10</td><td>28</td><td>24</td><td>44</td><td>0.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.02</td></tr></tbody></table>

Dierckxsens N, Mardulyn P, Smits G (2017) NOVOPlasty: De novo assembly of organelle genomes from whole genome data. Nucleic Acids Research 45(4):e18. https://doi.org/10.1093/nar/gkw955
<table><caption>Table D3 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><tr><td>Lathyrus_pubescens</td><td>17</td><td>44</td><td>36</td><td>52</td><td>0.02<br>0.02</td></tr><tr><td>Lophocereus_schottii</td><td>12</td><td>34</td><td>24</td><td>34</td><td>0.00<br>0.00</td></tr><tr><td>Pelargonium_nanum</td><td>11</td><td>32</td><td>36</td><td>128</td><td>0.33<br>0.84</td></tr><tr><td>Podocarpus_totara</td><td>13</td><td>40</td><td>26</td><td>40</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.18</td></tr></tbody></table>

François S, Andonov R, Lavenier D, et al (2018) Global Optimization for Scaffolding and Completing Genome Assemblies. Electronic Notes in Discrete Mathematics 64:185–194. https://doi.org/10.1016/j.endm.2018.01.020

Gurevich A, Saveliev V, Vyahhi N, et al (2013) QUAST: Quality assessment tool for genome assemblies. Bioinformatics 29(8):1072–1075. https://doi.org/10.1093/bioinformatics/btt086

Gusfield D (2019) The RNA-Folding Problem. In: Gusfield D (ed) Integer Linear Programming in Computational and Systems Biology: An Entry-Level Text and Course. Cambridge University Press, Cambridge, p 105–121, https://doi.org/10.1017/9781108377737.008

Huson DH, Reinert K, Myers EW (2002) The greedy path-merging algorithm for contig scaffolding. Journal of the ACM 49(5):603–615. https://doi.org/10.1145/585265.585267

Jin JJ, Yu WB, Yang JB, et al (2020) GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes. Genome Biology 21(1):241. https://doi.org/10.1186/s13059-020-02154-5

Mandric I, Zelikovsky A (2015) ScaffMatch: Scaffolding algorithm based on maximum weight matching. Bioinformatics 31(16):2632–2638. https://doi.org/10.1093/bioinformatics/btv211

Miller CE, Tucker AW, Zemlin RA (1960) Integer Programming Formulation of Traveling Salesman Problems. Journal of the ACM 7(4):326–329. https://doi.org/10.1145/321043.321046

Palmer JD (1983) Chloroplast DNA exists in two orientations. Nature 301(5895):92–93. https://doi.org/10.1038/301092a0

Palmer JD (1985) Comparative Organization of Chloroplast Genomes. Annual Review of Genetics 19(1):325–354. https://doi.org/10.1146/annurev.ge.19.120185.001545

Salmela L, Mäkinen V, Välimäki N, et al (2011) Fast scaffolding with small independent mixed integer programs. Bioinformatics 27(23):3259–3265. https://doi.org/10.1093/bioinformatics/btr562

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

Turmel M, Otis C, Lemieux C (2017) Divergent copies of the large inverted repeat in the chloroplast genomes of ulvophycean green algae. Scientific Reports 7(1):994. https://doi.org/10.1038/s41598-017-01144-1

Wang J, Wong GKS, Ni P, et al (2002) RePS: A Sequence Assembler That Masks Exact Repeats Identified from the Shotgun Data. Genome Research 12(5):824–831. https://doi.org/10.1101/gr.165102

Xiao-Ming Z, Junrui W, Li F, et al (2017) Inferring the evolutionary mechanism of the chloroplast genome size by comparing whole-chloroplast genome sequences in seed plants. Scientific Reports 7(1):1555. https://doi.org/10.1038/s41598-017-01518-5