• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
La molécule d'ADN est le support universel de l'information génomique chez tous les organismes vivants. Étudier in silico l'information qu'elle porte permet de mettre en exergue l'effet de son expression sur le fonctionnement des organismes, de les classifier, de connaitre l'histoire qui lie les êtres vivants entre eux. La molécule d'ADN est double brin. Chaque brin est un long mot formé de quatre lettres, les nucléotides : A, C, G, T . Les deux brins sont complémentaires : un nucléotide A fait face à T , C fait face à G -et vice-versa. Les nucléotides vont par paires, appelées comme paires de bases (pb). De plus, les deux brins sont transcrits dans des sens opposés. Ils sont par conséquent inverses-complémentaires. Pour connaitre sa composition en nucléotides, la molécule d'ADN est fragmentée en plusieurs morceaux chevauchants lors du séquençage. Ces fragments, appelés lectures sont assemblées pour obtenir la séquence nucléotidique d'un seul brin, l'autre brin s'obtenant en inversant et en complémentant la séquence du premier.

Les technologies de séquençage ne permettent toujours pas de séquencer en une seule fois la totalité d'une molécule ADN quand elle excède une certaine taille. À la fin des années 70, les recherches de Sanger permettent de séquencer des fragments dont les tailles des lectures approchent le millier de paires de bases. Ce processus est toutefois long et couteux. Les technologies de séquençage connaissent un développement important lors de la première décennie du 21 e siècle. Des millions de lectures sont produites en une seule exécution de séquençage. Elles sont courtes (de l'ordre de la centaine de paires de bases), et leur taux d'erreur de séquençage atteint au mieux 0,1 %, un ordre de grandeur en dessous de celui des technologies de Sanger. Depuis 2010, une troisième génération de technologies produit des longues lectures (de l'ordre de la dizaine de milliers de paires de bases) mais fortement bruitées (de quelques pourcents).

L'assemblage de génomes ou assemblage de fragments est le processus visant à reconstruire les plus longues sous-séquences des molécules ADN séquencées (sinon leur séquence entière) à partir de leurs lectures. Plusieurs facteurs rendent le problème de l'assemblage de fragments non trivial. Les erreurs de séquençage obligent à la correction au préalable des lectures, sinon à la production d'une séquence consensus à partir des lectures assemblées. Un plus grand nombre de courtes lectures augmente la consommation en mémoire et en temps des algorithmes d'assemblage. Les fragments séquencés proviennent de copies des deux brins. Par conséquent, chaque lecture doit être considérée sous ses deux orientations exclusives :

iii soit sa séquence inchangée participe à l'assemblage (orientation directe), soit sa séquence inverse-complémentaire (orientation inverse), sinon aucune des deux. Le plus difficile étant la résolution des régions répétées dans les génomes. Les occurrences d'une région répétée se distinguent en fonction des régions adjacentes. Plus les lectures sont longues, plus il est probable que les occurrences et leurs régions adjacentes soient couvertes. À l'inverse, si aucune lecture ne couvre entièrement les occurrences et leurs adjacences, leur ordre ne pourra pas être inféré. La résolution des régions répétées est ainsi fonction de la longueur des lectures et des longueurs des régions répétées -par exemple longues de plusieurs milliers de paires de bases dans le génome du mais.

Deux grandes étapes composent l'assemblage de fragments. (i) L'assemblage des lectures à partir de leurs chevauchements (le début d'une lecture s'aligne sur la fin d'une autre). Il en résulte un ensemble de plus longues séquences, les contigs, qui peuvent se retrouver en plusieurs copies dans le génome. (ii) L'échafaudage des contigs (scaffolding en anglais) consistant à les orienter et les ordonner pour former des échafaudages (scaffolds en anglais) grâce à des données supplémentaires. Les contigs d'un échafaudage peuvent être éloignés par une distance nucléotidique.

La thèse présentée ici s'est concentrée sur deux problématiques. La première porte sur l'analyse et la comparaison de structures de graphes utilisées tout au long des étapes de l'assemblage de fragments. Le deuxième consiste en une formulation originale du problème de l'échafaudage dédié aux génomes d'un organisme particulier, les chloroplastes, et à sa résolution.

Graphes de fragments ADN

Un graphe est une structure mathématique utilisée pour rendre compte de relations entre des objets. Ces derniers sont représentés par des sommets reliés par des arêtes (graphe non orienté) ou par des arcs (graphe orienté). L'assemblage des lectures s'appuie sur les chevauchements entre elles. Ainsi, il existe un lien permettant de passer d'une première lecture orientée vers une deuxième. De même, l'échafaudage des contigs dépend des liens entre deux contigs orientés. Ces liens proviennent souvent de paires de lectures séparées par une distance (connue ou non) s'alignant sur deux contigs (paired-end reads en anglais). Pour généraliser, nous utilisons le terme générique fragment. Dans le cas de l'assemblage de lectures, il désigne les lectures. Dans celui de l'échafaudage, il désigne les contigs. Dans ces deux cas, un graphe permet de représenter les liens entre deux fragments orientés.

Trois structures de graphes émergent de la littérature. La première représente chacune des deux orientations de chaque fragment par un sommet. Parce que les liens sont des paires ordonnées de fragments orientés, ils définissent alors des arcs dans un graphe orienté (GO). [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF] l'employa pour la première fois. Les objets y semblent redondants, car les sommets vont par pairs (deux orientations), ainsi que les arcs (les liens peuvent aussi être munis d'une orientation). [START_REF] Myers | Toward simplifying and accurately formulating fragment assembly[END_REF] propose de fusionner les deux orientations d'un fragment en un seul sommet, ainsi qu'un lien et son inverse en une seule arrête. L'arrête doit alors être pondérée par un couple d'orientations. La première correspond à celle du premier fragment, par exemple celui d'identifiant lexicographiquement inférieur. La seconde orientation est celle du fragment défini comme étant le deuxième. Le graphe est biorienté (GB). Enfin, pour leur méthode d'échafaudage, [START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF] simplifient la représentation en associant deux sommets pour chaque fragment, un pour chacune des deux extrémités d'un fragment : la tête et la queue. Ces deux sommets sont reliés par une arête (arrête-fragment). Passez de la tête à la queue correspond à choisir le fragment inverse. Ici, les liens sont des arêtes reliant l'extrémité d'un fragment à un autre (arrêtes-liens). Un chemin valide dans ce graphe non orienté (GN) commence et termine par une arête-fragment et alterne entre les deux types d'arêtes.

L'emploi de l'une ou l'autre structure de graphe varie entre les méthodes pour les étapes d'assemblage de fragments. Bien que leurs différences de conceptions soient parfois décrites, ces structures n'ont jamais été comparées théoriquement ni en pratique. Aucune base d'implémentation pour GO, GB et GN pour comparaison n'a donc été proposée. En réponse, nous avons proposé des implémentations s'appuyant sur des listes d'adjacence qui tirent bénéfice de la symétrie impliquée par l'inversecomplémentarité des deux brins. Les implémentations sont comparées en terme d'usage théorique en mémoire, et au travers du coût temporel des algorithmes pour réaliser des opérations élémentaires : itérer sur les liens, ajouter et supprimer un fragment ou un lien.

Bien que les coûts des implémentations ne diffèrent que de facteurs linéaires, nous concluons sur quelques recommandations. En particulier, nous avons mis en évidence que GB est à privilégier si la mémoire fait défaut, puisque deux fois moins de pointeurs sont requis, ou quand il s'agit de supprimer des fragments du graphe. L'implémentation ne décrivant que les successeurs dans GO est la meilleure pour itérer sur les voisins d'un fragment orienté, ajouter ou supprimer un lien. Une implémentation de GO, qui ne garde que les voisins des fragments orientés directs, est disponible pour Python3 sur PyPI sous le nom de revsymg.

Échafaudage de génomes de chloroplastes

Les chloroplastes sont des organismes qui vivent en endosymbiose avec les plantes et les algues. Ils profitent des mécanismes cellulaires de leur hôte afin de se répliquer et de maintenir leurs fonctions en échange de la production d'énergie à partir de la lumière du soleil via le processus métabolique de la photosynthèse. Les v RÉSUMÉ EN FRANÇAIS Échafaudage de chloroplastes chloroplastes possèdent leur propre génome en plus de celui de la plante (dans le noyau) et de celui des autres organites (comme des mitochondries). Les plantes sont des êtres pluricellulaires, et chacune des cellules de leurs feuilles possèdent plusieurs chloroplastes.

La structure et la répartition du génome des chloroplastes sont particulières. La structure la plus étudiée est un génome circulaire quadriparti, séparée en quatre régions : une Longue Région Unique (LRU), une Courte Région Unique (CRU) et une paire de régions inverse-complémentaires (RI). Chaque chloroplaste possède plusieurs copies de son génome, qui peuvent différer en formes. La présence des RI entraîne l'inversion-complémentaire d'une des régions uniques lors de la réplication qui s'opère sur les deux brins à partir d'un point d'origine. Ainsi, des formes multimériques réversibles coexistent dans un même chloroplaste.

Plusieurs approches ont déjà été proposées pour l'assemblage spécifique d'organites et de chloroplastes. Certaines sont des suites (pipeline en anglais) de méthodes génériques appliquées sur des données de séquençage filtrées [START_REF] Ankenbrand | chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data[END_REF]. D'autres profitent de la faible taille des génomes de chloroplastes et ont développé des algorithmes de graine-et-extension (seed-and-extend en anglais) visant un assemblage circulaire [START_REF] Coissac | From barcodes to genomes: Extending the concept of DNA barcoding[END_REF][START_REF] Dierckxsens | NOVOPlasty: De novo assembly of organelle genomes from whole genome data[END_REF]. [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] estiment la multiplicité des contigs afin qu'elle soit au plus proche des observations des alignements des lectures sur les contigs et que le graphe de contigs multipliés contienne un chemin circulaire. Les chemins circulaires représentant des RI identiques sont marqués comme solutions. Seule cette approche et celle de [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF] (un échafaudeur appliqué sur des données chloroplastiques) tiennent compte de solutions multiples liées aux formes multimériques. Toutes ces approches ne modélisent pas explicitement l'assemblage ou l'échafaudage d'un génome de chloroplaste (et de ses différentes formes), mais procèdent à des filtres a priori et a posteriori des assemblages.

Ici, nous proposons de nous appuyer sur les caractéristiques des structures de leur génome pour formuler l'échafaudage de génome de chloroplaste en un problème d'optimisation combinatoire. Cette formulation se centre sur la reconstruction des régions répétées et des régions uniques, pour ensuite déduire les différentes formes de génomes. Nous démontrons que ce problème est N P-complet. Nous le décomposons en trois sous-formulations (pour les RI, RU, et régions répétées directes), et pour chacune nous proposons un modèle linéaire en nombre entier. Ces modèles sont implémentés dans un unique programme qui les hiérarchise. Nos résultats sont encourageants sur des données synthétiques, choisies en fonction de la variété des difficultés génomiques structurales qu'elles présentent. La méthode est disponible sous le module PyPI Python3 khloraascaf, et les résultats sont reproductibles (voir https://khloraascaf-results.readthedocs.io/en/latest/).

vi

Conclusions et perspectives

L'assemblage de fragments est un problème abstrait difficile à formaliser dans le cas général. Découpé en sous-étapes, leurs approches de résolution dépendent de la nature des données en entrée et des spécificités des génomes des organismes à assembler. Les problèmes d'optimisation combinatoire sous-jacents gagnent en pertinence lorsque leur modélisation rend compte de ces complexités.

Le graphe est une structure mathématique utilisée tout le long du processus d'assemblage de fragments : des lectures jusqu'aux échafaudages. Trois structures émergent de la littérature. Bien qu'elles soient comparées dans leur conception abstraite, elles n'avaient jamais été comparées sur la base d'une implémentation de graphe. Nous proposons dans cette thèse une base commune de leur formalisation en proposant une implémentation sous la forme de listes d'adjacences. Chaque implémentation est mesurée selon sa place en mémoire et du coût temporel des algorithmes d'itération sur les données, d'ajout ou de suppression d'informations. Une des implémentations est disponible via un module Python3.

Nous envisageons d'implémenter les graphes sous Rust ou sous un langage de programmation compilé équivalent pour comparer les implémentations expérimentalement et confirmer les calculs théoriques. De même, puisque le choix d'une structure de graphe influence les modélisations ultérieures pour l'assemblage de fragments, nous souhaitons proposer, pour chaque structure, des modèles linéaires en nombres entiers de recherches de chemins ou de couvertures, les analyser et les comparer en théorie et en pratique.

Nous nous sommes ensuite penchés sur le problème de l'échafaudage dans le contexte des génomes de chloroplastes. La communauté travaille depuis maintenant une vingtaine d'année sur le problème de l'échafaudage dans le cas général, et depuis une dizaine d'années particulièrement sur les génomes de chloroplastes. Les méthodes diffèrent grandement sur les heuristiques utilisées ainsi que sur les stratégies de sous-échantillonnage des données concernant les outils les plus complets. L'intérêt du génome de chloroplaste pour son assemblage réside dans sa structure génomique particulière et dans l'existence de plusieurs formes au sein d'un même organite. Sa petite taille de génome se traduit par un jeu de données de taille relativement petite permettant ainsi de résoudre le problème de l'échafaudage par des méthodes exactes. Ainsi, nous posons une formulation du problème que nous résolvons au moyen d'une décomposition en sous-problèmes, modélisés et résolus par programmation linéaire en nombres entiers. À travers elle, nous couvrons de nombreuses spécificités du génome de chloroplaste, de l'agencement des régions génomiques jusqu'à l'existence de formes multimériques. Notre méthode produit des solutions multiples corrélées avec le phénomène d'inversion-complémentaire de certaines régions au cours de la réplication ADN. Enfin, nous avons implémenté ces vii approches en un module Python3, testé sur des données synthétiques aux difficultés hétérogènes. Les résultats sont reproductibles.

L'échafaudage n'étant qu'une sous-étape du processus d'assemblage, nous souhaitons injecter notre approche dans une suite permettant d'assembler les génomes de chloroplastes à partir des lectures. Dans un souci de comparaison pertinente, nous envisageons d'injecter notre partie dans la suite GetOrganelle, à l'endroit où l'on identifie une équivalence d'objectif, entrées et sorties. Cela nous permettra de nous comparer avec l'état de l'art. Les chloroplastes ne sont pas les seuls organites en endosymbiose avec la plante. Les mitochondries par exemple partagent des caractéristiques génomiques communes avec celles des chloroplastes. Nous envisageons d'investiguer dans quelle mesure notre approche peut-être adaptée pour l'assemblage de génomes mitochondriaux. This chapter presents general and basic knowledge on the subject of genome assembly. Deoxyribonucleic acid (DNA) and genomic information studies have known great development since the 20 th century. They impact many areas, from ecological understanding to economics, through agronomics and medicine.

viii

CONTENTS

We first give a quick overview on what DNA is and the information the molecule holds in Section 1. Section 2 describes the particularities of the chloroplast genome we exploit to formalise their assembly in Chapter IV. Then, we describe sequencing methods to obtain a DNA sequence through its fragmentation (Section 3) and the challenges it raises (Section 4). Section 5 gives the background of the approaches enabling to retrieve the sequence from fragments. Finally, we address the topics this thesis focuses on in Section 6.

DNA overview

DNA molecules are the support of the genomic information contained in each cell of living organisms. Studies in DNA have been developed since the beginning of the 50s, when X-ray pictures shot by Gosling while he was supervised by Franklin [START_REF] Franklin | The structure of sodium thymonucleate fibres. I. The influence of water content[END_REF] lead to the double-strand model proposed by [START_REF] Watson | The Structure of Dna[END_REF] Each strand is a sequence of four nucleotides, adenine (A), cytosine (C), guanine (G) and thymine (T ) that are organic molecules considered as elementary subunits.

We denote by Σ nuc = {A, C, G, T } the nucleotide alphabet, where A and T are complementary nucleotides as well as C and G. Each strand has a reading direction from the defined 5' extremity to the 3' one. One strand is the complement of the other, and the 5' 3' extremities are reversed. As an example, if AAT GCCA is a strand (or a DNA sequence), then T GGCAT T is its reverse-complement as illustrated in Figure 1. For the sake of clarity, the reverse-complement word is shortened to the reverse.

As the DNA molecule is double-stranded, the nucleotides are said to be paired, and the length of a DNA molecule is measured in base pairs (bp). The DNA molecules among the living organisms vary in forms and in length: from millions of base pairs for bacteria [START_REF] Trevors | Genome size in bacteria[END_REF] to billions for eukaryotes [START_REF] Kidwell | Transposable elements and the evolution of genome size in eukaryotes[END_REF]. This double-strand specificity contributes to chemical stability as well as enabling replication and repairing mechanisms. These mechanisms are the key to sharing genomic information through generations.

INTRODUCTION I

Figure 1 -DNA double-strand molecule 5' and 3' extremities give the reading direction of each DNA strand, illustrated by the two arrows. One strand is the reverse-complement (or the reverse) of the other. Thus, AAT GCCA is the reverse of T GGCAT T , and vice-versa. This DNA molecule's length is 7 bp.

A recipe for cellular mechanisms

DNA operates as both an archive and a recipe for cellular machinery, through hierarchical processes of synthesis. The DNA molecule can be divided into gene regions separated by intergenic regions. Particularly, enzymes read the DNA molecule from the 5' extremity to the 3' one, and synthesise RNA single-strand molecules from the genes (transcription stage). RNA molecules are then translated to proteins that can form protein-complexes performing different functions such as the catalytic function of enzymes. Especially during the translation stage of the RNA (denoted mRNA, standing for "messenger RNA") the nucleotide alphabet is translated to the amino-acid alphabet, where three consecutive nucleotides correspond to one amino-acid. Specific molecular biology research fields focus on each molecule and the synthesis processes. Various types of RNA cover different functions: messenger, participating in the translation stage, silencing genes etc.

On the motivations for knowing the whole DNA sequence

Uniquely studying proteins is insufficient to understand all the cellular mechanisms and their influences on phenotypes. Studying how DNA transitions into mRNA or tRNA completes knowledge of the DNA's 3D behaviour, helping researchers to understand life processes. Similarly, gene sequence analysis alone cannot fully explain the functioning of living organisms. Whole genome sequencing provides supplementary results compared to protein, gene or RNA studies. The interested reader is referred to the review by [START_REF] Rice | New Approaches for Genome Assembly and Scaffolding[END_REF] from which the following derives.

Without specifically studying the function of the genomic regions, a whole genome can serve as a reference sequence on which other sequences are aligned or mapped in order to infer history, timing or location of events based on edition distances. Gene expression (reading a gene results in the production of a mRNA) can be affected over large genomic regions: a gene can be silenced or over-expressed due to inhibitor/enhancer genomic regions. Indeed, while it is expected that gene expression is affected by genomic regions near to the gene in the DNA sequence, these regions can also be near in space because of DNA molecular folding mechanisms. Some phenotypes are multifactorial corresponding to polygenic traits (e.g. gene co-expression), and it can be valuable to know the location order of the gene of interest. Also related to gene expressions, during the meiosis stage (cell division of germ cells for the production of gametes) chromosome recombination sometimes occur. Finally, changes in chromosomes can result in the emergence of new species.

Therefore, having the complete genome sequence influences biochemistry, immunology, evolution and ecological sciences.

Various molecular conformations structure the genome

The DNA or genomic information is differently organised across living organisms. The genomic information can be organised along one molecule or split into several ones. It can also be multiplied in the genomic container (e.g. in the cell's nucleus for eukaryotes). When multiple copies coexist in the container, they can be grouped: they are paired for humans (diploid species), or grouped by more than two, e.g. by 6-uplet as for the wheat (hexaploid). The DNA molecules supporting the whole or a part of the genomic information also differ in topology: they can be linear (as for humans) or circular (as for the Escherichia coli bacteria). Figure 2 illustrates divers DNA molecule conformations.

Chloroplast genome

Chloroplasts are plant and algae cell organelles. They are the result of endosymbiosis between bacteria and eukaryotic cells [START_REF] Gould | Algae's complex origins[END_REF]. They produce carbon compounds from CO2 in the air through the photosynthesis biological process. Mitochondria are other organelles in endosymbiosis with plants and algae. In

INTRODUCTION I

1986, the first chloroplast genomes were sequenced and assembled. Those are chloroplasts of a liverwort, Marchantia polymorpha [START_REF] Ohyama | Chloroplast gene organization deduced from complete sequence of liverwort Marchantia polymorpha chloroplast DNA[END_REF], and of a tobacco plant, Nicotiana tabacum [START_REF] Shinozaki | The complete nucleotide sequence of the tobacco chloroplast genome: Its gene organization and expression[END_REF]. Knowing the sequence of organelle genomes allows performing ecological, phylogenetic and evolutionary studies [START_REF] Daniell | Chloroplast genomes: Diversity, evolution, and applications in genetic engineering[END_REF][START_REF] Sun | Evolutionary and phylogenetic aspects of the chloroplast genome of Chaenomeles species[END_REF][START_REF] Long | Complete chloroplast genomes and comparative analysis of Ligustrum species[END_REF].

Genome division

Plants are multicellular organisms. As eukaryotes, the genomic information of the plants is hosted in the nucleus of their cells, noted nucDNA. Chloroplasts and mitochondria, possess their own genomic material (respectively cpDNA and mtDNA), and each plant's cell can host several organelles possessing the same genomic material. Furthermore, each chloroplast has multiple copies of its genome [START_REF] Bendich | Why do chloroplasts and mitochondria contain so many copies of their genome?[END_REF][START_REF] Kobayashi | Detection and Localization of a Chloroplast-Encoded HU-Like Protein That Organizes Chloroplast Nucleoids[END_REF]. The number of copies have been shown to decrease with the leaf age [START_REF] Kumar | Changes in DNA damage, molecular integrity, and copy number for plastid DNA and mitochondrial DNA during maize development[END_REF]. Figure 3 The grey octagon container represents a plant's cell. The purple, blue and green oval shapes are respectively the plant cell's nucleus, mitochondria and chloroplasts. There are several mitochondria and chloroplasts in a same cell. Each of them possesses its own genomic material: in the nucleus, for this artificial example, the plant is diploid (2x = 4); each mitochondrion has several copies of its genome; each chloroplast has multiple forms of its own genome, linear, circular and multigenomes.

Genome multimeric forms

The molecular form of the copies of chloroplast genomes differ [START_REF] Bendich | Circular chloroplast chromosomes: The grand illusion[END_REF], but the circular form has been mainly studied. For example [START_REF] Seyer | A physical map of Nicotiana tabacum plastid DNA including the location of structural genes for ribosomal RNAs and the large subunit of ribulose bisphosphate carboxylase/oxygenase[END_REF] give an electron micrograph of a circular cpDNA from Nicotiana tabacum. Although the circular structures are a minority, and multigenomes can be found (linear genomes combined in a tree-like structure due to recombination-dependent replication), the circularity property is less sensitive to virus interactions and thus is dominantly transmitted to the next generation [START_REF] Bendich | Circular chloroplast chromosomes: The grand illusion[END_REF].

For this reason, in this manuscript we focus on the chloroplast circular genome copies. During the DNA molecular replication of one cpDNA molecule, flipflop inversions can occur so that the resulting replicated molecules' form differs. This transformation is due to the presence of inverted repeats, where the DNA subsequence between them is reversed. Inversions during the replication phase were also studied for bacteria, and they are considered as reversible operations. The copies of a genome, exact or obtained by inversions, are denoted as multimeric forms [START_REF] Palmer | Chloroplast DNA exists in two orientations[END_REF][START_REF] Deng | The chloroplast genome exists in multimeric forms[END_REF][START_REF] Wang | Long-Reads Reveal That the Chloroplast Genome Exists in Two Distinct Versions in Most Plants[END_REF]. Figure 4 schematises a flip-flop inversion for a chloroplast genome with inverted repeats. During the DNA replication of the chloroplast genome, one of the region between the inverted repeats can be reversed (here the red region SSC). This provokes the existence of several forms of the genome in the same chloroplast characterised as multimeric forms.

Flip-flop inversion

Repeats and single-copy regions

Genomic regions can be classified according to whether their sequence is found approximately elsewhere in the genome, regardless of their functions. In this case, the classification depends on two main factors:

-the minimum length of a sequence to consider it as a repeat (otherwise, a sequence of size one could be said to be repeated);

INTRODUCTION I

-an acceptable number of word editions (such as substitutions, insertions and deletions) to pass from one sequence to another.

The difficulty with these factors lies in the fact that they are left to the discretion of the researchers analysing them, and therefore suffer from a certain subjectivity.

For the chloroplast we will focus on two types of repeats.

Direct Repeat (DR)

The sequences are highly similar;

Inverted Repeat (IR) One sequence is the reverse of the other.

Figure 5 gives two examples of perfect repeats. One of the most studied chloroplast genome architectures is the circular DNA molecule whose quadripartite sequence includes a pair of highly similar (or identical) reversed nucleotide subsequences (IR) separated by a Long Single-Copy (LSC) and a Short Single-Copy (SSC) sequences [START_REF] Bock | Genomics of Chloroplasts and Mitochondria[END_REF]. Figure 6a illustrates IR-quadripartite chloroplast genome. For example, [START_REF] Thode | Comparative analyses of Mikania (Asteraceae: Eupatorieae) plastomes and impact of data partitioning and inference methods on phylogenetic relationships[END_REF] show that IRs in Mikania plastomes have length around 25 kbp. However, an example of the difficulties in discriminating what is a repeat is highlighted by [START_REF] Kim | Widespread Occurrence of Small Inversions in the Chloroplast Genomes of Land Plants[END_REF]. The latter studies quadripartite chloroplast genomes described as above, where inversions can be found in LSC because of small IRs. The structures are still described as quadripartite because the IRs in the LSC are smaller than the IRs considered in the structure. Quadripartite structures with DR, generally shorter than the IRs, can also be found [START_REF] Palmer | Comparative Organization of Chloroplast Genomes[END_REF], as well as more complex structure involving both type of repeats [START_REF] Tsai | Dispersed repetitive sequences in the chloroplast genome of Douglas-fir[END_REF]. Figure 6 summarises the main chloroplast genome structures found in the literature.

Genome evolution

Chloroplasts and mitochondria in endosymbiosis with eukaryotic cells have gradually become the organelles of plant cells. As a result, organelles no longer need specific genomic traits, leading to a loss of genes and a reduction of their genome size [START_REF] Xiao-Ming | Inferring the evolutionary mechanism of the chloroplast genome size by comparing wholechloroplast genome sequences in seed plants[END_REF]. Horizontal transfers (a transfer of genes not by heredity) were also studied from the plant to its organelles.

DNA sequencing

DNA sequencing aims to obtain the nucleotide sequence of a DNA molecule for downstream analyses. The sequence of only one DNA strand is sufficient because retrieving its complementary is immediate (c.f. Figure 1). Since the 1960s, it has been necessary to sequence subpart of molecules, corresponding to a few nucleotides. Although fragment lengths in 2023 could exceed the megabase range, sequencing the entire DNA molecule of most living organisms in one run remains unfeasible.

INTRODUCTION I

In some papers presenting sequencing results, the term sequencing refers both to the genome fragmentation and the assembly of the fragments. Here we separate the two procedures. Sequencing is the process of fragmenting DNA and obtaining their sequences using a sequencer. The fragment assembly solves the fragment puzzle with computational methods.

In the late 1960s, the available techniques were restricted to employing digestion enzymes which split the sequence at specific base positions. As a result, they output very short fragments, denoted as digests, sometimes of length two and generally limited to length eight. Sequencing DNA molecules of several kilobases was therefore impossible, and the digest sequencing technology was limited to sequence protein or RNA sequences.

The Shotgun sequencing approach

Since the 1970s, sequencing technologies have known continuous development and fast improvement leading to various data types. In the following, a read refers to a sequence obtained by these technologies, while a fragment refers to a partial or a whole sequence of a DNA molecule. The main approach to sequence a genome is proceeding to a shotgun sequencing: the reads are obtained from the sequencing of fragments randomly sampled from copies of the genome. This strategy was firstly proposed by [START_REF] Staden | A strategy of DNA sequencing employing computer programs[END_REF] and is illustrated in Figure 7. It is based on the fact that a fragment from a copy of the genome may overlap fragments from other copies. Since the reads are subsequences of the fragments, they also overlap, enabling to retrieve a continuous sequence.

The next sections summarise the main sequencing technologies, inspired by the reviews written in [START_REF] Mehdi | Generations of Sequencing Technologies: From First to Next Generation[END_REF] and [START_REF] Pervez | A Comprehensive Review of Performance of Next-Generation Sequencing Platforms[END_REF].

First generation: Sanger and BAC technologies

In the late 1970s Sanger's research succeeded in sequencing much longer fragments [START_REF] Sanger | DNA sequencing with chain-terminating inhibitors[END_REF]. Their length increases from 8 bp up to one thousand (400 bp to 900 bp). The process is denoted by Sequencing By Synthesis (SBS): it consists in completing a single-strand DNA molecules from a fragment (a template) with denatured nucleotides. Each nucleotide addition is marked along a gel thanks to electrophoresis. The gel can be seen as an n × 4 boolean matrix, where n is the fragment's length and each of the four columns corresponds to one nucleotide. If there is a band at the cell (i, j), then there is the j th nucleotide at the i th position of the fragment.

The Bacterial Artificial Chromosome (BAC) benefits from a bacterial DNA molecule as a vector to amplify inserted fragments of interest (O'Connor et al.,

Coverage

Duplication

+ Fragmentation

Fragment sequencing From the bottom to the top of the figure: the reads (in blue) are sequenced from longer fragment (in green) randomly sampled from copies of the genome (in orange). Based on the overlapping reads, the assembly process aims to retrieve the genome. The coverage counts the average number of times a base in the genome is covered by the reads.

1989). [START_REF] Venter | A new strategy for genome sequencing[END_REF] suggested sequencing them at both ends. Thanks to adaptors at both extremities of the fragments, two reads are sequenced from the 5' extremity of the two complementary strands (Denoted as paired-end data, in opposition to single-end, c.f. Figure 9). As the length of the vector is known, BAC technology provides an approximated distance between the two paired-read.

Although the fragments are long and the accuracy can reach 99.99 %, the Sanger technology suffers from low throughput and high costs. Since the beginning of the 21 th century it has been supplanted by high-throughput technologies, even if it is still used for validating DNA sequence and for target resequencing.

Second generation: high-throughput sequencing

As mentioned above, the first decade of the 21 st century witnessed the emergence of new sequencing technologies. Refered to as Next Generation Sequencing technologies (NGS), they produce millions of short reads at low prices allowing the sequencing of numerous organisms. The first error type is someone noted as a mismatch, and the two latters are categorised as indels.

Illumina sequencing (2000) The Illumina/Solexa sequencing technology is also an SBS technology. Adaptors are attached to the extremities of the double-strand fragments. The two strands are separated, the adaptors find their complementary on a flow cell, and the single-strand fragments are then amplified by Polymerase Chain Reaction (PCR). As modified nucleotides complete the strands, a laser emits a light signal that is detected by a camera and interpreted by a computer. The best machines are able to sequence millions of reads at a reduced cost, by providing a 0.1 % to 1 % error rate in favour of substitutions (Figure 8a). However, the read lengths are in the range of 100 bp to 300 bp.

sequencing (2005)

Similarly to Illumina technology, the 454 sequencing is labelled SBS, it needs PCR amplification and nucleotide additions are detected by light emission. It produces reads with lengths from 100 bp up to 700 bp with a 1 % error rate mainly caused by homopolymers errors. Homopolymers are repetitions of the same nucleotide, and their misdetection leads to insertions and deletions (indels) in the reads (Figures 8b and8c).

Ion Torrent (2011) SBS technology, it measures the PH variation to detect nucleotides. It produces reads with lengths from 200 bp up to 400 bp, with 1 % error type mainly due to indels.

Paired-end data Analogously to the BAC technology, it is possible to sequence short reads by pairs. For Illumina, the distances are short (100 bp -300 bp), while for mate-pair data (circularised fragments generating paired reads, [START_REF] Korbel | Paired-End Mapping Reveals Extensive Structural Variation in the Human Genome[END_REF], they reach the kilobases at the expense of accuracy.

Third generation: single molecule run time

The Third Generation Sequencing technologies (TGS) refer to Single Molecule Run Time technologies (SMRT) that have emerged in the last decades. They produce fewer but longer reads without PCR. The procedure is shorter and less costly. Although they first suffered from relatively high error rate with indel error types, these technologies have made significant advances in lowering the error rate in recent years.

Pacific biosciences SMRT sequencing (2010) Each DNA fragment passes into a cell that contains a DNA polymerase. When the single-strand fragment passes through the polymerase, fluorescent labelled nucleotides complete the strand and release light signals. The average read length is 10 kbp and can reach 60 kbp. The first error rate was measured to 13 % but more recent HiFi technology would decrease the rate to 0.5 % (according to PacBio, [START_REF] Hon | Highly accurate long-read HiFi sequencing data for five complex genomes[END_REF].

Oxford nanopore sequencing (2014) In Oxford Nanopore Technology (ONT), the DNA fragments pass though a pore (made by a protein complex). The passage of the molecule through the pore generates a variation in the ionic current. The variations are recorded and translated to the nucleotide alphabet. The average read length is 10 kbp and can reach 150 kbp, and even 4 Mbp for PromethION machine. Because controlling the speed at which the fragment enters the pore is challenging, the initial error rates are approximately 12 %. According to the firm, they can obtain reads with less than 1 % error rate. However, [START_REF] Delahaye | Sequencing DNA with nanopores: Troubles and biases[END_REF] measure it to 5 % -8 % but remark that the base-caller (the software interpreting ionic signals) is upgrading fast.

INTRODUCTION I

Supplementary sequencing data

The following new data type can help to assemble, finish or correct an assembly of the reads, but at a supplementary cost.

10x linked read 10x linked read is a sequencing technology permitting to label the reads that come from the same fragment. It generates short reads that can be assembled by groups.

Hi-C protocol Hi-C is a technology highlighting physical contact between two regions in the genomes. It generates Illumina paired reads coming from two close genomic regions due to chromosome conformation. While two close regions during one conformation can be far away from each other along the genome, it is possible to prevent the compaction of the DNA in such a way that contacts only occur between two regions close to each other in the sequence. Such contacts are not very accurate and provide a square matrix (each side the length of the genome) with a strong diagonal and a lot of noise around.

Optical map Finally, optical mapping is a process of partitioning copies of the genome at specific known subsequence site thanks to nicking enzyme. For each copy a specific enzyme is associated. As a result, the positions of several subsequences are known, which helps for mapping and ordering fragments along the genome.

The challenges of fragment assembly

The fragment assembly problem began in the late 1960s with the assembly of digests. As the fragments are just a few base lengths (up to 8 bp) it raises the issue of the unicity of their assembly. [START_REF] Shapiro | An Algorithm for Reconstructing Protein and RNA Sequences[END_REF] proposes an algorithm for RNA and protein assembly based on the digest overlap, and gives sufficient conditions for the unicity. Since 1977, more fragments are produced thanks to the Sanger sequencing technology, and has become necessary to use computers to assemble them. The first assemblers were proposed by [START_REF] Gingeras | Computer programs for the assembly of DNA sequences[END_REF] and [START_REF] Staden | A new computer method for the storage and manipulation of DNA gel reading data[END_REF]. They enable the first shotgun assembly project with the assembly of the 50 kbp virus λ using reads of length 200 bp [START_REF] Sanger | Nucleotide sequence of bacteriophage λ DNA[END_REF].

Read sequence alignment

Read sequences can be compared by sequence alignment algorithms. Figure 10 schematises three types of read alignment. (a) Overlapping sequences in different reads enable the algorithm to reconstruct the original DNA sequence. However, because of sequencing errors, two reads sequenced at the same genome region can differ in sequence, so it is necessary to find approximated overlaps. (b) Reads that are contained in others are useful to correct the sequencing errors by producing a consensus from the alignments. (c) Reads sharing a part of their sequence, and whose alignment does not fall into the previous categories, can be considered as detecting repeats because the flanking regions are not aligned.

Unknown fragment orientations

As the DNA is double-stranded, and the reads are sequenced from DNA fragments, two reads can be sequenced from the two different strands. This introduces the notion of sequence orientations. 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). As a consequence, (i) both the direct and its reverse must be considered when searching for read overlaps and (ii) at most one orientation must be chosen for the reads when assembling them. Figure 11 illustrates the fragment orientation issue.

Sequence similarity: single-copy or repeat?

The main challenge of genome assembly is the genome repeats resolution. The definition of a repeat depends on the length of the reads and the similarity between the repeats in the genome. Indeed, if reads are longer than the repeats, there are some reads that contains them with their non identical flanking sequences. Note that some repeats can exceed kilobases (e.g. long transposable elements in the

INTRODUCTION I Fragment 1
Fragment 2

Figure 11 -Unknown fragment's orientation. R aw is the set of sequenced reads. Reads u f and v f have been sequenced from two different fragments of the same genome, but they do not come from the same strand. u r and v r are respectively the reverse of u f and v f . An assembly should either consider u f followed by v r , or v f followed by u r . maize genome or in bacteria genomes [START_REF] Kidwell | Transposable elements and the evolution of genome size in eukaryotes[END_REF]. Raising the threshold to accept overlaps is not sufficient as some repeats are very similar in sequence, and because of sequencing error the alignments must be flexible such as not losing true overlaps. Figure 12 illustrates the genome repeats impact on the reads' assembly in the case of two direct repeats. Furthermore, the overlapping of the reads must also take into account inverted repeats because both orientations must be considered.

True sequence divergences or sequencing errors?

As previously mentioned, the repeats can slightly differ in sequence. Sequencing errors can hide the signal of a true difference. Although a substitution can be resolved by consensus, indel errors are more difficult to solve, and have the worse consequence as they can shift the protein reading phase (three nucleotides are translated into one amino acid). It is thus necessary to have a sufficient sequencing depth to counterbalance the errors. These errors are supposed to be randomly distributed along the genome (which is not really the case, e.g. for ONT). Identifying true biological signals is a prior step to further separating heterozygotes regions for polyploid organisms or splitting close ones shared by different species in metagenomic samples. 5 Fragment assembly approaches 5.1 The Shortest Common Superstring [START_REF] Peltola | SEQAID: A DNA sequence assembling program based on a mathematical model[END_REF] tackle the problem of fragment assembly with the Shortest Common Superstring (SCS) algorithm. Given an error ratio δ, each read must participate in the superstring F , in only one orientation, modified in sequence with an edit distance less than δ times the length of the read. F must be as short as possible, and this makes the problem mathematically non-trivial, although it is not necessarily a biologically meaningful requirement (parsimony principle).

Assembling the fragments by choosing the longest overlaps first gives a heuristic to solve the SCS but does not guarantee the optimal solution. [START_REF] Räihä | The Shortest Common Supersequence Problem over Binary Alphabet is NP-Complete[END_REF] have previously modelled SCS as the search for a Hamiltonian path in a directed graph (digraph) that represents overlaps. In that case, the fragment

INTRODUCTION I

assembly problem is N P-hard. [START_REF] Li | Towards a DNA sequencing theory (learning a string)[END_REF] provided hypotheses that support the SCS modelling: the DNA molecule is assumed to be a random nucleotide sequence and the reads are uniformly sampled. However, these hypotheses do not reflect genomic complexity. Finally, [START_REF] Myers | Toward simplifying and accurately formulating fragment assembly[END_REF] found that the SCS can wrongly make the genome shorter by merging the repeats.

Overlap-Layout-Consensus

To address the sequencing errors and the repeat issues, Kececioglu and Myers introduced the Overlap-Layout-Consensus (OLC) approach. It consists of three stages: (i) computing approximated overlaps between the reads; (ii) finding a layout on the overlaps to partition them into parts classified as single-copy and repeated regions; (iii) producing consensus sequences for each part.

The overlaps computed in stage (i) can be represented in a digraph denoted as the overlap graph. [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF] defined it as G = (V, E), where V is the set of oriented reads and E is the set of arcs (u, v), u and v are two overlapping oriented reads (Figure 13b).

The fragment orientation and the maximum-weight dovetail-chain branching problems

The premise formulations are given in [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF] and [START_REF] Kececioglu | Combinatorial algorithms for DNA sequence assembly[END_REF]. At this time, the orientation of the reads and the ordering of the oriented reads are separated in two global optimisation formulations: first, choose a subset of the overlaps that maximises the sum of their lengths and assign only one orientation for each read; then, find a non-branching path that maximises the overlap length. The authors observed that with known orientations and without errors, the latter corresponds to the SCS.

The string graph formulation To overcome the SCS shortened repeats, [START_REF] Myers | Toward simplifying and accurately formulating fragment assembly[END_REF] proposed the maximum likelihood -valid layout formulation. It consists in preserving the read coverage along the assembly. Indeed, the coverage rate of compressed repeat regions should be abnormally higher than that of single-copy regions in the assembly. This leads to the string graph layout: the overlap graph is transformed to a string graph. (i) The contained reads are removed from V (criticised in Jain 2023); (ii) the graph is simplified by transitive edge removal (in pseudo-linear time complexity, Myers 2005); (iii) the single paths are collapsed. It results in a string graph for which, theoretically, edges correspond to genomic regions and the vertices are the junctions between them. Each edge is weighted by the number of times it may participate in a walk representing the genome. Celera

I INTRODUCTION

Fragment assembly approaches 5 implemented this approach and enabled to assemble the genome of Drosophila melanogaster [START_REF] Myers | A whole-genome assembly of Drosophila[END_REF].

Local and variant approaches The string graph formulation has inherited many variants. Some of them are focused on data cleaning strategies and successive local corrections, as introduced by ARACHNE [START_REF] Batzoglou | ARACHNE: A whole-genome shotgun assembler[END_REF] for BAC data, and done in Canu (Koren et al. 2017, Celera's successor) for TGS data. Some adopt aggressive strategies such as the Best Overlap Graph (BOG), which keeps only the highest-quality extending overlap for each orientation of each read. Introduced in CABOG, a software dedicated to Roche 454 and Illumina paired-end reads [START_REF] Miller | Aggressive assembly of pyrosequencing reads with mates[END_REF], the BOG strategy is part of hifiasm [START_REF] Cheng | Haplotyperesolved de novo assembly using phased assembly graphs with hifiasm[END_REF] which is dedicated to haplotype resolving thanks to HiFi long reads.

OLC complexity In [START_REF] Myers | The fragment assembly string graph[END_REF], a valid assembly of the genome is equivalent to finding a cyclic walk respecting the number of times it passes through each edge in the string graph. According to the parsimony principle, this walk minimises the genome length the assembly produces. In this case, [START_REF] Medvedev | Computability of Models for Sequence Assembly[END_REF] proved the problem to be N P-hard. Nevertheless, due to the plethoras of heuristics in the current methods (e.g. in Canu and hifiasm), the actual complexity of the algorithms is highly correlated to the computation of the overlaps, a quadratic procedure.

De Bruijn Graph approach

New approaches have emerged to handle the millions of reads produced by NGS. Sequencing-by-Hybridisation (SBH) was the name of a left-behind sequencing technology that indicated the presence of sequences of size k among a set in the genome. It inspired the De Bruijn Graph (DBG) assembly approach. In 1989, Pevzner introduced the use of a partial de Bruijn graph (dbgraph) to assemble the genome. This approach first decomposes each read in words of length k (k-mers) overlapping by k -1 nucleotides. Then it builds a dbgraph G = (V, E) of all the k-mers of all the reads is built: each vertex in V corresponds to a k -1 overlap between two k-mers, and each edge in E is a k-mer. Each read corresponds to a walk in G, and the genome corresponds to a superwalk (a walk of walks) (Figure 13c).

As for the OLC approach, it consists in reducing G: (i) Each path of vertices that do not have more than one in-edge and one out-edge are collapsed in one edge (unitig); (ii) branching vertices are resolved with read mapping on the paths. [START_REF] Idury | A New Algorithm for DNA Sequence Assembly[END_REF] For the example, reads 3 and 7 are exactly the same and are contained in R. (b) From the overlap of the read (at least two bases), the string graph is built. Note that neither read 3 nor read 7 are removed because no one is longer than the other while they are mutually contained each other. (c) The reads are split into k -1 overlapping words of length k (k-mers). In the (partial) edge-centric dbgraph, vertices are the k -1 exact overlaps and edges are the k-mers. Here k = 3. For (b) and (c): unique paths are collapsed to form unitigs. The over-simplified illustrations for the OLC and respectively the DBG approaches suggest that the genome can be found by solving a Hamiltonian path through the simplified string graph resp. an Eulerian path through the simplified dbgraph.

than the search for a Hamiltonian path. Furthermore, it reduces the consensus phase to find paths through the k-mers.

Size of k and solid k-mers However, as the overlaps between the k-mers are exact, G is sensible to sequencing errors. In fact, each substitution in a k-mer forms a bubble (two paths starting from and ending to the same k-mers). For a small value of k, the graph is more connected but more branching, while for a larger value the graph branches out less but suffers from bubbles and can be partitioned into several connected components. One strategy is to compute the number of occurrences a k-mer appears in the set of reads, and only use those associated with a sufficiently large occurrence value to build the graph (solid k-mer, [START_REF] Pevzner | An Eulerian path approach to DNA fragment assembly[END_REF].

Local and variant approaches Many DBG variant algorithms have implemented their own heuristic for specific targets. Velvet and SOPAdenovo merge bubbles thanks to a Dijkstra-like Breadth-First Search (BFS) algorithm [START_REF] Zerbino | Velvet: Algorithms for de novo short read assembly using de Bruijn graphs[END_REF][START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF]. SPAdes applies the DBG approach with different values of k and iteratively merges the results [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to singlecell sequencing[END_REF]. While the latter software are suitable for short paired-end reads, the DBG principle inspires methods for long erroneous reads. Flye first generates a draft assembly without considering misassemblies, and detects repeats by aligning the draft against itself [START_REF] Kolmogorov | Assembly of long, error-prone reads using repeat graphs[END_REF]. While wtdbg2 follows the OLC approach, it does not compute the overlap precisely: it partitions the reads in K groups (bin) of 256 bp (K-bin) so that two K-bin are the same if they share enough k-mer [START_REF] Ruan | Fast and accurate long-read assembly with wtdbg2[END_REF]. Finally, mdBG, for Minimiser dbgraph, represents the reads as sequence of minimiser (in a minimiser alphabet), that are a subset of all possible k-mers [START_REF] Ekim | Minimizer-space de Bruijn graphs: Whole-genome assembly of long reads in minutes on a personal computer[END_REF]. [START_REF] Pevzner | L-Tuple DNA Sequencing: Computer Analysis[END_REF] has proposed two dbgraphs: the one described in this section is edge-centric (the k-mers are the edges), while the other one is vertex-centric (the k-mers are the vertices). Pevzner abandoned the latter because he associated it with the Hamiltonian path problem that is N P-complete and for which no algorithm had been proposed. On the other hand, finding an Eulerian path can be done in linear time on the number of vertices and edges under the Euler hypothesis on even degrees [START_REF] Hierholzer | Ueber die Möglichkeit, einen Linienzug ohne Wiederholung und ohne Unterbrechung zu umfahren[END_REF]. However, DBG corresponds to finding a super-walk. [START_REF] Medvedev | Maximum Likelihood Genome Assembly[END_REF] associate DBG with the Chinese postman problem (that echoes the SCS approach) on bidirected graphs, and propose an exact polynomial time algorithm. They next link DBG to a maximum-likelihood problem to avoid compressing the repeat at the opposite of the SCS, similar to the Myers' propositions concerning the near-constant coverage of the reads on the assembly. In the latter case, they prove DBG to be in P.

DBG complexity

INTRODUCTION I

Breaking down the fragment assembly problem

Because of sequencing errors and the difficulties of defining repeats in OLC or dbgraphs, formulating the fragment assembly problem as a global optimisation problem sounds unreachable. The literature suggests splitting the fragment assembly problem into several stages. The development of stage-dedicated methods and software has enabled researchers to create their pipelines depending on the available sequencing data. Here we present a list of hierarchical stages:

Filtering and correcting the sequencing data Reads associated with a low sequencing quality score can be partially or entirely removed from the read set. Multiple alignment techniques can also produce a self-correction of the reads before assembling them.

Assembling the reads into contigs Contigs are nucleotide sequences longer than the read. While unitigs correspond to unambiguous paths in the string and the dbgraphs, contigs result from resolving branches thanks to read alignment or paired-end data.

Scaffolding the contigs Scaffolds are sequences of oriented contigs with potential gaps between them. In Chapter II, Section 2, we describe several approaches for orienting and ordering the contigs.

Gap-filling the scaffolds The gaps between the oriented contig in each scaffold are filled thanks to overlapping reads joining the two contigs' extremities.

Haplotyping Generating consensus sequence on unitig in the OLC approach does not handle heterozygosity for polyploid organisms. Similarly, it does not address the difference between near-genome regions shared between two organisms in metagenomic samples. Read alignment can differentiate a sequencing error from a true signal. Finally, each group is assigned to homologous chromosomes or specific organisms (phasing). Concerning the DBG approach, the separation can use k-mer counting between the paths in bubbles.

Generating a consensus sequence Multiple alignment is the main approach to produce consensus sequences. The final nucleotide bases result from a weighted vote.

Evaluating the assembly The assembly evaluations focus on the number of contigs, their length distribution and their cumulative length. A reference genome enables mapping the contigs to measure mismatch and indel mapping error rates. It may be necessary to split the contigs in order to compare them to the reference: the number of breakpoints denotes the number of misassemblies. One can also compute the number of genes present in the assembly.

Addressed research topics

In this thesis we focus mainly on two aspects of genome assembly.

Graph structure for read assembly and scaffolding stages

Chapter II, Section 1, and Chapter III concern graph data structure for storing and iterating over links between oriented DNA fragments. In the read assembly stage, the links are the overlaps between the reads in the OLC approach (overlap graph). In the scaffolding stage, they often correspond to paired-end information between two reads, translated to link between oriented contigs (scaffolding graph). The overlap graph or the scaffolding graph are fundamental data structure in read assembly and scaffolding stages. Furthermore, their descriptions and their representations influence methods.

In the literature, the graph structures are described mathematically in the best cases, very often commented and illustrated, and sometimes suffer from confused description. To the best of our knowledge, no one has compared the impact of the different representations and especially their implementations.

In this thesis we analyse different representations and propose suitable implementations that we organise in a graph and implementation map. We then design associated algorithms to iterate over the vertices and the edges, and to dynamically maintain them.

Scaffolding of multimeric chloroplast genome forms

Chapter II, Sections 2 and 3, and Chapter IV focus on the assembly and the scaffolding of chloroplast genomes. We first summarise the scaffolding problem formulations and strategies in the general case, mixed with the read assembly stage or as independent stages. We then describe the whole fragment assembly process approaches dedicated to chloroplast genome from the literature.

As we consider that the dedicated approach does not appropriately handle the particularities of the chloroplast genome structure, we translate the biological knowledge into mathematical properties. We then formulate the scaffolding problem in the case of chloroplasts and model it as a global optimisation problem. We prove the decision version of the problem to be N P-complete. Furthermore, our formulation enables us to tackle the presence of multimeric genome forms in chloroplasts. Finally, we have implemented our approach and have tested it with synthetic data in order to measure the time complexity in practice and the robustness. Succession relationships between oriented reads are the keys to assembling them into contigs. Those between oriented contigs are necessary to assemble them into scaffolds. Those between oriented scaffolds enable a chromosome-scale assembly.

II STATE-OF-THE-ART

All over the fragment assembly process, from the read assembly to the gap-filling through the scaffolding stages, it is fundamental to use suitable data structures to store and query these succession relationships between oriented fragments. Most of the OLC approaches for the read assembly stage, and most of the scaffolding methods use graph structures to address the above issues. In this section, we present the underlying graph structures. We denote the graphs representing the fragments by abstraction of their sequence as fragment graphs, at the opposite of the dbgraph for which the vertices and the edges represent words. For example, two fragments with identical sequences may have several vertices in the fragment graph, while they will be compressed in the same path in the dbgraph.

Notations and fundamental definitions

Fragment set

A DNA fragment can either be a read resulting from the sequencing, a contig from the assembly of the reads, or a scaffold from the scaffolding of the contigs. The below definition serves to rigorously define an oriented fragment afterwards.

Definition 1.1: Unoriented fragment set

Denote by F un the set of unoriented fragments, in bijection with the set of integer label Σ Fun = 0, |F un | , and denote by unf id : F un → → Σ Fun the labelling function. Let Σ nuc = {A, T, G, C} be the nucleotide alphabet and unseq : F un → Σ nuc + the function that associates a fragment with its nucleotide sequence.

Because of the double-strand sequencing, illustrated in Figure 11, each fragment must be considered in two orientations: the original nucleotide sequence (forward), and the reverse-complemented nucleotide sequence (reverse):

Definition 1.2: Fragment orientation set

Denote by {f, r} the set of orientations, where f and r stand for the forward and the reverse orientation, respectively. Those orientations are complementary. The overline function • gives the complementary orientation (its reverse), i.e. f = r and r = f . The set of oriented fragments results from the Cartesian product of F un and {f, r}:

Definition 1.3: Oriented fragment sets

Let F = F un × {f, r} = F f F r be the set of oriented fragments, associated with the projection functions π Fun : F → → F un and f or : F → → {f, r}, Σ F = 0, 2|F un | be the set of fragment integer identifiers and f id : F → → Σ F be the associated labelling function. Let denote by seq : F → Σ nuc + the function that gives the nucleotide sequence of an oriented fragment. It holds that:

-F f is the set of forward fragments, such that for each a f ∈ F f : -f or(a f ) = f ; -seq(a f ) = unseq(π Fun (a f )).
-F r is the set of reverse fragments, such that for each a r ∈ F r :

f or(a r ) = r; seq(a r ) = unseq(π Fun (a r )).

Where • gives the reverse-complement of a given nucleotide sequence.

Unless stated otherwise, we will use the term fragment to denote an oriented fragment and precise when talking about an unoriented one.

Definition 1.4: The fragment reverse operation

For a given fragment a ∈ F, a denotes its reverse, where:

-f id(a) -f id(a) = 1; -f or(a) = f or(a) (a

Link set

We denote by the links the succession relationships between the oriented fragments. Figure 14 gives three examples of links exploited in the fragment assembly stages.

Property 1.1: Link reverse symmetry

One link between two fragments implies one link between the reversed fragments:

∀ (u, v) ∈ F 2 , m L ((u, v)) = m L ( v, u ) i.e. (u, v) ∈ L ⇐⇒ v, u ∈ L
Analogously to Definition 1.4, Definition 1.6 extends the reverse operation to the links:

Definition 1.6: The reverse operation for links

For a given link l ∈ L, where linkp(l) = (u, v), l denotes its reverse such that

linkp(l) = v, u = (u, v) = linkp(l). F un Σ nuc + {f, r} Σ Fun F un × {f, r} = F = F f F r Σ F F f L ⊂ F 2 F r Σ L L N >0 unseq unf id f id f or seq π Fun m L lid linkp

Directed graph (DG): oriented fragments based

The directed graph (digraph) structure (DG) was the first idea to represent fragments and their links. As far as we know, this structure was first mentioned by [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF].

For the read assembly stage, DG is described in FALCON [START_REF] Chin | Phased diploid genome assembly with single-molecule real-time sequencing[END_REF], HINGE [START_REF] Kamath | HINGE: Long-read assembly achieves optimal repeat resolution[END_REF], Shasta [START_REF] Shafin | Nanopore sequencing and the Shasta toolkit enable efficient de novo assembly of eleven human genomes[END_REF] and hifiasm [START_REF] Cheng | Haplotyperesolved de novo assembly using phased assembly graphs with hifiasm[END_REF]. Concerning the scaffolding stage, GAT refers to an unitig graph [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF].

Definition 1.7: Oriented fragments based directed multigraph (DG)

Let DG = (V , E, Φ) a directed multigraph (multidigraph) such that:

-V is the set of vertices, where each v ∈ V represents one of the two orientations for a fragment in F;

-E is the multiset of edges, where each e ∈ E represents one link in L;

-Φ : E → {(u, v) | u, v ∈ V , u = v}
is the incidence function that associates an edge with an ordered pair of vertices.

Denote by f rag : V → → F the bijective function that associates each vertex with one and only one fragment, and by link : E → → L the bijective function that associates each edge with one and only one link. Since one vertex is associated with only one fragment, Definition 1.8 extends the reverse operations given in Definitions 1.4 and 1.6 for DG.

a) (a f , b f ) ∈ L (b) (b f , a f ) ∈ L (c) (a f , b r ) ∈ L (d) (b r , a f ) ∈ L

Definition 1.8: The reverse operation in DG

For each vertex v ∈ V and edge e ∈ E, v ∈ V and e ∈ E denote their reverse. From the bijective properties of functions f rag and link, it holds that:

-f rag(v) = f rag(v); -link(e) = link(e).
Definition 1.9 specifies the conditions for a valid path in DG.

Definition 1.9: Path in DG

A path p = (v 0 , v 1 , . . . , v n-1 ) ∈ V n , n ∈ N in DG = (V , E, Φ) is valid if the following two properties hold: Contiguity ∀ i ∈ 0, n -1 , ∃ e ∈ E | Φ(e) = (v i , v i+1 ) Exclusive orientation ∀ v ∈ p, v / ∈ p
Note that a path in DG satisfies the same properties as for a generic multidigraph, except that a vertex and its reverse cannot both participate in the path. Figure 17 gives an example of a valid path in DG.

Figure 17 -A path in DG. p = v f u f x r is a valid path coloured in red. p = x f u r v r
is its reverse and is also valid. However, p and p are exclusive.

Finally, Property 1.3 defines the reverse symmetry in DG, that inherits from Property 1.1 and illustrated in Figure 16.

Property 1.3: Reverse symmetry in DG -∀ a ∈ F, ∃! v ∈ V | f rag(v) = a. a -∀ l ∈ L, ∃! e ∈ E | link(e) = l. a ∃! x ∈ X means "there is a unique x in X".

Proof

From Definition 1.7, f rag : V → → F and link : E → → L are bijective.

Bidirected graph (BG): oriented walk based

The bidirected graph (bigraph) was introduced by [START_REF] Edmonds | Matching: A well-solved class of integer linear programs[END_REF]. [START_REF] Myers | Toward simplifying and accurately formulating fragment assembly[END_REF] is the first to suggest a bigraph structure (BG) to store the links exploiting the reverse symmetry of the fragment and the link sets F and L.

The key idea is to represent both a fragment and its reverse by only one vertex. Because the links are ordered pairs of oriented fragments, the edges are bidirected. Definition 1.10 provides a formal description, simpler than the one given in [START_REF] Gritsenko | GRASS: A generic algorithm for scaffolding next-generation sequencing assemblies[END_REF] and suggested by [START_REF] Myers | Toward simplifying and accurately formulating fragment assembly[END_REF].

For the read assembly stage, BG is described in Bambus [START_REF] Pop | Hierarchical Scaffolding With Bambus[END_REF], Minimus [START_REF] Sommer | Minimus: A fast, lightweight genome assembler[END_REF] and in Edena [START_REF] Hernandez | De novo bacterial genome sequencing: Millions of very short reads assembled on a desktop computer[END_REF]. Concerning the scaffolding stage, BG is found in SOPRA [START_REF] Dayarian | SOPRA: Scaffolding algorithm for paired reads via statistical optimization[END_REF], Bambus 2 [START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF], MIP [START_REF] Salmela | Fast scaffolding with small independent mixed integer programs[END_REF], Opera [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF], SSPACE [START_REF] Boetzer | Scaffolding pre-assembled contigs using SSPACE[END_REF], GRASS [START_REF] Gritsenko | GRASS: A generic algorithm for scaffolding next-generation sequencing assemblies[END_REF], BOSS [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF], LACHESIS [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF], ScaffoldScaffolder [START_REF] Bodily | ScaffoldScaffolder: Solving contig orientation via bidirected to directed graph reduction[END_REF] and SLR [START_REF] Luo | SLR: A scaffolding algorithm based on long reads and contig classification[END_REF].

Definition 1.10: Oriented walk based bidirected multigraph (BG)

Let BG = (V , E, Φ, attre) a bidirected multigraph (multibigraph) such that:

-V is the set of vertices, where each v ∈ V represents the two orientations for a fragment in F;

-E is the multiset of edges, where each e ∈ E represents one link and its reverse in L;

-Φ : E → {{u, v} | u, v ∈ V , u = v} is the incidence function that associates
an edge with a pair of vertices;

attre : E → {f, r} 2 gives the orientations of the vertices in the associated link, and in the lexicographic order of the unoriented fragment identifiers.

Denote by unf rag : V → → F un the bijective function that associates each vertex with one and only one unoriented fragment, and by:

link : {(u, v, e) | e ∈ E, Φ(e) = {u, v}} → → L
the bijective function that associates each edge and its ordered vertices with one and only one link.

While visually each edge has two extremities in the multibigraph, Definition 1.10 describes the graph as an undirected multigraph (multiungraph). Figure 18 provides the link representation in BG, where each subfigure gives the multibigraph and the multiungraph view.

(a) (a f , b f ) ∈ L (b) (b f , a f ) ∈ L (c) (a f , b r ) ∈ L (d) (b r , a f ) ∈ L Figure 18 -Link cases in BG.
Fragments a f , b f and their reverse a r , b r can link in four manners. Their associated vertices are u and v. For each subfigure, the first graph has bidirected edges and the second presents the same information but with undirected edges where the attributes are given by the function attre. Here we assume that unf id(a) < unf id(b). In (a) and (c), attre gives the orientations in the order of the fragments in the described link. In (b) and (d), attre gives the orientations in the order of the fragment in the reverse of the described link ({r, r} for (a r , b r ) and {r, f } for (a r , b f ), respectively).

Property 1.4: Sizes of BG

For BG = (V , E, Φ, attre), |V | = 1 2 |F| and |E| = 1 2 |L|.
Less explicitly than in Definition 1.8, Definition 1.11 extends the reverse operations given in Definitions 1.4 and1 Definition 1.12 specifies the conditions for a path in BG to be valid, illustrated in Figure 19.

Definition 1.12: Path in BG

A path p = (v 0 , v 1 , . . . , v n-1 ) ∈ V n , n ∈ N in BG = (V , E, Φ, attre) is valid if the following two properties hold: Contiguity ∀ i ∈ 0, n -1 , ∃ e i , e i+1 ∈ E s.t.: -Φ(e i ) = {v i-1 , v i } ∧ Φ(e i+1 ) = {v i , v i+1 } -link(v i-1 , v i , e i )[1] a = link(v i , v i+1 , e i+1 )[0] Exclusive orientation ∀ v i , v j ∈ p, v i = v j =⇒ i = j a The notation x[i] denotes the i th element in x.
Finally, Property 1.5 defines the reverse symmetry in BG, that inherits from Property 1.1. 1.4 Undirected graph (UG): tail-head fragments based [START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF] first describe the undirect graph (ungraph) structure (U G) for the fragment and their links, in the context of contig scaffolding. The idea is to keep the merge of the two strands as in BG, but with a simpler description of the orientations for the fragments. Each fragment is split by a tail and a head. Both the tail and the head are vertices, and they are connected by an edge. Passing through first the tail and then the head corresponds to choosing the fragment in its forward orientation, while passing through first the head and then the tail corresponds to choosing it in its reverse orientation. This new type of edges are called fragment-edges, at the opposite of link-edges that correspond to the links. A path in the graph must alternate between fragment-edges and link-edges.

Property 1.5: Reverse symmetry in BG

∀ l ∈ L, ∃! e ∈ E, Φ(e) = {u, v} | link(u, v, e) = l link(v, u, e) = l.
The structure mainly appears in method for the scaffolding stage, as in SCARPA [START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF], Scaftools [START_REF] Briot | An integer linear programming approach for genome scaffolding[END_REF], BESST [START_REF] Sahlin | BESST -Efficient scaffolding of large fragmented assemblies[END_REF], fragScaff [START_REF] Adey | In vitro, long-range sequence information for de novo genome assembly via transposase contiguity[END_REF], in [START_REF] Chateau | A complexity and approximation framework for the maximization scaffolding problem[END_REF] and [START_REF] Weller | Exact approaches for scaffolding[END_REF], ScaffMatch [START_REF] Mandric | ScaffMatch: Scaffolding algorithm based on maximum weight matching[END_REF], SALSA2 [START_REF] Ghurye | Integrating Hi-C links with assembly graphs for chromosome-scale assembly[END_REF], LRScaf [START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF], in [START_REF] Davot | On a greedy approach for genome scaffolding[END_REF] and [START_REF] Aganezov | Orienting Ordered Scaffolds: Complexity and Algorithms[END_REF]. Concerning the read assembly stage, U G can be found in [START_REF] Myers | The fragment assembly string graph[END_REF] and [START_REF] Li | Minimap and miniasm: Fast mapping and de novo assembly for noisy long sequences[END_REF].

Definition 1.13: Tail-head fragments based undirected multigraph (U G)

Let U G = (V , E F , E L , Φ L ) be an multiungraph such that:

-V is the set of vertices, where each v ∈ V represents the tail or the head of a fragment;

-E F = {{v t , v h } | v t , v h ∈ V , v t = v h }
is the set of fragment-edges (v t for the tail, v h for the head), where each e ∈ E F represents both the forward and the reverse orientation of one fragment;

-E L is the multiset of link-edges, where each e ∈ E L represents both one link and its reverse;

-Φ L : E L → {{u, v} | u, v ∈ V , u = v}
is the incidence function that associates a link-edge with two vertices that are not connected by a fragment-edge.

Denote by:

f rag : {(u, v) | {u, v} ∈ E F } → → F
the bijective function that associates each ordered fragment-edge with one and only one (oriented) fragment, and by:

link : {(u, v, e) | e ∈ E L , Φ L (e) = {u, v}} → → L
the bijective function that associates each edge and its ordered vertices with one and only one link.

To complete Definition 1.13, as in [START_REF] Chateau | A complexity and approximation framework for the maximization scaffolding problem[END_REF], we can define the set E F of fragment-edges as a perfect matching of the vertices. Figure 20 illustrates each link case in U G.

Property 1.6: Sizes of U G

For U G = (V , E F , E L , Φ L ), |V | = |F|, E F = 1 2 |F| and |E L | = 1 2 |L|.
Analogously to Definitions 1.8 and 1.11, Definition 1.14 formalises the reverse operation in U G.

(a) (a f , b f ) ∈ L (b) (b f , a f ) ∈ L (c) (a f , b r ) ∈ L (d) (b r , a f ) ∈ L Figure 20 -Link cases in U G. u t ,
u h and v t , v h are the tail and the head of fragments a and b, respectively. The double edges correspond to fragment-edges while the plain edges correspond to link-edges. Each subfigure illustrates one link case in U G.

Definition 1.14: The reverse operation in

U G -Let {u, v} ∈ E F . Hence f rag(v, u) = f rag(u, v). -Let e ∈ E L , Φ L (e) = {u, v}. Hence link(v, u, e) = link(u, v, e).
Definition 1.15 specifies the conditions for a path in U G to be valid, illustrated in Figure 21.

Definition 1.15: Path in

U G A path p = (v 0 , v 1 , . . . , v 2k-1 ) ∈ V 2k , k ∈ N in U G = (V , E F , E L , Φ L ) is valid if the following two properties hold:
Contiguity fragment-edges and link-edges alternate in the path:

-odd-numbered edges in the path are fragment-edges:

∀ i ∈ 0, k , {v 2i , v 2i+1 } ∈ E F
-even-numbered edges in the path are link-edges:

∀ i ∈ 1, k , ∃ e ∈ E L | Φ L = {v 2i-1 , v 2i } Exclusive orientation ∀ v i , v j ∈ p, v i = v j =⇒ i = j
Note that walking in U G is not similar that walking in a classic ungraph, because a valid walk must alternate between fragment-edges and link-edges. Figure 21 gives an example of a valid path in U G.

Figure 21 -A path in U G.

A valid path is coloured in red. As U G is undirected, the resulting path can be read in two possible ways: from vertex v t or from vertex x t . Let p v = v t v h u t u h x h x t the one that begins from v t , and p x = x t x h u h u t v h v t the one that begins from x t . Note that p v is valid if and only if p x is valid.

Finally, Property 1.7 defines the reverse symmetry in U G, that inherits from Property 1.1.

Property 1.7: Reverse symmetry in

U G -∀ a ∈ F, ∃! {u, v} ∈ E F s.t.: f rag(u, v) = a f rag(v, u) = a -∀ l ∈ L, ∃! e ∈ E L , Φ L (e) = {u, v} s.t.: link(u, v, e) = l link(v, u, e) = l Proof From Definition 1.13, f rag : {(u, v) | {u, v} ∈ E F } → → F and link : {(u, v, e) | e ∈ E L , Φ L (e) = {u, v}} → → L are bijective.

Quick reference

For a quick overview of the different fragment graph structures, the readers can refer to Figure 22.

DG BG UG

Figure 22 -Overview of all the fragment graph structures. Let a and b two fragments. Each vertical dashed line separates two graph structures. Each colour is associated with an link case (and its reverse):

(a f , b f ) ∈ L in blue, (b f , a f ) ∈ L in green, (a f , b r ) ∈ L in yellow and (b r , a f ) ∈ L in violet. DG Directed fragment graph.
Each fragment is represented by two vertices for the two orientations. BG Bidirected fragment graph. Each fragment is represented by only one vertex. The first column gives the bidirected view while the second provides the undirected graph view with edge attributes (described at the bottom). UG Undirected fragment graph. Each fragment is represented by two vertices for its tail and its head, connected by a fragment-edge (double edges). Link-edges are simple edges.

Scaffolding the contigs

Recall from Chapter I, Section 5.4 that the first step of fragment assembly is to produce contigs from the reads. The second stage, the scaffolding, aims to orient and order the contigs. It results in a set of scaffolds, such that each of them is an order of oriented contigs separated by gaps, with possible known nucleotide length.

We first categorise the type of scaffolding input data in Section 2.1 and their subsampling strategies in Section 2.2. In Sections 2.3 to 2.5 we describe the scaffolding formulations proposed in the literature and the approaches to solve them, exactly or heuristically.

As the scaffolding stage was originally designed to use the mate-pair and pairedend data (see Chapter I, Section 3 and Figure 9), the latter sections are especially focusing on methods that take them as input. Nevertheless, the way they address the scaffolding problem formulations enables to generalise the orienting and the ordering of the contigs, under the assumption that different sequencing data can provide the same link information, as discussed in the first section. Table 1 lists the scaffolding methods described in this section. The three graph representations (DG, BG, U G) can be compared in Figure 22.

Scaffolding input data

We propose to discriminate the scaffolding input data according to two properties. The first one is a link property: the data consist in a (multi)set of ordered pairs of oriented contigs (c.f. Section 1.1.2). The second one is a nucleotide distance property: in addition to their order, the oriented contigs are spaced appart by an estimated nucleotide distance.

As the contigs are the read assembly results, it is necessary to align the reads (or the sequencing data) against the contigs. In fact, especially with the DBG read assembly approach, the reads permitting to build a contig are unknown. Figure 23 illustrates the two properties with a sequencing data centric point of view, while Table 2 provides a property centric view.

Paired-end and Hi-C reads The reads come by pairs and are sequenced from the two different strands. The reads are mapped against the contigs. Consider r i , r j to be two paired-end reads such that r i maps on contig c i and r j maps on contig c j . Therefore, oriented contig c i is before c j . Two cases: i. The two reads map two different contigs. The orientation of c i is given by the one of r i in the mapping, and the orientation of c j is the reverse of this of r j in the mapping. In the case of paired-end read and mate-pair reads (not for Hi-C), the distance between the oriented contig is determined according the mapping coordinates of the reads on the contigs and the distance between the paired reads.

ii. The two reads map the same contig. To be consistent, the reads must map in different orientations, otherwise it may indicate a misassembly of the two Table 1 -Non-exhaustive categorised list of methods and approaches for the scaffolding stage. Each row corresponds to a method dedicated to the scaffolding stage in the fragment assembly process. In the Graph column, DG, BG and U G respectively stand for directed, bidirected and undirected graphs. strands. For paired-end data, if the distance between the reads and the absolute difference of the two mapping coordinates do not coincide, then it also suggest a misassembly (insertion or deletion).

Dedicated input data

Long reads and reference The order and the orientations of the contigs is determined according their mapping against a longer sequence (a long read or a reference).

Optical map The genome is cut multiple times by a specific restriction site. Each cut produces a map that consists in a list of distances between the same nucleotide site. The same process is artificially produced for the contigs. Then, the merged maps of a contig is aligned with the merged maps given by the optical map. This produce an order between the oriented contig.

Linked reads Each pair of reads is tagged with a barcode. Two pairs with the same barcode means they come from the same fragment. It enables to cluster the contigs, without orienting or ordering them. One way to find an order between two oriented contigs is to split each contig in two, and check which half parts of the contigs are covered by the reads with the same barcode.

Table 2 -Scaffolding input data properties. Each sequencing data is provided with at most the two properties, link and nucleotide distance.

Sequencing data Link Distance

Paired-end Mate pair Hi-C Long reads Reference Optical mapping Linked reads

Paired-end and mate pair data provide the same type of data except that the distances are larger for the second one but less confident. The barcoding feature of linked read itself does not provide distance information. However, each paired-end read is barcoded.

Subsampling the input data

Several difficulties arise for the scaffolding problem. First, the number of links between the contigs is the main factor for the complexity of orienting and ordering them. Merging several links into one drops the computational time and memory use. Moreover, the presence of sequencing error suggests filtering the result of the alignments, or the sequencing data themselves. Even after a filtering step, artefact alignments may occur because of repeats in the genome, and consequently may produce links resulting in misassemblies if they are chosen. It implies that c i forward is followed by c j forward, separated by a nucleotide distance d . (b) This case is the same as (a) except that for Hi-C data, the distance d is unknown, implying an unknown distance d . (c) c i and c j map on two overlapping long reads lr i and lr j . It implies that c i forward is followed by c j forward, separated by an approximated nucleotide distance d . (d) G is a reference genome (e.g. of a species near to the organism of interest). Mapping c i and c j on G provides an orientation of the contigs and an order with a distance d . (e) OP M is an optical map with three markers (yellow square, green diamond and purple circle). Finding the same markers orders in c i and c j enables to estimate their orientation and their order with a distance d . (f) r i , r j , and r k , r l are two paired-end reads with the same barcode (yellow circle), so they are sequenced from the same fragment. Because r i , r j map on the head of c i , and r k , r l map on the tail of c j , it is possible to orient and order c i and c j .

Bundling the links

When several alignments produce the same ordered pair of oriented contigs, they can be bundled. Thus, only one link represents them. This is the strategy that [START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF], [START_REF] Dayarian | SOPRA: Scaffolding algorithm for paired reads via statistical optimization[END_REF], [START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF], [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF], [START_REF] Gritsenko | GRASS: A generic algorithm for scaffolding next-generation sequencing assemblies[END_REF], [START_REF] Salmela | Fast scaffolding with small independent mixed integer programs[END_REF], [START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF], [START_REF] Briot | An integer linear programming approach for genome scaffolding[END_REF], [START_REF] Mandric | ScaffMatch: Scaffolding algorithm based on maximum weight matching[END_REF] and [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] adopt with paired-end data, [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] and [START_REF] Ghurye | Integrating Hi-C links with assembly graphs for chromosome-scale assembly[END_REF] for Hi-C data, and [START_REF] Adey | In vitro, long-range sequence information for de novo genome assembly via transposase contiguity[END_REF] and [START_REF] Coombe | ARKS: Chromosome-scale scaffolding of human genome drafts with linked read kmers[END_REF] for linked reads. The link representing the bundle is weighted proportionally to the number of links supporting the bundle. In the case of distances, the mean and the standard deviation is computed.

Removing paired-end reads

Detecting artefact data Low mapping score are ignored, as [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] does. If only one read of the pair is aligned, the whole pair is ignored [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF][START_REF] Coombe | ARKS: Chromosome-scale scaffolding of human genome drafts with linked read kmers[END_REF]. [START_REF] Mandric | ScaffMatch: Scaffolding algorithm based on maximum weight matching[END_REF] and [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] remove paired-end alignments from the bundle if the distance they provide between the contigs is over a confident interval. [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] compute the distance between the reads if the two contigs are concatenated, and remove the alignment of the pair if the distance is too long. Finally, because of the nature of Hi-C paired-end reads, [START_REF] Ghurye | Integrating Hi-C links with assembly graphs for chromosome-scale assembly[END_REF] trim the reads that are aligned after a ligation.

Detecting repeats When one of the two paired reads maps on distinct contigs, [START_REF] Mandric | ScaffMatch: Scaffolding algorithm based on maximum weight matching[END_REF], [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] and [START_REF] Ghurye | Integrating Hi-C links with assembly graphs for chromosome-scale assembly[END_REF] do not consider the pair. [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] ignore the paired-end reads when they map to a highly covered contig.

Removing contigs

Detecting artefact data [START_REF] Dayarian | SOPRA: Scaffolding algorithm for paired reads via statistical optimization[END_REF] remove a contig on which the two reads of a pair are mapped on it with the same orientations. They also compute the standard deviation of all the distances in all the contigs and remove those for which the reads' distance is out of the deviation. For Hi-C data, [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] remove a contig if it does not contain enough restriction sites (only for their clustering step). In the long read scaffolding method of [START_REF] Luo | SLR: A scaffolding algorithm based on long reads and contig classification[END_REF], short contigs are ignored. Finally, [START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF] only consider contigs that are contained in long reads.

Detecting repeats A contig is dismissed if its alignment coverage is over a threshold of beyond a distribution's confidence interval [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF][START_REF] Sahlin | BESST -Efficient scaffolding of large fragmented assemblies[END_REF][START_REF] Mandric | ScaffMatch: Scaffolding algorithm based on maximum weight matching[END_REF][START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF]. For Hi-C data, [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] identify a contig as a repeat if the number of restriction sites is high. [START_REF] Luo | SLR: A scaffolding algorithm based on long reads and contig classification[END_REF] remove a contig if it maps in the middle of different long reads, and if its neighbour contigs differ. [START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF] computes all the shortest path between two contigs: a contig appearing in too many paths is dismissed.

Removing links

Detecting artefact data A bundle with a low confidence is not considered. The lower threshold can be a fixed value, or based on the distribution of the confidence over all the bundles [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF][START_REF] Gritsenko | GRASS: A generic algorithm for scaffolding next-generation sequencing assemblies[END_REF][START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF][START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF][START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF]. Note that in [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF] this threshold is the result of a simulation. If the distribution of the distances on a link significantly differs from the distribution over all the bundles, then the link is dismissed [START_REF] Sahlin | BESST -Efficient scaffolding of large fragmented assemblies[END_REF][START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF] In [START_REF] Ghurye | Integrating Hi-C links with assembly graphs for chromosome-scale assembly[END_REF], the confidence of neighbours links is compared, and only the "best buddy weighted" link is kept (a O(|L|) procedure). In addition, as they take in input an assembly graph, the authors remove Hi-C links if there is no link in the assembly graph that match the contig orientations. [START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF] dismiss a contig alignment if the long read is contained in the contig. 2019) compute all the shortest path between two contigs and keep only a link if its ratio of presence in the paths is significantly greater than the presence of the other links involving the two same contigs (a O(|L|(|C| + |L|)) procedure). Finally, if two contigs are aligned against the same end coordinate of a long read, the best mapping is kept in [START_REF] Luo | SLR: A scaffolding algorithm based on long reads and contig classification[END_REF].

Reducing redundant links

Partitioning the instances

Biconnected components can be solved independently For BG and U G fragment graphs (Sections 1.3 and 1.4), finding articulation edges and removing them enables to partition the graph into biconnected components. [START_REF] Dayarian | SOPRA: Scaffolding algorithm for paired reads via statistical optimization[END_REF] show that each component can be solved independently without a loss of optimality. [START_REF] Salmela | Fast scaffolding with small independent mixed integer programs[END_REF] and [START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF] also apply this divide and conquer strategy. Finding all the articulation points can be done in O(|C| 2 )

Clustering according the number of chromosomes In the [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] Hi-C scaffolding method, the number of choromosomes N is an input. Based on the links' weight, the contig set is partionned in N parts with a hierarchical agglomerative clustering approach. Note that the clustering does not take into consideration the relative orientations of the contigs.

Orienting the contigs

Here we present the methods that focus on the contig orientation problem independently of their ordering. [START_REF] Kececioglu | Combinatorial algorithms for DNA sequence assembly[END_REF] have already defined it in the context of the read assembly. It consists in removing links or contigs in order to find only two orientation assignments, where one is the reverse of the other. A link set over a contig set is defined as an orientation-valid layout if it respects the latter definition.

⊂ = L ⊂ f L ⊂ r such that: -L ⊂ f = |L ⊂ r | and ∀ l ∈ L ⊂ f , l ∈ L ⊂ r ; -o and o are respectively consistent in L ⊂ f and L ⊂ r , i.e. ∀ l ∈ L ⊂ f where linkp(l) = (a, b), and so l ∈ L ⊂ r where linkp(l) = b, a , ∃! i, j ∈ 0, n , i = j such that: -a = (c i , o i ) ⇐⇒ a = (c i , o i ); -b = (c j , o j ) ⇐⇒ b = (c j , o j ).

Maximising the sum of used link weights

As originally proposed by [START_REF] Kececioglu | Combinatorial algorithms for DNA sequence assembly[END_REF], the orientation problem consists in removing links in order to get an orientation-valid layout that maximising the sum of the remaining links' weight. Here, no one of the contigs is removed, i.e. C = C. Equivalently, it consists in minimising the sum of removed links' weight. [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF] reduced this formulation to the maximum-weight cut problem whose decision version is N P-complete. [START_REF] Pop | Hierarchical Scaffolding With Bambus[END_REF] consider the weight proportional to the number of supporting paired-end reads in the bundle, or inversely proportional to the nucleotide distance of the bundle. [START_REF] Dayarian | SOPRA: Scaffolding algorithm for paired reads via statistical optimization[END_REF] aim to minimise a cost function (simply, the sign of the weight is inverse comparing to the other methods). In Bambus 2 [START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF] (as for its predecessor Bambus 2), the problem consists in removing the minimum number of reversal links, i.e. links where the two contigs in the pair are not identically oriented. [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] solve the problem iteratively by increasing cleaned subgraph. The first subgraph contains the heaviest links, and after resolving it for the orientation, lighter links are added etc. [START_REF] Luo | SLR: A scaffolding algorithm based on long reads and contig classification[END_REF] maximise the sum of remaining links' weight, where the weights correspond to the number of aligned long reads.

2.3.2 Remove the minimum number of contigs and links to avoid odd reversal cycles by minimising the weight sum

Maximum log-likelihood function

At the opposite of the previous methods, [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] first order the contig. Then, they orient the contig according to a maximum log-likelihood function, where each edge weight equals to the inverse of the number of Hi-C links normalised by the number of fragment sites per contig.

Ordering the oriented contigs

Once an orientation assignment is determined for the contigs, the next step is ordering them and positioning them when nucleotide distances are provided.

2.4.1

Maximising the sum of used distances bunches' weights [START_REF] Pop | Hierarchical Scaffolding With Bambus[END_REF] reduce the ordering problem to the optimal linear arrangement problem. It consists in assigning a coordinate on an axe to each contig, a N P-hard problem. As for the orientation problem, they consider to be the best weights either the number of supporting links in a bundle, or the smallest distances. However, they do not position the contigs by satisfying the distances, as done in [START_REF] Luo | SLR: A scaffolding algorithm based on long reads and contig classification[END_REF] for long read based scaffolding. [START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF] refer to the feedback arc set problem (N P-hard, [START_REF] Karp | Reducibility among Combinatorial Problems[END_REF]) and define the contig ordering problem as finding a minimal set of edges whose removal leads to a directed acyclic graph. [START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF] position the contigs in confidence intervals of three times the standard deviation.

Minimum spanning tree

For Hi-C data, [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] first order the contig. Recall that each link weight equals to the inverse of the number of Hi-C links normalised by the number of fragment sites per contig. They define the ordering problem as a sequence of process: (i) finding a minimum weight spanning tree; (ii) extracting the longest path (a trunk) to get a first order; (iii) ordering the branching contigs according to the trunk, from the longest to the smallest, with maximising the number of links in the bundle. [START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF] define the positioning problem as a sub-step of the scaffolding problem following the ordering of the contigs. According to the distance distribution in each bundle, they aim to give a nucleotide position for each contig, such that it minimises the distance approximation in each link. [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF] maximise a quadratic likelihood function on the distances to estimate the gaps in each scaffold. A unique solution can be found thanks to the Goldfarb-Idnani active-set dual method in polynomial time [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF].

Contig positioning only

Orienting and ordering the contigs simultaneously

Already for the read assembly stage, [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF] has suggested combining both the orienting and the ordering/positioning problems in one step. The purpose is to seek a global optimisation approach. Here we enumerate the scaffolding methods that follow this approach.

2.5.1 Maximising the sum of distance bundles' weights to order linearly the contigs [START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF] aim to maximise the sum of weights of happy mate-edges in U G. A mate-edge (a bundle of paired-end links) is happy if the orientations of its contigs respect it, and if the interval of the contigs' coordinates is in a confidence interval of three times the standard deviation from the distances median provided by the edges. They prove the corresponding decision problem to be N P-complete by reducing the bandwidth problem to it. [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF] focus on a bounded version of the graph bandwidth problem, for which the complexity is in O(|C| w |L| p+1 ) where w is the fixed width, and there are at most p discordant edges in a scaffold.

In their Mixed-Integer Linear Programming (MILP) approach, [START_REF] Salmela | Fast scaffolding with small independent mixed integer programs[END_REF] model if a link is kept by real variables between 0 and 1 instead of binaries. A value that approaches 1 indicates that two contigs' positions are close to the distance median provided by the link, while a value that approaches 0 indicates that the contigs' positions are far from the median. In the objective function, the links' weights are multiplied by these reals. [START_REF] Gritsenko | GRASS: A generic algorithm for scaffolding next-generation sequencing assemblies[END_REF] combine three goals in a single objective function: minimising the cost of orientation, minimising the orientation and the distance cost, and minimising the orientation and the order cost.

2.5.2

The heaviest matching, paths and cycles [START_REF] Briot | An integer linear programming approach for genome scaffolding[END_REF] maximise the weight of the remaining links such that at most one link connects each contig's end in U G. This approach can result in several paths. Building on this, [START_REF] Weller | Exact approaches for scaffolding[END_REF] are looking for at most σ p paths and σ c cycles (where σ p , σ c are parameters) that cover all the contig-edges in U G. They follow the scaffolding definition given in [START_REF] Chateau | Complexity and Polynomial-Time Approximation Algorithms around the Scaffolding Problem[END_REF], where the authors study the complexity of the problem on several classes of graphs and according to the number of contigs (n). In their paper, they are looking for exactly σ p paths and σ c cycles. Except in the case where n = σ p + 2σ n , for which the problem is in P, the scaffolding problem (decision version) remains N P-complete.

In the global Integer Linear Programming (ILP) formulation of [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF], the objective function aims to find the longest nucleotide sequence that maximises the number of satisfied distances. It consists in searching for an elementary path in DG. Combining these two objectives in one linear function echoes the SCS criticisms formulated in the case of read assembly (see Chapter I, Section 5.1). Note that the paired-end links are not especially bundled, and each of them are two links in DG. At the opposite of the other approaches, here the contigs are unitigs associated with a multiplicity, which is an upper bound of its use. The model is enough general to integrate other weight definitions (both on the overlaps and on the unitigs).

2.5.3

The maximum (weighted) matching and spanning tree [START_REF] Mandric | ScaffMatch: Scaffolding algorithm based on maximum weight matching[END_REF] reduce the scaffolding problem to the Maximum-Weight Acyclic 2-Matching (MWA2M) problem on the link-edges in U G. The latter one's decision version is N P-complete.

In their Hi-C scaffolding method, Ghurye et al. ( 2019) iteratively find the maximum-weighted matching of the link-edges. Then, they remove the lightest links to eliminate the cycle in the matching to order and orient the covered contigs. Finally, the covered vertices are removed and they re-iterate on the resulting subgraph. The whole procedure is in O(|L||C| 3 ).

Concerning the scaffolding with linked read, [START_REF] Adey | In vitro, long-range sequence information for de novo genome assembly via transposase contiguity[END_REF] adopt a similar strategy as in [START_REF] Burton | Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions[END_REF] for Hi-C data by finding a maximum-weighted minimum spanning tree, extract the longest path (trunk), and add the branching contigs to the linear ordering.

Finally, for long reads, [START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF] only extract valid simple paths in U G and return them as scaffolds.

Solving approaches

Greedy approaches

As in the majority of the cases the scaffolding problem (orienting and ordering separately or in the same time) is N P-hard, in practice some authors develop greedy algorithm to solve large instances.

Mate-pair in decreasing order of weight The main greedy strategy consists in prioritising the heaviest links. [START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF] are the first to present the greedy path merging algorithm. [START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF] adopt this strategy for the orienting and the ordering of the contigs separately. They use the same algorithm as in [START_REF] Kececioglu | Combinatorial algorithms for DNA sequence assembly[END_REF]: a O(|C| + |L|) complexity algorithm on BG. [START_REF] Mandric | ScaffMatch: Scaffolding algorithm based on maximum weight matching[END_REF] propose a greedy heuristic for the MWA2M problem that chooses the heaviest feasible link-edge in U G. A link-edge is said feasible if it does not make a vertex degree higher than two and does not form cycles with the previously chosen edges. Their algorithm is in O(|C|log(|C|)) with a max heap implementation. In the continuity of [START_REF] Chateau | A complexity and approximation framework for the maximization scaffolding problem[END_REF], [START_REF] Davot | On a greedy approach for genome scaffolding[END_REF] propose a O(|C||L|σ 2 c ) 3-approximation algorithm for connected cluster graphs (U G), where σ c is the upper bound of cycles.

Seed-and-extend As for the read assembly stage, Seed-and-Extend (S&E) algorithms have been developed for the scaffolding stage. They ensure to form long scaffolds through the extension of the current paths. [START_REF] Boetzer | Scaffolding pre-assembled contigs using SSPACE[END_REF] start with the largest contig and extend it with links supporting at least a given threshold. In the case of alternative links, they try to find an ordering of all the alternative adjacent contigs with the distance attributes. If it fails, they compute a ratio of the weight of the alternatives, and choose the best link if its ratio is above a fixed threshold. As in [START_REF] Pop | Hierarchical Scaffolding With Bambus[END_REF], [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] first order the longest contigs and connect them with a linear BFS procedure. [START_REF] Dayarian | SOPRA: Scaffolding algorithm for paired reads via statistical optimization[END_REF] use a simulated annealing which is a Monte Carlo approach. By decreasing a parameter (the temperature) in the sample's weight, if the energy of the system reaches a value close to the minimum, most of the constraints are solved. This approach is applied for both the orientation and for the ordering, independently.

Simulated annealing

Spanning tree As described in the previous sections, [START_REF] Roy | SLIQ: Simple Linear Inequalities for Efficient Contig Scaffolding[END_REF] implement a spanning tree algorithm for the orienting of the contigs while [START_REF] Adey | In vitro, long-range sequence information for de novo genome assembly via transposase contiguity[END_REF] implement it for both the orienting and the ordering of the contigs. 2.6.2 Fix parameters dynamic programming [START_REF] Gao | Opera: Reconstructing Optimal Genomic Scaffolds with High-Throughput Paired-End Sequences[END_REF] propose a O(|C| w |L| p+1 ) procedure, where w is the width and p is the upper bound of discordant edges in BG. They increase the parameters until a solution is found. [START_REF] Weller | Exact approaches for scaffolding[END_REF] give an exact algorithm in O(|C|σ 2 p ) for trees (U G graph), where σ p is the maximum number of paths to cover with the matching. By extension, they describe an exact tree decomposition algorithm in O(tw tw σ p σ c |C|), where tw is the tree width for the decomposition, and σ c is equivalent to σ p for cycles.

Majority voting on predictions

Mathematical programming

Mathematical programming enables to model and to solve complex optimisation problems by using dedicated solvers.

Mixed Integer Linear Programming (MILP) [START_REF] Salmela | Fast scaffolding with small independent mixed integer programs[END_REF] use a MILP approach on BG, where the link-choice variables are not binaries but are reals between 0 and 1. It enables to model the uncertainty of the link during the positioning of the contigs. Recall that they use MILP on each bi-connected component found by removing articulation vertices. [START_REF] Briot | An integer linear programming approach for genome scaffolding[END_REF] however model the link choice variables as binaries and connect a vertex with at most one link in U G. [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] maximise the sum of weights of used links for orienting the contigs on iteratively increasing subgraphs. [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF] model in one objective function the longest path and the maximum number of satisfied constraints on distances. [START_REF] Luo | SLR: A scaffolding algorithm based on long reads and contig classification[END_REF] apply the same strategy as in [START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] on both the orienting and the ordering problems.

Mixed Integer Quadratic Programming (MIQP) [START_REF] Gritsenko | GRASS: A generic algorithm for scaffolding next-generation sequencing assemblies[END_REF] are the only ones that use Mixed-Integer Quadratic Programming (MIQP) by merging orienting and ordering the contigs in one objective function.

Linear Programming (LP) [START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF] use linear programming to position the oriented and ordered contigs to minimise the positioning approximation according to the distance distributions on the links.

Chloroplast genome assembly

As the chloroplast is one of several living organelles in the same plant cell, the DNA sequencing process produces reads from chloroplasts, other organelles and nuclear DNA. Consequently, methods have been developed to produce a separate assembly of the chloroplasts. They mainly consist in an assembly pipeline and adopt various strategies in the filtering of the data during the whole assembly process.

In the following, all the described methods, listed in Table 3, are based on short paired-end read data. Recently, [START_REF] Freudenthal | A systematic comparison of chloroplast genome assembly tools[END_REF] have compared some of them on several criteria, from user-friendly aspects to assembly qualities. Here, we decompose the dedicated chloroplast genome methods that solve the whole or a part of the fragment assembly process.

Three major issues arise from the literature: (i) How to extract the chloroplast data from the mixture of genomic material in the sequencing output? (ii) What is the most suitable assembly formulation for the chloroplast genomes? (iii) How to address the multimeric forms in the assembly results? Section 3.1 covers point (i), Sections 3.2 and 3.3 cover point (ii) and Section 3.4 covers point (iii).

Chloroplast sequence extraction

Although there are routine data filters, as in [START_REF] Bakker | Herbarium genomics: Plastome sequence assembly from a range of herbarium specimens using an Iterative Organelle Genome Assembly pipeline[END_REF] and McKain and afinit (2017), here we discuss specific approaches to identifying chloroplast data Table 3 -Chloroplast genome assembly approaches. Approach: the main contribution given in the paper. By Seed-and-Extend we mean a read assembly method. Filter: an arbitrary score (from 0, no filter, to 3 bullets) equivalent to the number of filter steps and the variability of their strategies. Multimerism: if the method is chloroplast multimeric forms aware, i.e. if it addresses the issue of several chloroplast forms in the sequencing data. Cannot find a description of the method, software manual not found.

(original and transformed) and extracting them.

Filtering the reads

De novo filtering Before any assembly of the raw data, [START_REF] Ankenbrand | chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data[END_REF] assume that the chloroplast (and mitochondria) reads outnumber the nuclear reads. Therefore, they keep reads for which the k-mer distribution is higher comparing to the other reads (thanks to Jellyfish, Marçais and Kingsford 2011).

Genome reference-based filtering McKain and afinit (2017) and [START_REF] Sancho | Comparative plastome genomics and phylogenomics of Brachypodium: Flowering time signatures, introgression and recombination in recently diverged ecotypes[END_REF] align the reads against one reference chloroplast genome, [START_REF] Ankenbrand | chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data[END_REF] against several closely related species, [START_REF] Bakker | Herbarium genomics: Plastome sequence assembly from a range of herbarium specimens using an Iterative Organelle Genome Assembly pipeline[END_REF] against several but not necessarily closely related. [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] in addition align the reads against organelle subsequences. In the latter, they denote the reads they extract as seed reads. Note that in [START_REF] Sancho | Comparative plastome genomics and phylogenomics of Brachypodium: Flowering time signatures, introgression and recombination in recently diverged ecotypes[END_REF] the use of a reference genome is optional.

Baiting new reads from assembly Some methods continue to filter the data after assembling a subset of the reads into contigs. [START_REF] Bakker | Herbarium genomics: Plastome sequence assembly from a range of herbarium specimens using an Iterative Organelle Genome Assembly pipeline[END_REF] iteratively align the previously unaligned reads to a running read assembly. Then they compute a new assembly on the larger read subset and repeat the process until no new read is extracted. [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] extend the seed contigs (from the assembly of the seed reads) with the remaining reads.

Filtering the contigs

Similarly to the reads, the contig filters follow two approaches.

De novo filtering McKain and afinit (2017) and [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] keep the contigs with a high read mapping coverage. Especially, the latter authors partition the contig set in parts representing each organelle type according to a coverage distribution described by a Gaussian mixture.

Label database filtering In addition, [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] align each contig to a label database, where each sequence (e.g. of a gene, a protein or a conserved region) is associated with a specific organelle (e.g. chloroplasts). They only keep the contigs that match the label database.

Chloroplast reads assembly

The assembly of the read is the first stage of the fragment assembly process.

However, as we have seen in the previous subsection, assembling the reads often enables to bait more reads and thus raise the chances to have a complete chloroplast assembly.

De Bruijn graph approach

Recall that DBG is a dbgraph read assembly approach.

Baiting assemblies As mentioned before, [START_REF] Bakker | Herbarium genomics: Plastome sequence assembly from a range of herbarium specimens using an Iterative Organelle Genome Assembly pipeline[END_REF] iteratively produce several read assemblies with several sizes of k-mer (recall Chapter I, Section 5.3) thanks to SOPAdenovo2 [START_REF] Luo | SOAPdenovo2: An empirically improved memory-efficient short-read de novo assembler[END_REF]. At each iteration they keep the assemblies where the N50 criterion is high, and bait new reads on selected assemblies (i.e. contigs).

Assembling the reads [START_REF] Sancho | Comparative plastome genomics and phylogenomics of Brachypodium: Flowering time signatures, introgression and recombination in recently diverged ecotypes[END_REF] use Velvet [START_REF] Zerbino | Velvet: Algorithms for de novo short read assembly using de Bruijn graphs[END_REF], whereas [START_REF] Bakker | Herbarium genomics: Plastome sequence assembly from a range of herbarium specimens using an Iterative Organelle Genome Assembly pipeline[END_REF], [START_REF] Ankenbrand | chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data[END_REF] and [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] prefer SPAdes that internally tests different size of k-mer [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to singlecell sequencing[END_REF]. Unlinke the latter that use paired-end for the read assembly, McKain and afinit (2017) use SPAdes without the paired-end option. Indeed, paired-end data are further exploited during the scaffolding stage.

Seed-and-extend

The S&E approach is a local assembling approach adapted for "simple" genomes, particularly with few repeats. [START_REF] Coissac | From barcodes to genomes: Extending the concept of DNA barcoding[END_REF] implement a S&E whose description could not be found at the time of writing this thesis. The S&E strategy in [START_REF] Dierckxsens | NOVOPlasty: De novo assembly of organelle genomes from whole genome data[END_REF] consists in: (i) extending the seed until a circularity is found; (ii) verifying the paired-end distance constraint at each extension; (iii) identifying repeats when extensions are not overlapping. In the case it cannot resolve the repeats based on more stringent extension parameters, the extension is cut, and it begins a new seed.

Chloroplast scaffolding

The scaffolding stage follows the assembly of the reads. Few chloroplast assembly tools tackle this stage.

Paired-end scaffolding McKain and afinit (2017) and [START_REF] Sancho | Comparative plastome genomics and phylogenomics of Brachypodium: Flowering time signatures, introgression and recombination in recently diverged ecotypes[END_REF] use the paired-end based scaffolder SSPACE [START_REF] Boetzer | Scaffolding pre-assembled contigs using SSPACE[END_REF]. Finally, the latter fill the gap in the scaffolds with GapFiller [START_REF] Boetzer | Toward almost closed genomes with GapFiller[END_REF]. [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF] test their scaffolder on chloroplast data. As described in Section 2.5.1, they aim to find the longest DNA sequence by finding a unitig path that maximises the number of satisfied distance constraints.

Estimating contig multiplicities to obtain a circular pattern [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] aim to determine the multiplicity of each contig such that Eulerian circuits are found in the resulting graph. They use a MIQP aiming to minimise the squared difference between the observed and the variable multiplicities. Each equation expresses the multiplicity of a contig's end with those of its adjacent contigs' ends. They compute the observed multiplicities according to the contigs' coverages.

Chloroplast assembly validation

Validation of assemblies can be carried out either by conventional evaluation processes, or by processes more specific to chloroplasts. Furthermore, this last stage can enable to detect multiple solutions especially due to the multimeric forms of the chloroplast genome.

De novo evaluation [START_REF] Bakker | Herbarium genomics: Plastome sequence assembly from a range of herbarium specimens using an Iterative Organelle Genome Assembly pipeline[END_REF] select the best assemblies according to the Assembly Likelihood Estimation (ALE) [START_REF] Clark | ALE: A generic assembly likelihood evaluation framework for assessing the accuracy of genome and metagenome assemblies[END_REF]. The assembly candidates come from both the final assemblies during the iterative read baiting step and the final step.

Reference genome evaluation [START_REF] Ankenbrand | chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data[END_REF] keep the contigs that align sufficiently against reference genomes.

Multiple solutions [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF] detect in the DG path, representing the scaffold, pairs of subpaths p = (v 0 , . . . , v n-1 ), p = v 0 , . . . , v n-1 such that:

∀ i ∈ 0, n , f rag(v i ) = f rag(v n-1-i ),
i.e. they detect subpaths that represent inverted repeats (IR). They cut the scaffold at the beginning and at the end of these subpaths in order to consider alternative paths representing the multimeric forms. They prove the alternative paths to be optimal. Jin et al. ( 2020) produce all the paths consuming all the multiplied contigs. They keep the paths that represent iso-IR, i.e. paths that have the same property as in [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF]. 22). For the most interested ones, Section 2 provides algorithms to explore the graph, to add or to delete vertices and edges. Finally we conclude on the memory and the time costs of the proposed implementations and their associated algorithms. For a quick overview, Figure 30 is a map that groups the graph implementations according to their structures, and summarises their properties.

III FRAGMENT GRAPH IMPLEMENTATIONS AND COMPARISON

Implementations

The implementations should satisfy both querying and dynamic updating requirements:

Querying requirements

-given an (oriented) fragment, getting all the fragments linking it;

-given a link, answering true if it is represented in the graph.

Dynamic updating requirements

-adding a fragment/link; -deleting a fragment/link.

A square matrix can represent the fragments and their links. The matrix is sparse in practice, e.g. in the read assembly stage, the reads are overlapping only if they are sequenced from closed genomic region or from repeats. A Compressed Sparse Row (CSR) or a Compressed Sparse Column (CSC) is suitable to address the sparsity. It also enables fast link querying. However, dynamic operations are costly. Thus, we implement the graphs with adjacency lists. They are suitable for both the sparsity and the efficiency of addition or deletion operations. However, as the number of neighbours is not a priori known, adjacency lists require the use of pointers to memory address, that are costly in memory.

The implementations require suitable graph labelling strategies: the vertices and the edges are labelled with indices, such that their ranges are minimised according to our purpose. For the sake of clarity, in the next sections we consider all the fragments and links are represented. Since for each graph structure there are bijective functions to the fragment set F and to the link set L, there must be bijective functions from the implementations to the graph structure. In that case, we will ensure the implementations represent the data.

Common notations

Canonical link

We described what a forward and a reverse fragments are (F = F f F r ) but we have not yet split the link set into two. Analogously to the fragment set, we partition the link set in two, L = L f L r . L f contains the links for which in the pair, the fragment with the smallest Σ F label is in forward orientation, i.e. L f = l ∈ L | f or(a) = f , f id(a) = min x∈linkp(l) f id(x) . Similarly to the unoriented fragment label set Σ Fun , let Σ Lcan = 0, 1 2 |L| be the label set of canonical links and canlid : L → → Σ Lcan its associated function, such that ∀ l ∈ L, canlid(l) = canlid(l). In this chapter, a link label represents both a link and its reverse, i.e. a canonical link. The order of the fragments enables distinguishing if the link label corresponds to a link in L f or in L r .

Adjacency list notation The adjacency lists are represented by functions f : N → X, where X can be the set of neighbour lists (themselves represented by functions) or the set of neighbours. From a programming point view, f [i] is the object at the i th position in the list f . Here we denote this object by f (i) except in the pseudocodes. We denote by N V and N v the adjacency lists for all the vertices and for the particular vertex v, respectively. N -and N + stand for the predecessor and the successor lists.

In/out degrees For a vertex v ∈ V , we denote by

d v = d - v + d +
v its total degree equals to the sum of the number of its predecessors and successors.

Octet-isation For n ∈ N, oct(n) = log 2 n 8 gives the number of octets necessary to memorise an integer between 0 and n -1.

Directed graph (DG): oriented fragments based

Section 1.1.1 describes the first implementation for DG, while Sections 1.1.2 and 1.1.3 propose clever implementations thanks to the reverse symmetry. Recall from Definition 1.7 that DG = (V , E, Φ).

All oriented fragments directed graph (DGA)

DGA is a classical multidigraph implementation with both the predecessor and the successor lists. For each vertex (oriented fragment), we associate a unique label. We do the same for the edges (the links). Figure 25 illustrates each link case in DGA.

LEGEND

Figure 25 -implementation.

Each squared mathematical formula provides the link case and its reverse (in parentheses). The left and the right lists are the predecessor and the successor lists for all the vertices.

For each v ∈ V , vind(v) = v f if f rag(v) ∈ F f , v r otherwise.
For each link and its reverse, e is their canonical label.

Definition 1.1: labels and functions

Let vind : V → → Σ V and eind : E → → Σ E , where Σ V = Σ F and Σ E = Σ Lcan are the label sets for the vertices and the edges, respectively. The predecessor and the successor lists are defined by two functions:

predl : Σ V → N - V = N - v , ∀ v ∈ V | N - v : 0, d - v → Σ V × Σ E succl : Σ V → N + V = N + v , ∀ v ∈ V | N + v : 0, d + v → Σ V × Σ E For each vertex v ∈ V , N - v : 0, d - v → Σ V × Σ E and N + v : 0, d + v → Σ V × Σ E
give its predecessors and its successors, respectively. The following property holds:

∀ e ∈ E, Φ(e) = (u, v), ∃! n ∈ 0, d + u , ∃! m ∈ 0, d - v s.t. succl(vind(u))(n) = (vind(v), eind(e)) predl(vind(v))(m) = (vind(u), eind(e))
Definition 1.2 provides the logic between the label functions and the reverse properties of the vertices and the edges.

Definition 1.2: DGA reverse labels

Given two fragments a f ∈ F f and a r = a f ∈ F r and their associated vertices

v f and v r (f rag(v f ) = a f , f rag(v r ) = a r ), the labels of v f and v r respect vind(v r ) -vind(v f ) = 1.
For each edge e ∈ E, its label respects eind(e) = eind(e). where P is the memory size of a memory address.

Proof

There is one pointer for the predecessor and the successor lists, and two for each fragment, i.e. 2P(|F| + 1). Then, for each link, for each predecessor and for each successor, the neighbour's label and the edge's label are provided, i.e. 2|L| oct(|F|) + oct 1 2 |L| .

Oriented fragments' successors directed graph (DGS)

DGS corresponds to DGA without the predecessor lists. The property at the end of Definition 1.3 show how to get the predecessors of a vertex (note the d + v ). The reverse of the successors of the reverse are the predecessors of the vertex. Figure 26 illustrates each link case in DGS.

LEGEND

Figure 26 -DGS implementation.

Each squared mathematical formula provides the link case and its reverse (in parentheses). The list bellow is the successor list for all the vertices. The labels u f , u r , v f and v r are the ones of the vertices associated to two forward and reverse fragments. For each link and its reverse, e is their canonical label.

Definition 1.3: DGS labels and functions

Let vind : V → → Σ V and eind : E → → Σ E , where Σ V = Σ F and Σ E = Σ Lcan are the label sets for the vertices and the edges, respectively. The successor list is defined by:

succl : Σ V → N + V = N + v , ∀ v ∈ V | N + v : 0, d + v → Σ V × Σ E For each vertex v ∈ V , N - v : 0, d - v → Σ V × Σ E and N + v : 0, d + v → Σ V × Σ E
give its predecessors and its successors, respectively. The following property holds:

∀ e ∈ E, Φ(e) = (u, v), ∃! n ∈ 0, d + u , ∃! m ∈ 0, d + v s.t. succl(vind(u))(n) = (vind(v), eind(e))
succl(vind(v))(m) = vind(u), eind(e) Definition 1.4 provides the logic between the label functions and the reverse properties of the vertices and the edges.

Definition 1.4: DGS reverse labels

Given two fragments a f ∈ F f and a r = a f ∈ F r and their associated vertices

v f and v r (f rag(v f ) = a f , f rag(v r ) = a r ), the labels of v f and v r respect vind(v r ) -vind(v f ) = 1.
For each edge e ∈ E, its label respects eind(e) = eind(e).

Proposition 1.2 gives the amount of octet DGS consumes.

Proposition 1.2: DGS memory consumption

The memory size M em(DGS) of DGS (in octets) is equal to:

M em(DGS) = P(|F| + 1) + |L| oct(|F|) + oct 1 2 |L|
where P is the memory size of a memory address.

Proof

There is one pointer for the successor lists, and one for each fragment, i.e. P(|F| + 1). Then, for each link, for each successor, the neighbour's label and the edge's label are provided, i.e. |L| oct(|F|) + oct 1 2 |L| .

Forward fragments directed graph (DGF)

DGF is DGA only for the forward fragments. The neighbours of a reverse, are the reverse neighbours of the forward. The vertex label range is divided by two, so the adjacency lists must provide both the label of the neighbours and their orientation.

Figure 27 illustrates each link case in DGF .

LEGEND

Figure 27 -DGF implementation. Each squared mathematical formula provides the link case and its reverse (in parentheses). The left and the right lists are the predecessor and the successor lists for the forward vertices. A vertex and its reverse have the same label (v). For each link and its reverse, e is their canonical label. The orientation set {f, r} is binarised (f = 0, r = 1).

Definition 1.5: DGF labels and functions

Let vind : V → → Σ V f and eind : E → → Σ E , where Σ V f = Σ Fun and Σ E = Σ Lcan are the label sets for the forward vertices and the edges, respectively. The predecessor and the successor lists are defined by two functions:

predl : Σ V → N - V f = N - v , ∀ v ∈ V f | N - v : 0, d - v → Σ V f × {f, r} × Σ E succl : Σ V → N + V f = N + v , ∀ v ∈ V f | N + v : 0, d + v → Σ V f × {f, r} × Σ E For each forward vertex v ∈ V f , N - v : 0, d - v → Σ V f × {f, r} × Σ E and N + v : 0, d + v → Σ V f × {f, r} × Σ E
give its predecessors and its successors, respectively. The following property holds:

∀ e ∈ E, Φ(e) = (u, v), f rag(u) ∈ F f =⇒ ∃! n ∈ 0, d + u s.t. succl(vind(u))(n) = (vind(v), f or(f rag(v)), eind(e)) f rag(u) ∈ F r =⇒ ∃! n ∈ 0, d + u s.t. predl(vind(u))(n) = vind(v), f or(f rag(v)), eind(e) f rag(v) ∈ F f =⇒ ∃! m ∈ 0, d - v s.t. predl(vind(v))(m) = (vind(u), f or(f rag(u)), eind(e)) f rag(v) ∈ F r =⇒ ∃! m ∈ 0, d - v s.t. succl(vind(v))(m) = vind(u), f or(f rag(u)), eind(e)
Definition 1.6 provides the logic between the label functions and the reverse properties of the vertices and the edges.

Definition 1.6: DGF reverse labels

For each vertex v ∈ V , its label respects vind(v) = vind(v). For each edge e ∈ E, its label respects eind(e) = eind(e).

Proposition 1.3 gives the amount of octet DGF consumes.

Proposition 1.3: DGF memory consumption

The memory size M em(DGF ) of DGF (in octets) is equal to:

M em(DGF ) = P(|F| + 2) + |L| oct(|F|) + oct 1 2 |L|
where P is the memory size of a memory address.

Proof

There is one pointer for the predecessor and for the successor lists, and two for each forward fragment, i.e. P(|F| + 2). Then, for each link, for each forward predecessors and successors, the neighbour's label, combined with a binary value (f or r), and the edge's label are provided, i.e. |L| oct(|F|) + oct 1 2 |L| .

1.2 Bidirected graph (BG): oriented walk based Section 1.2.1 describes the first implementation for BG and Section 1.2.2 describes how to transform BGU to DGF . Recall from Definition 1.10 that BG = (V , E, Φ, attre).

Unoriented fragments bidirected graph (BGU)

BGU is a classical multiungraph implementation with neighbour lists. The attributes on the edges are stored in a attribute list. For each vertex (unoriented fragment), we associate a unique index. We do the same for the edges (a links and its reverse). Figure 28 illustrates each link case in BGU .

Definition 1.7: BGU labels and functions

Let vind : V → → Σ V and eind : E → → Σ E , where Σ V = Σ Fun and Σ E = Σ Lcan are respectively the label sets for the vertices and the edges. The neighbour and the edge attribute lists are defined by two functions:

nborl : Σ V → N V = {N v , ∀ v ∈ V | N v : 0, d v → Σ V × Σ E } attrel : Σ E → {f, r} 2
For each vertex v ∈ V , N v : 0, d v → Σ V × Σ E gives its neighbours. The following property holds:

∀ e ∈ E, Φ(e) = {u, v}, ∃! n ∈ 0, d u , ∃! m ∈ 0, d v s.t. nborl(vind(u))(n) = (vind(v), eind(e)) nborl(vind(v))(m) = (vind(u), eind(e))
attrel(eind(e)) = attre(e) Proposition 1.4 gives the amount of octet BGU consumes.

Proposition 1.4: BGU memory consumption

The memory size M em(BGU ) of BGU (in octets) is equal to:

M em(BGU ) = P 1 2 |F| + 1 + |L| oct 1 2 |F| + oct 1 2 |L| + oct(|L|) LEGEND Figure 28 -BGU implementation.
Each squared mathematical formula provides the link case and its reverse (in parentheses). The left and the right lists are the neighbour for all the vertices and the edge attribute list, respectively. For each fragment and its reverse v is their canonical index. For each link and its reverse, e is their canonical index. Eattr provides the orientations associated with the vertices, in the lexicographical order of their label (here u < v).

where P is the memory size of a memory address.

Proof

There is one pointer for each vertex and one for the edge attribute list, i.e. P 1 2 |F| + 1 . Each vertex is associated with its neighbours and the edge's label. Each edge is stored twice (|L|). So it gives |L| oct 1 2 |F| + oct 1 2 |L| . Finally, Eattr provides for each edge two binary values (one bit per orientation), i.e. oct(|L|).

Transformation to DGF

It is possible to transform BGU to DGF . For each u ∈ Σ V , for each (v, e) ∈ Σ V ×Σ E in cod nborl(u) (cod is the codomain of function nborl) such that u < v implying attrel(e) = (uor, vor):

-if uor = f , then append (v, vor, e) to the successor list N + u ; -else if uor = r, then append (v, vor, e) to the predecessor list N - u ; -if vor = f , then append (u, uor, e) to the predecessor list N - v ; -else if vor = r, then append (u, uor, e) to the successor list N + v .

Undirected graph (UG): tail-head fragments based

Section 1.3.1 describes the first implementation for U G and Section 1.3.2 describes how to transform U GA to DGS. Recall from Definition 1.13 that U G = (V , E F , E L , Φ L ).

All oriented fragments undirected graph (UGA)

U GA is a classical multiungraph implementation with neighbour lists. Each fragment is associated with two vertices corresponding to the two extremities (5' and 3'). We provide a unique label for each vertex. There is one neighbour list for the fragment-edges and another one for the link-edges. Figure 29 illustrates each link case in U GA.

Definition 1.8: UGA labels and functions

Let vind : V → → Σ V , and eind : E L → → Σ E a , where Σ V = Σ F and Σ E = Σ Lcan are the label sets for the vertices and the link-edges, respectively. The fragmentneighbour function is denoted by nborlf : Σ V → → Σ V . For each vertex v ∈ V , nborlf (v) gives the other extremity of the fragment. The link-neighbour list is defined by:

nborll : Σ V → N L V = N L v , ∀ v ∈ V | N L v : 0, d L v → Σ V × Σ E For each vertex v ∈ V , N L v : 0, d L v → Σ V × Σ E provides its link-neighbours.
The following property holds:

∀ e ∈ E L , Φ L (e) = {u, v}, ∃! n ∈ 0, d L u , ∃! m ∈ 0, d L v s.t. nborll(vind(u))(n) = (vind(v), eind(e)) nborll(vind(v))(m) = (vind(u), eind(e)) a
We do not provide edge labels for fragment-edges. Definition 1.9 provides the logic between the label functions and the reverse properties of the vertices. Definition 1.9: UGA reverse labels Given a fragments a f ∈ F f and its reverse a r = a f ∈ F r and their associated vertices

{v t , v h } ∈ E F (f rag(v t , v h ) = a f , f rag(v h , v t ) = a r ), the label of v t LEGEND Figure 29 -UGA implementation.
Each squared mathematical formula provides the link case and its reverse (in parentheses). The violet list is not stored in memory because it does not depend of the link case. It corresponds to nborlf function. The N L

V green list provides the link-neighbours. Vertices v t and v h correspond to the tail and the head of a fragment. For each link and its reverse, e is their canonical index. and v h respect vind(v h ) -vind(v t ) = 1.

Proposition 1.5 gives the amount of octet U GA consumes.

Proposition 1.5: UGA memory consumption

The memory size M em(U GA) of U GA (in octets) is equal to:

M em(U GA) = P(|F| + 1) + |L| oct(|F|) + oct 1 2 |L|
where P is the memory size of a memory address.

Proof

There is one pointer for the link-neighbour lists and one for each fragment extremity, i.e. P(|F| + 1). Then, for each link, for each extremity, the neighbour's label and the edge's label are provided, i.e. |L| oct(|F|) + oct 1 2 |L| .

Transformation to DGS

U GA can easily be adapted to obtain DGS. For each v ∈ Σ V , for each n ∈ 0, d L v where nborll(v)(n) = (w, e), replace w by nborlf (w). In a valid path, a link-edge is followed by a fragment-edge, which is unique for each vertex (see Definition 1.15). Finally, to exactly match DGS, the order of the index for the extremities is reversed, i.e. vind(v t ) -vind(v h ) = 1.

Fragment graph map

Figure 30 -Graph structures and their implementations. Blue boxes: they correspond to graph structures. They are described by the type of graph they represent, the theoritical number vertices and edges. Green boxes: they correspond to discussed graph implementations. After their name, the second line gives the adjacency list(s) definition and the last line gives the length of the adjacency list(s) (i.e. the number of vertices for which the neighbours are given). In this section, we give elementary operation algorithms for DGS, DGF and BGU graph implementations. To respect the computing notations, N V , N - V and N + V are the neighbour, predecessor and successor lists, respectively, codomains of the functions nborl, predl and succl. Similarly, for any vertex v we denote by N v , N - v and N + v the neighbour, predecessor and successor lists of v, respectively. Iterating over one of these lists is equivalent to increasing an integer k from 0 to the length of the list minus one and peek the k th element.

DG: oriented Fragments based

Directed

Graph structure chronology

Classical digraph

Subfunctions

The following subfunction algorithms are often subparts of the next algorithms. Algorithm 1 is usefull for DGS. Algorithms 2 and 3 are usefull when deleting vertices or edges.

Iterating over the predecessors

Given a fragment, we aim to return the fragments linking it upstream associated with the canonical link labels.

Iterating over the successors

Given a fragment, we aim to return the fragments linking it downstream associated with the canonical link labels.

Adding a vertex

The following algorithms allocate a place for new vertices by providing their label. N V .append(empty_list( ))

3:

return |N V | -1

Adding an edge

The following algorithms append new edges by providing their label. ind_edges is the number of edges. It is always even. ind_edges is the number of edges. It is always even.

3: if uor = f then u f → v f or u f → v r 4: N + u .append(v, vor, ind_edges) 5: else u f ← v f or u f ← v r 6: N - u .append(v, 1 -vor, ind_edges) 7: if vor = f then u f → v f or u r → v f 8: N - v .append(u, uor, ind_edges) 9: else u f ← v f or u r ← v f 10: N + v .append(u,

Deleting a vertex

The following algorithms delete a vertex (and its reverse where appropriate) and all their edges.

Algorithm 16: Delete a vertex for DGS Require: A vertex label v such that v | 2. Ensure: The edges containing the vertex or its reverse are deleted from the neighbour lists. The length of N + V is reduced by two.

if v = N + V -2 then 19:
It is the last, just pop it and its reverse 20:

N + V .pop( ) 21: N + V .pop( ) 22:
else 23:

Replace it and its reverse by the last and its reverse 24: Delete the whole vertex 23:

N + v_rev ← N + V .pop( ) 25: N + v ← N + V .
if v = N + V -1 then 24:
It is the last, just pop it and its reverse 25:

N - V .pop( ) 26: N + V .pop( ) 27:
else 28:

Replace it and its reverse by the last and its reverse 29:

N - v ← N - V .pop( ) 30: N + v ← N + V .pop( ) 31:
Update the vertex index 32:

for (u, uor, e) ∈ N - v do 33:

if uor = f then 34:

adj_ind ← get_index(N + u , e, 1)

35:

N + u [adj_ind] ← (v, f , e) 36:
else 37:

adj_ind ← get_index(N - u , e, 1)

38:

N - u [adj_ind] ← (v, r, e) 39:
for (w, wor, e) ∈ N + v do 40:

if wor = f then 41: adj_ind ← get_index(N - w , e , 1) 42:
N The last (the greatest) identifier becomes a fewer and breaks the identifier order 20:

- w [adj_ind] ← (v, f , e) 43:
Eattr[e] ← (1 -Eattr[e][1], 1 -Eattr[e][0])

Deleting an edge

The following algorithms delete an edge (and its reverse where appropriate).

Algorithm 19: Delete an edge for DGS

Require: Two vertex labels u and v and a edge label. Ensure: The edges are deleted for the two vertices.

1: function dgs_delete_edge(u, v, e)

2:

Remove v from u's successors 3:

adj_ind ← get_index(N + u , e, 1)

4:

delete_list_entry(N + u , adj_ind)

5:

Remove u from v's predecessors. adj_ind ← get_index(list_succs, e, 2)

8:

delete_list_entry(list_succs, adj_ind)

9:

Remove u from v's predecessors. Remove v from u's neighbours.

3:

adj_ind ← get_index(N u , e, 1)

4:

delete_list_entry(N u , adj_ind)

5:

Remove u from v's neighbours.

6:

adj_ind ← get_index(N v , e, 1)

7:

delete_list_entry(N v , adj_ind)

8:

card_edges ← card_edges -1

Time costs

In this section we provide the time costs for each of the algorithms in Section 2, in the best, worst and average cases. In the average cases, we consider all the events to be equiprobable. Table 4 details the cost for basic operations. In order to easily compare the algorithms' theoritical speed, we assume there is an even integer k such that each fragment has k upstream and k downstream fragments linking it. It holds if |F| >> k and |L| >> k.

Table 5 gives the algorithmic costs of subfunctions in Algorithms 1 and 2. Table 6 provides the algorithmic costs for Algorithms 4 to 9. Table 7 gives the

|F |-1 |F | k 2 + 36 |F |-2 3 |F | k + 1-|F | |F | DGF 24k + 4 20k 2 + 28k + 6 10 |F |-1 |F | k 2 + 28 |F |-4 7 |F | k + 4 |F |-1 |F | BGU 22k + 3 40k 2 + 38k + 7 20 |F |-1 |F | k 2 + 24 |F |-1 2 |F | k + 11 2 |F |-14 11 |F |

Conclusions and perspectives

Graphs are mathematical structures that are useful to store links between two fragments. Two main stages of the fragment assembly employ graphs to handle the input data or to output the result. In the read assembly stage, especially in the OLC approach, the fragments correspond to the reads and the links are the

FRAGMENT GRAPHS III

Table 10 -Comparison between DGS, DGS and BGU. Each implementation is associated with starred scores equal to 1 (bad, in red), 2 (medium, in orange) or 3 (good, in green). Two implementations can have the same score if their memory or algorithmic costs are near-equal, and if the third one's cost is far away from them. For each elementary graph operation, for each of the three cases (best, worst, average), the stars correspond to a score overlaps between two reads. In the scaffolding stage, the fragments represent the contigs and links can come from paired-end read or long read alignments against the contigs.

DGS DGF BGU

Memory consumption

Iterating over the predecessors

The literature describes three graph structures: a digraph, a bigraph and an ungraph. While the first one highlights visually the double-strand sequencing, its weakness is that it requires two vertices for each read [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF]. [START_REF] Myers | Toward simplifying and accurately formulating fragment assembly[END_REF] is the first to employ a bigraph to store overlaps. The key idea is to aggregate the two orientations of a read into only one vertex, a link and its reverse into only one edge. Finally, the ungraph structure associates one vertex for each of the two extremities of a fragment to simplify the graph traversal [START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF]. The ungraph handles two sets of edges. A fragment-edge connects the two vertices corresponding to the two extremities of a fragment (the tail and the head). Two vertices representing the extremities of a fragment are connected by a fragment-edge. The second one is the multiset of link-edges. A link-edge represents both a link and its reverse.

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

We have retained three implementations: DGS and DGF for the multidigraph and BGU for the multibigraph. After we have compared their memory consumption, we theoretically measure the cost function for each of their algorithm on elementary operations, such as iterating over the neighbours, adding or deleting a vertex or an edge. We come to the conclusion that if memory is the critical issue, then the BGU implementation should be preferred. Otherwise, we recommend DGS and DGF for iterating over the neighbours and for deleting an edge. BGU is preferable for adding vertices or edges while DGF is the best for deleting vertices. To conclude, DGF proves to be well-balanced and ideally tailored for dynamic graph operations. DGF is available in a Python3 package named revsymg1 and is easily installable via PyPI. In this chapter, we focus on the scaffolding problem specific to chloroplast genomes. We first recall the specificities of these genomes as described in the introduction (Chapter I, Section 2), highlight the state-of-the-art in chloroplast genome assembly (Chapter II, Section 3) and present our approach.

1 Introduction

Chloroplast genome specificities

Chloroplasts are plants' organelles derived from the integration of a cyanobacterium in an eukaryotic host. They conduct photosynthesis, a process to convert light energy into chemical energy. Over the evolution time, the chloroplast genome has reduced in length and loosed in terms of complexity [START_REF] Xiao-Ming | Inferring the evolutionary mechanism of the chloroplast genome size by comparing wholechloroplast genome sequences in seed plants[END_REF]. 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. The circular genome is partitionned in four regions: two identical (or highly similar) nucleotide subsequences, separated by two single-copies (Long Single-Copy, LSC, and Short Single-Copy, SSC) [START_REF] Bock | Genomics of Chloroplasts and Mitochondria[END_REF][START_REF] Palmer | Comparative Organization of Chloroplast Genomes[END_REF]. 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 6 illustrates the common chloroplast genome architectures. Furthermore, each chloroplast has multiple copies of its genome, and the molecular forms of the copies differ (multimeric genome forms, [START_REF] Bendich | Circular chloroplast chromosomes: The grand illusion[END_REF]). This phenomenon can be induced by flip-flop inversion: one subsequence is reversecomplemented (reversed) during the DNA replication. This inversion is provoked by the existence of facing inverted repeats on either side of the reversed subsequence.

State-of-the-art

Although there are chloroplast genome assemblers and scaffolders, they do not fully exploit the chloroplast genome's specificities. Some of them are pipelines of generic methods applied on cleaned input dataset [START_REF] Ankenbrand | chloroExtractor: Extraction and assembly of the chloroplast genome from whole genome shotgun data[END_REF], or based on locally approaches as seed-and-extend algorithms [START_REF] Coissac | From barcodes to genomes: Extending the concept of DNA barcoding[END_REF][START_REF] Dierckxsens | NOVOPlasty: De novo assembly of organelle genomes from whole genome data[END_REF]. The authors of GetOrganelle [START_REF] Jin | GetOrganelle: A fast and versatile toolkit for accurate de novo assembly of organelle genomes[END_REF] statistically compute the contigs' multiplicities in order to they are the closest as possible to the mapping coverage by the reads.

Concerning the handling of multimeric forms, GetOrganelle returns several solutions and explores in post-process the corresponding architectures. In [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF] the flip-flop inversion breakpoints are detected in a post-scaffoldingprocess, and the optimality of the new solutions is proven in the case of the use of distances.

Our approach

We raise the following two questions: (i) Is it possible to scaffold chloroplast genomes without using any distances? (ii) How to include the multimerism into the scaffolding problem formulation?

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

We first describe the input data for our approach and provide mathematical definitions (Section 2). Based on the above three assumptions, we propose a new formulation for scaffolding chloroplast genomes without requiring any distances (Section 3). Our approach is region-driven and focuses on retrieving the repeats first. We model the optimisation problem on a digraph where we apply several ILP (Sections 4 and 5). We then detail how we combine the ILP solutions (Section 6). The ILP's solutions and the digraph correspond to an oriented contig order representing one genome's form (Section 7). We partition this order 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 N P-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).

Input data and notation

Set of contigs C

Contigs are words in the nucleotides DNA alphabet Σ nuc + . A contig can occur in the genome up to an integer called multiplicity. Function mult : C → 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 R ≥0 given by the function wex. The weight indicates the number of times a contig aligns with chloroplast sequences from a given set (from related or unrelated species). Finally, one contig in this set is defined as the starter (s) that must uniquely participate in the genome (mult(s) = 1). The starter is a contig whose sequence matches a sequence shared by most chloroplast genomes in a single-copy. Table 11a gives an example of contig set.

Set of links L

Each link is an ordered pair of oriented contigs. We denote by L ⊂ (C × {f, r}) 2 the link set1 . The nature of the double-strands DNA requires that ∀ (c, d) ∈ L, (d, c) ∈ L, where c and d denote the oriented contigs c and d in their reverse orientation, respectively (note that c = c). The links between two oriented contigs c and d are valid for all occurrences of c and d respecting the same orientations. Table 11b gives an example of link set.

Mathematically defining genomic regions

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

Definition 2.1: Sequence of oriented contigs

Let SOC = (c 0 , c 1 , . . . , c n-1 ) be a sequence of oriented contigs:

-∀ i ∈ 0, n -1 , (c i , c i+1 ) ∈ L; -∀ c ∈ C, c i ∈SOC|c i =c 1 ≤ mult(c).
Based on the biological knowledge, we address the dedicated chloroplast scaffolding problem as a region-driven scaffolding, such that specific regions fit into a circular structure. We identify three types of regions to scaffold: the Inverted Repeat (IR), the Direct Repeat (DR) and the Single-Copy (SC).

Definition 2.2: Region

A region r = (c 0 , c 1 , . . . , c n-1 ) is a sequence of oriented contigs. Each region is oriented. Let r = c n-1 , . . . , c 1 , 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, r also respects Definition 2.1.

Definition 2.3: Direct repeat (DR)

A DR is a couple of regions (dr i , dr j ) where dr i = dr j .

Definition 2.4: Inverted repeat (IR)

An IR is a couple of regions (ir k , ir l ) where ir l = ir k .

Definition 2.5: Repeat

A repeat is the generic term to denote either DR or IR. The length replen(R) of a repeat R = (r i , r j ) equals the length of r i (or r j , as

|r i | = |r j |).

Definition 2.6: Single-copy region (SC)

A Single-Copy (SC) is a region that is not part of a repeat.

Definition 2.7: Region weight

The weight rwex(r) of a region r is defined as rwex(r) = c∈r wex(c).

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

Definition 2.8: Sequence of oriented regions

Consider a sequence SOR = (r 0 , r 1 , . . . , r n-1 ):

-∀ i ∈ 0, n -1 , (r i [|r i | -1] a , r i+1 [0]) ∈ L; -∀ c ∈ C, r∈SOR c i ∈r|c i =c 1 ≤ mult(c).
a The notation x[i] denotes the i th element in x.

Chloroplast scaffolding problem formulations

Solving the repeats is the most challenging task. A formulation that does not restrict the occurrences can lead to misassemblies where the results are longer than 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 only the longest repeats of the same size.

However, in the case of repeat degeneration (e.g. two subsequences inside of the two regions of an identified repeat differ) finding the longest repeats is not an appropriate model. Figure 31 illustrates the impact of degenerations on quadripartite structures. Indeed, in Figure 31b, we cannot guaranty to find both IR 1 and IR 2 , but perhaps only the longest. We address this issue by considering the cumulative length of the longest repeats by constraining their order if they are of the same type. Note that some IR losses have been reported in the chloroplast genomes of green algae [START_REF] Turmel | Divergent copies of the large inverted repeat in the chloroplast genomes of ulvophycean green algae[END_REF]. highlight the impact of a their degeneration on their structures. We give the bellow region orders according to the LSC arrow's direction. In (b), the degeneration results in two IRs: IR 1 = (i, j) and IR 2 = (k, l), such that i is before k and l precedes j. In (d) it results in two DRs: DR 1 = (i, j) and IR 2 = (k, l), such that i is before k and j precedes l.

Definition 3.1: Chloroplast scaffolding problem CHSP

Given a set of contigs with their multiplicities and their weights, a starting contig and a link set. The aim is to obtain a circular sequence of oriented regions maximising the cumulative length of the longest repeats, joined by the maximum-weighted single-copies.

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

DRP and IRP constrain the number of occurrences to the structure of pairs of repetitions. Indeed, each repeat type defines a valid repeat structure. The problems consist in finding the longest (valid) repeated nucleotide sequences, and maximising their number.

Definition 3.2: Chloroplast direct repeat scaffolding problem DRP

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

-it maximises the cumulative length of the longest DR, joined by regions of any kind;

-for any couple of DR (i, j) and (k, l) found in SOR such that their respective positions σ respect σ(i) < σ(j), σ(k) < σ(l), and σ(i) < σ(k), then:

-

σ(i), σ(j) ∩ σ(k), σ(l) = ∅;
-or σ(k), σ(j) ⊂ σ(i), σ(l) .

Definition 3.3: Chloroplast inverted repeat scaffolding problem IRP

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

-it maximises the cumulative length of the longest IR, joined by regions of any kind;

-for any couple of IR (i, j) and (k, l) found in SOR such that that their respective positions σ respect σ(i) < σ(j), σ(k) < σ(l) and σ(i) < σ(k), then:

-

σ(i), σ(j) ∩ σ(k), σ(l) = ∅; -or σ(k), σ(l) ⊂ σ(i), σ(j) .
Figure 32 provides examples of valid oriented contig positioning for each common genome structure (Figure 6). Although Figure 36 illustrates the authorised and forbidden positions for the latter defined repeated fragments, it is also applicable for the DRP and IRP regions' position cases.

Definition 3.4: Chloroplast single-copy scaffolding problem SCP

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

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

Definition 3.5: Chloroplast scaffolding problem succession

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

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

Definition 3.6: Hierarchical problem succession

The form of each solution of CHSP satisfies one of the two problem successions: DRP-IRP-SCP (h 1 ) and IRP-DRP-SCP (h 2 ).

The next question is how to prioritise DRP and IRP? We propose resolving the order by comparing the scores defined in Section 5.5: if DRP score is better than this of IRP, then the retained succession will be DRP-IRP-SCP, otherwise it will be IRP-DRP-SCP. In the equality case, we discriminate at a further step of the hierarchical successions. The process is detailed in Section 6.

Finally, each hierarchical problem succession produces a circular sequence of oriented regions. From the obtained sequence it is possible to extract a set of ordered pairs of oriented regions. This procedure allows the building of several circular sequences of oriented regions of the same length. Each of them represents one possible chloroplast genome form. This all-equivalent-form process is described in Section 8.

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.

Graph structure

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

Definition 4.1: Multiplied Doubled Contig Graph -M DCG

Given a set of contigs C, their multiplicities and the link set L, the multiplied doubled contig graph M DCG = (V, E, vwex) is defined such that a :

-V = {v f,0 , . . . , v f,n-1 , v r,0 , . . . , v r,n-1 | c ∈ C, n = mult(c)} is the set of all
the forward an reverse occurrences of all the contigs (|V | = 2 c∈C mult(c)).

The vertices are associated with four functions:

contig : V → → C provides the contig associated with a vertex;

vor : V → → {f, r} provides the orientation of the contig;

vocc : V → N >0 provides the occurrence number of the contig;

-vwex : V → R ≥0 provides the weight of each vertex such that ∀ v ∈ V, vwex(v) = wex(contig(v))
.

-E = {(u, v) ∈ V 2 | ((contig(u), vor(u)), (contig(v), vor(v))) ∈ L} is the set of multiplied links (|E| = (c,d)∈L mult(c)mult(d)).
a M DCG is a digraph, not a multidigraph as DG.

Figure 33 illustrates the M DCG representing the example data given in Table 11. Each vertex is associated with an occurrence of an oriented contig, and each contig is represented by an even number of vertices. For example, vertex labelled c r,1 means that it comes from contig c = contig(c r,1 ), in its reverse orientation (vor(c r,1 ) = r), and in its second occurrence (vocc(c r,1 ) = 1). The colours are the same as the ones in Figure 6a 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.

Repeated fragment sets

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

Definition 4.2: Repeated fragment

A repeated fragment is an unordered pair of vertices such that one of the corresponding oriented contig belongs to the first region of a repeat, while the other belongs to the second region. The vertices are associated with the same contig but their occurrences differ, i.e. for each repeat

(r i , r j ), ∃ u, v ∈ V, c ∈ C where contig(u) = contig(v) = c and vocc(u) = vocc(v) such that (c, vor(u)) ∈ r i and (c, vor(v)) ∈ r j .
For example, (c f,0 , c r,1 ) is a repeated fragment for the IR in Figure 33. We then precise the set of repeated fragments for each repeat type. Denote by R = {c ∈ C | mult(c) > 1} the set of contigs candidate to be part of repeats. For the sake of clarity, for each vertex v ∈ V , we note

ctg v = contig(v), or v = vor(v), occ v = vocc(v), wex v = vwex(v) and mult v = mult(contig(v)).
We also assume there is an arbitrary strict total order on

C, i.e. ∀ c, d ∈ C, c = d ⇐⇒ c < d c > d.

Definition 4.3: Set of direct fragments DirF

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

DirF = c∈R                    (i, j) ∈ V 2 s.t. ctg i = ctg j = c ∧ or i = or j ∈ {f, r} ∧ occ i = occ j -1 = 2k 0 ≤ k < mult(c) 2                   

Definition 4.4: Set of inverted fragments InvF

An inverted fragment (u, v) ∈ V 2 is a repeated fragment, such that the orienta-

(u f,0 , u f,1 ) ∈ DirF so dirf rag (u f,0 ) = dirf rag (u f,1 ) = (u f,0 , u f,1 ). |P DirF | = 12
, and, e.g.

((u f,0 , u f,1 ), (u r,2 , u r,3 )) ∈ P DirF . (b) |InvF | = 3, and, e.g. (u f,2 , u r,3 ) ∈ InvF so invf rag (u f,2 ) = invf rag (u r,3 ) = (u f,2 , u r,3 ). |P InvF | = 3, and, e.g. ((u f,2 , u r,3 ), (v f,0 , v r,1 )) ∈ P InvF .
tions of u and v differ.

InvF = c∈R                    (i, j) ∈ V 2 s.t. ctg i = ctg j = c ∧ or i = f ∧ or j = r ∧ occ i = occ j -1 = 2k 0 ≤ k < mult(c) 2                   
Figure 34 illustrates DirF and InvF sets. In addition, we add two functions to retrieve the repeated fragments from the vertices in Θ(1): dirf rag : V ⊂ V → DirF and invf rag : V ⊂ V → InvF (abstracted with the repf rag function, c.f. Appendix 1 for their definitions).

Furthermore, Definitions 3.2 and 3.3 contrain the order between the repeated fragments. Hence, they respectively require to compare pairs of direct/inverted fragments, that must be defined: Definition 4.5: Set of pairs of direct fragments

P DirF =          (i, j), (k, l) ∈ DirF 2 s.t. ctg j < ctg k ∨ ctg j = ctg k ∧ occ j < occ k         
Definition 4.6: Set of pairs of inverted fragments

P InvF =          (i, j), (k, l) ∈ InvF 2 s.t. ctg j < ctg k ∨ ctg j = ctg k ∧ occ j < occ k         
Figure 34 illustrates P DirF and P InvF sets. The constraints defining DirF , InvF , ADirF , AInvF , P DirF and P InvF are explained in details in Appendix 2 where we proof they are the smallest sets enabling to find all the distinct solutions.

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

Definition 4.7: Set of adjacent repeated fragments

An adjacent repeated fragment is an edge (u, v) ∈ E such that u and v participate in two distinct repeated fragments.

Definition 4.8: Set of adjacent direct fragments

An adjacent direct fragment is an edge between two direct fragments. Let

((c, f ), (d, f )) ∈ 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

∀ (u, v) ∈ E, v, u ∈ E. (a) diradj (u f,0 , v f,0 ) = (u f,1 , v f,1 ); (b) invadj (u f,0 , v f,0 ) = (v r,1 , u r,1 ).
ADirF be the set of adjacent direct fragments:

ADirF =                            (u, v) ∈ E s.t. ctg u = ctg v ∧ occ u = 2k 0 ≤ k < mult u 2 ∧ occ v = 2k 0 ≤ k < mult v 2                                                   (u, v) ∈ E s.t. ctg u = ctg v ∧ (or u = f ∨ or v = f ) ∧ occ u = 2k ∧ occ v = 2k 0 ≤ k < k < mult u 2                       
Definition 4.9: Set of adjacent inverted fragments An adjacent inverted fragment is an edge between two inverted fragments. Let

Circuit constraints

The following set of constraints defines a valid circuit of oriented contig in M DCG, and is defined with a flow formulation as in [START_REF] Andonov | Complete assembly of circular and chloroplast genomes based on global optimization[END_REF] instead of using Miller-Tucker-Zemlin constraints to avoid cycles [START_REF] Miller | Integer Programming Formulation of Traveling Salesman Problems[END_REF].

Binary variables

-x e encodes if the edge e ∈ E participates in the circuit.

Continuous variables

-i v ∈ [0, 1] encodes if the vertex v ∈ V \ s, s participates in the circuit.
Although it is a continuous variable, it acts as a binary one as proven in [START_REF] François | Global Optimization for Scaffolding and Completing Genome Assemblies[END_REF].

-f e ∈ R ≥0 is the positive flow on the participating edge e ∈ 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 ≤ f e ≤ Mx e ∀ e ∈ E (C1) v∈N + s x sv = v∈N - s x vs = 1 (C2) v∈N + s f sv = 1 (C3) x vs = 0 ∀ v ∈ N - s (C4) x sv = 0 ∀ v ∈ N + s (C5) ∀ v ∈ V \ s, s : i v + i v ≤ 1 (C6) u∈N - v x uv ≤ i v ≤ w∈N + v x vw (C7) w∈N + v f vw - u∈N - v f uv = i v (C8) x e ∈ {0, 1} ∀ e ∈ E i v ∈ [0, 1] ∀ v ∈ V \ s, s f e ∈ R ≥0 ∀ e ∈ E

Repeated regions constraints

The following constraints are general to define ILPs for DRP and IRP. Definitions 3.2 and 3.3 define the repeated regions according to the positions of the oriented contig in them. It follows that some order of the vertices in the pairs of repeated fragments are allowed, and some others are forbidden. We decide to write the constraints for the forbidden cases because they are fewer than the allowed ones. To model the forbidden orders between 4 vertices, we compare the positions between two. Specifically for IRP, modelling the forbidden orders echoes the approach for the RNA folding problem [START_REF] Gusfield | The RNA-Folding Problem[END_REF], except that the positions of the RNA's nucleotides are known.

According to Definitions 3.2 and 3.3, and given P DirF and P InvF (Definitions 4.5 and 4.6), denote by F orbidDR and F orbidIR the sets of forbidden quartet vertices for the DRs and IRs, respectively:

F orbidDR = (i, k, l, j) (k, i, j, l) i, j ∈ p, i = j ∧ k, l ∈ q, k = l ∀ (p, q) ∈ P DirF F orbidIR = (i, k, j, l) (k, i, l, j) i, j ∈ p, i = j ∧ k, l ∈ q, k = l ∀ (p, q) ∈ P InvF
Figure 36 illustrates the authorised and forbidden positions for DRP and IRP.

To know if we are in the forbidden cases described in these two sets, we propose to compare the vertices two-by-two. Denote by AlphaDR and AlphaIR the sets containing the couples of vertices to be compared to determine the forbidden cases 106 respectively associated with F orbidDR and F orbidIR sets, such that:

AlphaDR =      (i, j), (k, l), (i, k), (j, l), (i, l), (j, k) s.t. ((i, j), (k, l)) ∈ P DirF      AlphaIR = (i, k), (i, l), (j, k), (j, l) s.t. ((i, j), (k, l)) ∈ P InvF
In the following, the sets of repeated fragments and these for the forbidden orders are abstracted to generalise DRP and IRP. Table 12 gives the correspondence of the sets depending on the problem to solve.

Binary variables

-m ij encodes if the repeated fragment (i, j) ∈ RepF is a part of a repeat.

-isadj e encodes if two repeated fragments connected by the edge e ∈ ARepF (and repadj(e) ∈ E) are adjacent in the circuit.

-f orbid ijkl encodes whether we are in the forbidden vertices order (i, j, k, l) ∈ F orbid.

-α uv encodes whether the vertex u is before the vertex v in the circuit. Since

α uv = 1 -α vu , even if (v, u) / ∈ Alpha, for clarity we write α vu instead of 1 -α uv .
Continuous variables

-i v ∈ [0, 1] encodes if the vertex v ∈ V \ s,
s participates in the circuit, and acts as a binary variable.

-f e ∈ R ≥0 is the positive flow on the participating edge e ∈ E in the circuit (zero otherwise). We use the exiting flow to define the position pos

(v) of a vertex v ∈ V , i.e. pos(v) = w∈N + v f w .
The vertices of participating repeated fragments must be in the circuit (Constraints C9 and C10). Constraints C11 to C13 implement with linear constraints the α uv definition. Constraints C14 to C16 implement the f orbid ijkl definition. Constraints C17 to C20 implement the isadj e definition.

CRepeat constraints

Add the set of constraints CCircuit

∀ (i, j) ∈ RepF : m ij ≤ i i (C9) m ij ≤ i j (C10) ∀ (u, v) ∈ Alpha : pos(v) -pos(u) ≤ Mα uv (C11) pos(u) -pos(v) ≤ M(1 -α uv ) (C12) pos(u) + pos(v) ≥ α uv (C13) ∀ (i, j, k, l) ∈ F orbid : 3f orbid ijkl ≤ α ij + α jk + α kl (C14) 2 + f orbid ijkl ≥ α ij + α jk + α kl (C15)
∀ (p, q) ∈ P RepF :

m p + m q ≤ 2 - (i,j,k,l) ∈F orbid s.t. (p,q) f orbid ijkl (C16) ∀ (u, v) ∈ ARepF : isadj uv ≤ x uv (C17) isadj uv ≤ x repadj(u,v) (C18) isadj uv ≤ m repf rag(u) (C19) isadj uv ≤ m repf rag(v) (C20) m p ∈ {0, 1} ∀ p ∈ RepF isadj e ∈ {0, 1} ∀ e ∈ ARepF f orbid ijkl ∈ {0, 1} ∀ (i, j, k, l) ∈ F orbid α uv ∈ {0, 1} ∀ (u, v) ∈ Alpha

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.

CF ixRegions constraints

∀ (u, v) ∈ DirF * ∪ InvF * : i u = 1 (C21) i v = 1 (C22) ∀ (u, v) ∈ ADirF * ∪ AInvF * : x uv = 1 (C23) x diradj(u,v) = 1 ∀ (u, v) ∈ ADirF * (C24) x invadj(u,v) = 1 ∀ (u, v) ∈ AInvF * (C25)

Speed-up constraints

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

ConsOcc =      (u, v) ∈ V 2 ctg u = ctg v ∧ or u = or v = f ∧ occ u = occ v -1     
Also, denote by ConsRepF the set of consecutive repeated fragments, such that:

ConsRepF =                ((i, j), (k, l)) ∈ RepF 2 s.t. ctg i = ctg j = ctg k = ctg l ∧ or i = or k ∧ or j = or l ∧ occ i = occ k -2 ∧ occ j = occ l -2                i v + i v ≤ i u + i u ∀ (u, v) ∈ ConsOcc (C26) m q ≤ m p ∀ (p, q) ∈ ConsRepF (C27)

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 , P RepF , ARepF , Alpha, F orbid and ConsRepF according to the repeats the problems scaffold. We aim to maximise the cumulative length (the number of repeated fragments) of the longest repeats (maximising the number of their adjacent repeated fragments). Note that the objective value corresponds to r∈Repeats (2replen(r) -1), where Repeats is the set of repeats.

DRP

Hierarchical problem succession

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

Definition 6.1: Hierarchical problem succession solutions

Denote by h 1 , h 2 the two hierarchical problem successions DRP-IRP-SCP and IRP-DRP-SCP. For each h ∈ {h 1 , h 2 }, denote by Z * h ∈ R ≥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 = s ∀ k ∈ 0, 2 , Z * s [k] = max h∈{h 1 ,h 2 } Z * h [k]
Note that 0 ≤ |S| ≤ 2, and Definition 6.1 is stable for any problem with an objective value equal to zero. For example, Z * h 1 [START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF] = 0 means that there is no inverted repeat. For an easier interpretation of the architecture, we adopt a problem code combination, summarised in Table 13.

Table 13 -Problem code combinations

Z * h h 1 h 2 [0] [1] [2] DRP-IRP-SCP IRP-DRP-SCP 0 0 -sc sc > 0 0 -dr-sc ir-sc > 0 > 0 -dr-ir-sc ir-dr-sc
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 ∈ 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 , . . . , r m+n-1 ) contains the forward regions, i.e. it is a (m+n)-uplet of m + n oriented contig sequences.

-SOR = ((rid 0 , ror 0 ), . . . , (rid 2m+n-1 , ror 2m+n-1 )) ∈ (N × {f, r}) 2m+n is a linearised circular sequence of oriented regions. For each i ∈ 0, 2m + n , if

ror i = f then SOR[i] represents the forward region COR[rid i ], else if ror i = r then SOR[i] represents the reverse COR[rid i ].
Figure 37 illustrates the regions extracted by Algorithm 25 in the toy example.

The key idea of such an extraction algorithm is to start from the starting vertex s and walk over the chosen edges. During the walk, for each vertex we need to identify the type of its region, and then add the vertex to the corresponding region (to the corresponding region identifier rid r ). Algorithm 22 aims to determine the The illustrated CHSP solution consists of the red circuit (a f,0 , c r,1 , b f,0 , d f,0 , b r,1 , c f,0 ) and of the two choosen inverted fragments InvF * = {(c f,0 , c r,1 ), (b f,0 , b r,1 )} (vertex pairs linked by the green dashed lines). Algorithm 25 returns four regions: m = 1 inverted repeat (the two green arrows pointing downwards) and n = 2 single-copies (the purple and the red arrows). COR = ((a f ), (c r , b f ), (d f )) and SOR = (0 f , 1 f , 2 f , 1 r ), where x y = (x, y) for clarity.

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

To build the repeated regions, we use a First In First Out (FIFO) for the DRs, and a Last In First Out (LIFO) for the IRs. Each queue rep_queue, is associated with three methods: rep_queue.put(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 22 is in Θ(1). Algorithm 23 is in O(|SC s |), where |SC s | is the number of contigs composing the single-copy region that contains the starting vertex. Algorithm 24 is in Θ(1), and so Algorithm 25 is in O(|V | + |E|).

Algorithm 22: Get the region type for a given vertex

Require: A vertex v, the chosen direct fragments DirF * and the chosen inverted fragments InvF * . Ensure: Returns the region type of the vertex. init_v ← initial_vertex() 37 results in the graph visualised in Figure 33 will give a region graph that shows an LSC-IR-SSC-IR architecture as the one given in Figure 6a. Two oriented region sequences (genome forms) are obtained by finding the eulerian circuits: -Ereg is the multiset of links between two oriented regions in SOR (including between the last and the first regions). Note that ∀ e ∈ Ereg, e ∈ Ereg denotes its reverse where Φ(e) = (u, v), Φ(e) = v, u (|Ereg| = 2|SOR|).

0 f → 1 f → 2 f → 1 r and 0 f → 1 f → 2 r → 1
-Φ : Ereg → {(u, v) | u, v ∈ V reg} 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), ∃! e ∈ Ereg | Φ(e) = (vreg -1 (r i ), vreg -1 (r j )).

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

N P-completeness

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

Definition 9.1: IRP decision problem (DIRP)

Given a set of contigs C, their multiplicities, a set of links L, a starting contig s, an integer k ∈ N, is there a valid sequence of oriented regions for IRP with ir∈IR (2replen(ir) -1) ≥ k, where IR is the solution set of inverted repeats? Proposition 9.1: DIRP is in N P Given the input of DIRP, a sequence of oriented regions SOR, the sequence of oriented contigs for each region COR, an integer k ∈ N. There is a polynomial time algorithm that checks if the given solution is valid, and if its accumulative lenght of the longest repeats equals at least k. Transform SOR in the sequence of oriented contigs SOC using COR. Verify the links and count the occurrences 7:

Proof

(prev_c, prev_or) ← SOC 8:

Count the contigs and verify the links

Proof

By reduction from the longest path decision problem from vertex s to vertex t (LPST P), known to be N P-complete [START_REF] Schrijver | Combinatorial Optimization: Polyhedra and Efficiency[END_REF]. Consider an instance I ∈ LPST P composed of a directed graph G = (V, E), two vertices s, t ∈ V and an integer k ∈ N (the hypothetical number of vertices between s and t in the longest path). We shall build an instance transform function mdgf such that I ∈ LPST P ⇐⇒ mdgf (I) ∈ DIRP. All four subgraphs G or,i in Figure 39b can be seen as copies of G V \{s,t} = V \ {s, t}, E V \{s,t} (the subgraph induced by the vertex set V \ {s, t}) in Figure 39a. For each or ∈ {f, r}, for each i ∈ {0, 1}, G or,i = V or,i , E or,i such that:

Function mdgf transforms graph G to graph G = (V , E ), vertex s ∈ V to vertices s f , s f ∈ V , t ∈ V to t f , t r ∈ V , and k to k = 2k -1.
-there is a bijective function vtrans or,i : V \ s, t → → V or,i , where ∀ v ∈ V or,i , vor(v) = or and vocc(v) = i;

-there is a bijective function etrans or,i : E V \{s,t} → → E or,i where ∀ (u, v) ∈ E V \{s,t} , ∃! (u , v ) ∈ E or,i such that:

vtrans or,i (u), vtrans or,i (v) = u , v if or = f v , u if or = r
There is a bijective function vsttrans st : {s, t} → → {(s f , s f ), (t f , t r )}. There is a function esttrans : {(s, w) ∈ E} ∪ {(u, t) ∈ E} → E such that:

122 -∀ (s, w) ∈ E: esttrans(s, w) ={(s f , vtrans f,i (w)) ∈ V 2 | i ∈ {0, 1}} ∪ {(vtrans r,i (w), s f ) ∈ V 2 | i ∈ {0, 1}} -∀ (u, t) ∈ E: esttrans(u, t) ={(vtrans f,i (u), t f ), (t f , vtrans r,i (u)) ∈ V 2 | i ∈ {0, 1}} ∪ {(vtrans f,i (u), t r ), (t r , vtrans r,i (u)) ∈ V 2 | i ∈ {0, 1}}
It is straightforward to see that there exists an algorithm in O(|V | + |E|) that computes this transform function.

The sets InvF , P InvF and AInvF are built based on G or,i graphs. Inverted fragments are visualised by green dashed vertical lines in Figure 39b, where

InvF = {(i, j) ∈ V f,0 × V r,1 | vtrans f,0 (i) = vtrans r,1 (j)}.
As the adjacent inverted fragments associate only vertices in V f,0 with those in V r,1 , the path that maximises the number of contiguous inverted fragments exits s f , goes through G f,0 to t f (or t r , it does not matter), and passes through G r,1 to s f . Since G f,0 is a copy of G V \{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. The lenght of this IR ( 12 (k + 1) = k) equals the hypothetical lenght of the longest path in LPST P.

To conclude, as there is a linear time complexity transform function mdgf such that I ∈ LPST P ⇐⇒ mdgf (I) ∈ DIRP, DIRP is at least N P-hard.

From Propositions 9.1 and 9.2 we conclude that DIRP is N P-complete.

Numerical results

We develop khloraascaf2 , a python package that computes the scaffolding of chloroplast contigs. It can either use Gurobi solver or CBC. All the following runs have been executed on a Linux laptop computer (32GB RAM, Intel ® Core™ i7-10610U CPU @ 1.80GHz ×8). Each time, Gurobi solver is selected.

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 N P-complete. Furthermore, the multimeric forms of the chloroplast genome is very often caused by the presence of an inverted repeat, that reverses the region(s) between it.

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

Perfect artificial data

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

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

Noisy artificial data

Here we test the behaviour of the solver when we add noise to the perfect artificial data. In that case, for each generated contig, the multiplicity has 25 % chance to be overestimated by one (that increases InvF , P InvF 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 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 N P-complete demonstration for a more general class of graphs.

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.

Input data generation

Here we briefly describe our protocol for input data generation 5 . 200 chloroplast genomes (selected in CpGDB 6 ) were downloaded from the NCBI 7 . For each of them, Table 14 -Gurobi solver metrics on perfect artificial growing data. |V|, |E|, |SC|, |IR| and |L| respectively stand for the number of vertices, edges, contigs in each single-copy region, contigs in each region of the inverted repeat, and links; Time: the presolve time plus the relaxation time (above) and the B&B time (below); Opt.: the linear relaxation bound U B (above) and the integer optimal value Opt (below); % Gap: the MIP gap equals 100 × (U B-Opt) U B

; Nodes: number of explored B&B nodes; Iter.: number of iterations for the LP relaxation (above) and for the B&B phase (below). a set of reads was generated. The contigs were generated with Minia (Chikhi and Rizk, 2012), a de Bruijn graph assembly approach. The links correspond to k-mer paths in the resulting compacted de Bruijn graph (cDBG) that connect two contigs. Finally, 31 instances were selected for which various difficulties have been detected, e.g. extra-links in the link set or combination of repeats.

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

To obtain the multiplicity of a given contig c, we sum the length of the alignments of the reads mapping on c (we denote by M apR 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: Where align read c is the sequence representing the alignment of the read read on the contig c. Its length equals the number of nucleotides of read that match on c (identity or substitution). Then the multiplicity mult(c) of c is obtained by normalising its coverage cov c by this of the starter s, cov s : mult(c) = max(1, cov c /cov s -0.1 ). As the multiplicity is an upper-bound of the usage of contig, we prefer to round up the normalisation only if the decimal part is greater than 0.1.

cov c = read∈M apRc align read c 126
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.

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 [START_REF] Gurevich | QUAST: Quality assessment tool for genome assemblies[END_REF] to evaluate the sequences associated to each genome form. As the genome reference is known, Quast tries to find the minimum number of differences (relocation, inversion, indels) between the reference and the sequence we provide. Three metrics are chosen to evaluate the best genome form for each problem succession:

-the genome fraction of the reference; -the number of misassemblies; -the number of local misassemblies.

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

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

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

Initial version

Table 16 provides all the metrics defined above for the 31 instances. The instances are solved very quickly (solver times < 4.5 sec., Table 19). khloraascaf successfully founds the sequence of the oriented contigs in 20 of them and retrieves the sequence of the oriented regions in 28 of them8 . 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: a maximum of the longest repeats The sequence of oriented regions for the Cucumix_hystrix's reference contains the following sub-sequences:

. . . , IR, IR , . . . , IR , SC, IR . As IRP is defined as finding a maximum number of the longest repeats, the "find the longest" part merges IR -IR and consequently not retrieves SC, normally inserted between IR and IR.

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

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

Nucleotide sequences misassemblies To avoid confusion, we always use nucl. seq. to denote "nucleotide sequence" to contrast with sequences of contigs / regions. In the following we analyse the instances presenting an unexpected behaviour for the Quast metrics, regarding if the sequences are found. Supplementary Table 18 permits verifying if the contigs that participate in the sequence have been correctly assembled. Except for Lathyrus_pubescens where one local misassembly is due to the missing contig in the sequence, all the (local) missassemblies in Commiphora_foliacea, Eucommia_ulmoides, Juniperus_scopulorum, Musa_ornata, Sciadopitys_verticillata, Taxus_baccata, Welwitschia_mirabilis provided by Quast are found in the nucl. seq. of links. 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).

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 bellow9 :

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.

Table 17 provides all the metrics obtained by running khloraascaf. The instances are also solved very quickly (solver times < 3 sec., Table 20). It truthfully finds all the sequences for the modified synthetic data except for Agathis_dammara instance. Although all the repeats (both the direct and the inverted ones) have been 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, La-thyrus_pubescens, Pelargonium_nanum and Triosteum_pinnatifidum just concern the nucl. seq. of links that is not a khloraascaf issue10 .

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

As our dedicated chloroplast scaffolding definition is a region-scaffolding-driven, the region graph is a natural data structure to reveal the multimeric genome forms that can exist 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 N P-complete in the general case, even though numerical results on perfect artificial data suggest there is a class of M DCG 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, khloraascaf11 , 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 N P-complete, the small size of the input data enables to quickly solve optimally CHSP.

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

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

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

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

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

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

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

Fragment graph

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

The literature describes three graph structures: a digraph, a bigraph and an ungraph. While the first one highlights visually the double-strand sequencing, its weakness is that it requires two vertices for each read [START_REF] Kececioglu | Exact and Approximation Algorithms for DNA Sequence Reconstruction[END_REF]. [START_REF] Myers | Toward simplifying and accurately formulating fragment assembly[END_REF] is the first to employ a bigraph to store overlaps. The key idea is to aggregate the two orientations of a read into only one vertex, a link and its reverse into only one edge. Finally, the ungraph structure associates one vertex for each of the two extremities of a fragment to simplify the graph traversal [START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF]. The ungraph handles two sets of edges. A fragment-edge connects the two vertices corresponding to the two extremities of a fragment (the tail and the head). Two vertices representing the extremities of a fragment are connected by a fragment-edge. The second one is the multiset of link-edges. A link-edge represents both a link and its reverse.

Thesis contribution

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

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

Short-term future work

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

Long-term future work

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

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

The GFA2 and PAF3 file formats, which store assembly graphs and matches between pairs of oriented fragments, could serve as a source of inspiration for disk storage methods. They are reminiscent of bidirected or undirected graphs. One strategy might be to keep the whole graph in this form on disk and convert a part of interest into a DGS or DGF in RAM.

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

Although this work illustrates the use of fragment graphs for the assembly process, other areas involving nucleotide sequences may benefit from this work. For instance, genome graphs only represent variation between genomes of individuals of the same species, but the presence of complementary reverse regions implies double-stranded behaviour. In this context, one could examine the adaptation of fragment graphs.

Chloroplast genome scaffolding

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

Thesis contribution

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

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

Moreover, we prove the Chloroplast Scaffolding Problem (CHSP) to be N Pcomplete in the general case, even though numerical results on perfect artificial data suggest there is a class of M DCG 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, khloraascaf4 , 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 N P-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.

Short-term future work

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

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

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

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

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

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

Long-term future work

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

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

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

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

diradj : E → E (u, v) → (u , v ) ctg u = ctg u or u = or u occ u = occ u + 1 if occ u | 2 occ u -1 otherwise ctg v = ctg v or v = or v occ v = occ v + 1 if occ v | 2 occ v -1 otherwise invadj : E → E (u, v) → (v , u ) ctg v = ctg v or v = 1 -or v occ v = 2 occ v 2 + (1 -or v ) ctg u = ctg u or u = 1 -or u occ u = 2 occ u 2 + (1 -or u )
2 Reduction of the repeated fragment sets

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

Repeated fragment set reductions

Here we consider the reduction operations for DirF and InvF , and conclude on their minimality. Only two vertices i, j ∈ 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.

A2

Reminder of pairs of direct fragments set definition:

P DirF =          (i, j), (k, l) ∈ DirF 2 s.t. ctg j < ctg k ∨ ctg j = ctg k ∧ occ j < occ k         
Proposition 2.7: Commutative reduction for P DirF ∀ ((i, j), (k, l)) ∈ P DirF , ((k, l), (i, j)) / ∈ P DirF Proof Let ((i, j), (k, l)) ∈ P DirF :

-if ctg j < ctg k thus ctg l > ctg i so ((k, l), (i, j)) / ∈ P DirF -else ctg j = ctg k ∧ occ j < occ k thus occ l > occ i so ((k, l), (i, j)) / ∈ P DirF Proposition 2.8

Pair of direct fragment set P DirF is a minimal set for DRP.

Proof

By reductio ad absurdum: if P DirF can be minimised to a P DirF set, then there exists a pair of direct fragments (p, q) ∈ P DirF | (p, q) / ∈ P DirF . Let p = (i, j) and q = (k, l), and let q 2 = (m, n) | (p, q 2 ) ∈ P DirF ∧ (q, q 2 ) ∈ P DirF . i, j, k, l, m, n vertices can simultaneously participate in the solution.

Let build a counter-example: suppose that we have the order i k l j m n. This implies the direct fragments p and q fall into a forbidden case for DRP, but this is not detected by Constraints C14 and C15 because (p, q) / ∈ P DirF . 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.

A6

However, for each order relation between the occurrences (occ u ≤ occ v and occ u ≥ 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 ≤ occ v . The reasoning is analogous for the second case.

Let a, b and c be three contigs in C such that all their identifiers are different, a mult = 2, b mult = mult(c) = 1, and (a f , b f ) ∈ L, (b f , a f ) ∈ L and (a f , c f ) ∈ L.

Let f a 1 = (i, j), f a 2 = (i , j ) be the two direct fragments associated to a and f b = (k, l), f c = (m, n) the two respectively associated to b and c. The idea behind the proof is to know if the adjacent direct fragments f a 1 f b f a 2 f c can participate into a solution path P ath. In such a case, P ath is sub-defined by edges i → k → i → m (and j → l → j → n). Remark that k can permute with l, and m with n.

According to ADirF set definition, (i, k) ∈ ADirF , (k, i ) ∈ ADirF and (i , m) ∈ ADirF . If we constrain the occurrences by the relation occ u ≤ occ v , it implies that:

occ i ≤ occ k (= occ l -1) occ l -1 ≤ occ i -1 occ i ≤ occ m
As f a 1 = f a 2 , occ i = 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) ≤ occ i (= 2) ≤ occ n (= 1), so it is absurd; else occ i > occ i -1 then when we focus on the first edges path, you can permute i and i and you fall into the same absurd implication as above.

Proposition 2.13: Occurences reduction for equal identifiers in ADirF

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

A10

∀ (i, j) ∈ InvF , ∀ (k, l) ∈ InvF | ctg i = ctg k ∧ occ i < occ k :

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

Proof

The reciprocals are trivial as AInvF ⊂ E. Let (i, j) ∈ InvF , (k, l) ∈ 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) ∈ E. The others follow the same logics.

-(i, k) ∈ AInvF as ∃ n ∈ 0, i mult Quast metric descriptions can be found at https://quast.sourceforge.net/do cs/manual.html#sec3.1. Here we adapt the description of the ones we use in this paper.

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

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

-flanking sequences overlap on more than 1 kbp;

-flanking sequences align to different strands.

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

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

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

Genome fraction (%) is the percentage of aligned bases in the reference genome. A base in the reference genome is aligned if there is at least one contig with at least one alignment to this base. Contigs from repetitive regions may map to multiple places, and thus may be counted multiple times 4 Supplementary results [50][51][52][92][93][94]99,100,104,119,A19 f Forward orientation 27,30,31,33,37,39,61,66,67,70,75,77,79,81,82, -30, 32, 34, 36, 38, 60, 61, 63, 65, 67-69, 72, 84, 86 L f (Multi)set of forward links 61, see also L & L r L Underlying set of the multiset L 28-30, 33, 37, 39, see also L & F L r (Multi)set of reverse links 61, see also L & L f L (Multi)set of links 28-30, 32-36, 38, 45, 46, 50-52, 60, 61, 63, 65, 67-70, 72, 84, 92-94, 99, 103, 119, 121, 126, 127, A20, G3, see also L & F {f, r} Set of orientations 27,28,33,46,66,68,[75][76][77]93,99,100,112,117,122,A3 r Reverse orientation 27,30,31,33,37,39,66,67,70,79,82, see also reverse & {f, r} R aw Set of raw reads 15

v1 scaffolding benchmark

• Reverse operation 27-29, 31, 34, 37, 46, 47, 57, 61, 63, 65, 67, 70, 71, 74, 75, 79-81, 83, 93, 94, 103-106, 108, 110-112, 117, 118, 124, 129 38,40,41,46,[48][49][50][51][52][53]70, see also undirected graph V Vertex set 30-34, 36, 37, 48, 61, 62, 64, 66-68, 70, 71, 74, 78, 80-83, G4, 21, 22, 39, 40, 48-50, 52, 56, 87, 92, 132, 135, 136, 138 , 22, 30-33, 35-37, 48, 50, 51, 60-69, 71, 72, 74, 78-84, 88, 136, 137, see also graph & vertex F forward Original nucleotide sequence orientation 14, 27, 35, 36, 61, 64-67, 75, 77, 78, 92, see also reverse fragment Generic nucleotide sequence 2, 9-14, 22, 26-30, 32, 33, 35-37, 60, 61, 63-65, 67-72, 75-79, 83, 84, 86-88, 135, 136 fragment assembly Generic term for the DNA fragment assembly, can also denote the whole process of assembling the reads to obtain the whole DNA sequences 9, 13,16,22,39,53,135 , 30, 32, 35, 36, 40, 60-72, 75, 76, 86-88, 93, 136, see also fragment link-edge Component of a graph that connects two vertices 35-37, 48, 50, 51, 70, 72, 88, 136, -15, 17, 18, 20, 21, 26, 39, 40, 53-55, 60, 86-88, 135, 136, see also sequencer read assembly Fragment assembly stage to obtain longer nucleotide sequence (contigs) from the reads 22,30,35,40,49,50,135, see also fragment assembly & contig G11 GLOSSARY reverse Reverse-complement nucleotide sequence orientation 2, 3, 14, 27-33, 35, 36, 61-66, 68, 69, 71, 75, 77, 80-83, 88, 92, 136, , 22, 30-33, 35, 36, 38, 48, 51, 53, 60-72, 74-84, 88, 136, 137 Dans ce manuscrit, nous comparons des structures de graphes représentant des relations de successions entre des séquences ADN orientées, utiles à différentes phases de l'assemblage. Puis, nous nous penchons sur le problème de l'échaffaudage dédié aux génomes de chloroplastes en proposant une nouvelle formulation, une résolution exacte et une implémentation.

Title: DNA fragment assembly: graph structures and chloroplast genome scaffolding Keywords: Genome assembly, Integer linear programming, Genome repeats, Multimeric forms Abstract: To obtain the nucleotide sequence of a DNA molecule, the molecule is fragmented using sequencing technology and the fragments are assembled. These fragments are called reads. They are subject to sequencing errors and must be considered in two orientations: that of their original DNA strand, or the reverse-complementary for the other strand. Assembly is based on pairwise overlaps between oriented reads and consists of three phases: assembling the reads to obtain contigs (sequences longer than the reads), scaffolding the contigs to obtain scaffolds (or-ders of oriented contigs), and completing the scaffolds (finding the nucleotide sequences separating the oriented contigs in the scaffolds).

In this manuscript, we compare graph structures representing succession relations between oriented DNA sequences, useful at different phases of assembly. Then, we address the scaffolding problem dedicated to chloroplast genomes by proposing a new formulation, an exact resolution and an implementation.

Figure 2 -

 2 Figure 2 -Different DNA conformations Each figure represents chromosomes, which are compressed DNA (double-stranded) molecules. ax = b means that b chromosomes are grouped in a-uplet: (a) haploid linear DNA 1x = 1; (b) haploid circular DNA 1x = 1; (c) diploid linear DNA 2x = 2.

Figure 3 -

 3 Figure 3 -Simplified DNA distribution in a plant's cell.The grey octagon container represents a plant's cell. The purple, blue and green oval shapes are respectively the plant cell's nucleus, mitochondria and chloroplasts. There are several mitochondria and chloroplasts in a same cell. Each of them possesses its own genomic material: in the nucleus, for this artificial example, the plant is diploid (2x = 4); each mitochondrion has several copies of its genome; each chloroplast has multiple forms of its own genome, linear, circular and multigenomes.

Figure 4 -

 4 Figure 4 -Inverted repeats causing multimeric forms: flip-flop inversion during the DNA replication.During the DNA replication of the chloroplast genome, one of the region between the inverted repeats can be reversed (here the red region SSC). This provokes the existence of several forms of the genome in the same chloroplast characterised as multimeric forms.

Figure 5 -

 5 Figure 5 -DNA repeats.A blue subsequence is the reverse-complement of a green subsequence, i.e. GCAA is the reverse of T T GC. (a) GCAA is direct repeated in the top strand, thus T T GC is direct repeated in the bottom strand. (b) GCAA is inverted repeated in the top strand, thus T T GC is inverted repeated in the bottom strand.

Figure 6 -

 6 Figure 6 -Common chloroplast genome's architectures. The most studied chloroplast genome form is circular and very often quadripartite. For each of the figures (a), (b) and (c), coloured arrows represent nucleotide sequences. LSC and SSC stand for Long Single-Copy and Short Single-Copy (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 Repeat 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 and illustrates that 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.

Figure 7 -

 7 Figure 7 -The shotgun sequencing approach.From the bottom to the top of the figure: the reads (in blue) are sequenced from longer fragment (in green) randomly sampled from copies of the genome (in orange). Based on the overlapping reads, the assembly process aims to retrieve the genome. The coverage counts the average number of times a base in the genome is covered by the reads.

Figure 8 -

 8 Figure 8 -Sequencing error types. Three types of errors can occur during the sequencing of DNA fragments. For each subfigure, the top sequence is the original/template sequence and the bottom sequence is the sequenced one. (a) Substitution: one nucleotide is changed to another one; (b) insertion: there is one extra nucleotide; (c) deletion: there is one missing nucleotide.The first error type is someone noted as a mismatch, and the two latters are categorised as indels.

Figure 9 -

 9 Figure 9 -Paired-end sequencing. The two main grey lines are the two strands of a fragment. The two reads r 1 and r 2 are sequenced by pairs distanced by d nucleotides (inner distance).

Figure 10 -

 10 Figure 10 -Sequence alignments. Sequence alignments distributed into three categories. u and v are two sequences. The arrows give their orientation (tail: 3', head: 5'). The blue lines between them correspond to aligned nucleotides. (a) Overlap: the head of u is aligned with the tail of v. (b) Containment: u is containing v. (c) Share: a subpart of u is aligned with a subpart of v, and at least one of the subpart does not begin from or finish to an extremity (here each subpart respects the statement).

Figure 12 -

 12 Figure 12 -Impact of the genome repeats on the fragment sequencing. The two subfigures show the same genome sequenced respectively with long and short reads. The genome possesses a repeated region R (in red) and three single-copies (resp. yellow, blue and green). (a) Reads b and d are covering the repeats entirely: the assembly gives abcde. (b) Reads c and f are subsequence of the repeats: the assembly can give abcdef gh, or abf decgh (if the repeats are not exactly the same, switching c and f can misassemble the genome), or even abcgh that shortens the genome.

Figure 13 -

 13 Figure 13 -OLC versus DBG fragment assembly approaches. (a) The genome possesses a repeated region R. The sequencing produces nine reads.For the example, reads 3 and 7 are exactly the same and are contained in R. (b) From the overlap of the read (at least two bases), the string graph is built. Note that neither read 3 nor read 7 are removed because no one is longer than the other while they are mutually contained each other. (c) The reads are split into k -1 overlapping words of length k (k-mers). In the (partial) edge-centric dbgraph, vertices are the k -1 exact overlaps and edges are the k-mers. Here k = 3. For (b) and (c): unique paths are collapsed to form unitigs. The over-simplified illustrations for the OLC and respectively the DBG approaches suggest that the genome can be found by solving a Hamiltonian path through the simplified string graph resp. an Eulerian path through the simplified dbgraph.

  and a have complementary orientations in {f, r}); seq(a) = seq(a) (the nucleotide sequence of a is the reverse of this of a).

Figure 15

 15 Figure15illustrates the fragment sets and their associated functions.

Definition 1 .

 1 5: Link setDenote by L = (L , m L ) the multiset of ordered pairs of fragments, where L ⊂ F 2 is a finite set, and m L : L → N >0 gives the multiplicity of each element of L in L. Also let linkp : L → → L be the function that gives the fragment couple associated with a link. Denote by Σ L = 0, |L| the set of link integer labels and its associated function lid : L → → Σ L .

Figure 14 -

 14 Figure 14 -Examples of links Three examples of links. u and v are two fragments. The arrows give their orientation (tail: 5', head: 3'). (a) The head of u overlaps the tail of v. (b) Paired-end: u is linked with v, separated by a distance approximately of d nucleotides. (c) Hi-C: u is linked with v, separated by an unknown distance.

Figure 15 -

 15 Figure 15 -Fragment and link sets with functions Top left, in green: unoriented fragment sets. Middle, in blue: oriented fragment sets. Bottom, in yellow: Link (multi)sets.

Figure 15

 15 Figure15illustrates the link sets and their associated functions. In Sections 1.2 to 1.4 we describe the three graph structures for the fragments and their links. They are loopless multigraphs.

Figure 16

 16 Figure16illustrates the link case representations in DG.

(

Figure 16 -

 16 Figure 16 -Link cases in DG. Fragments a f , b f and their reverse a r , b r can link in four manners. Their associated vertices are u f , v f , u r and v r , respectively. Each subfigure illustrates one link in DG.

  .6 for BG. Definition 1.11: The reverse operation in BG Let e ∈ E be an edge, Φ(e) = {u, v}, a = unf rag(u) and b = unf rag(v). Assume that unf id(a) < unf id(b). As function link is bijective, it holds that: link(v, u, e) = link(u, v, e); attre(e) = (f or(a), f or(b)) and attre(e) 0 • • 0 = f or(b), f or(a) .

  a

a

  is the exclusive disjunction notation.

Figure 19 -

 19 Figure 19 -A path in BG.A valid path is coloured in red. It can be read in two possible ways: from vertex v or from vertex x. Let p v = vux be the one that begins from v and p x = xuv be the one that begins from x. Note that p v is valid if and only if p x is valid. In (a) edges are bidirected while in (b) they are undirected.

  and ARKS are pre-scaffolder methods. Suggested in supplementary method.

Figure 23 -

 23 Figure23-Scaffolding input data. For each subfigure, c i and c j are two contigs and the vertical grey lines represent nucleotide alignment. (a) The two paired-reads r i and r j map on two different contigs. It implies that c i forward is followed by c j forward, separated by a nucleotide distance d . (b) This case is the same as (a) except that for Hi-C data, the distance d is unknown, implying an unknown distance d . (c) c i and c j map on two overlapping long reads lr i and lr j . It implies that c i forward is followed by c j forward, separated by an approximated nucleotide distance d . (d) G is a reference genome (e.g. of a species near to the organism of interest). Mapping c i and c j on G provides an orientation of the contigs and an order with a distance d . (e) OP M is an optical map with three markers (yellow square, green diamond and purple circle). Finding the same markers orders in c i and c j enables to estimate their orientation and their order with a distance d . (f) r i , r j , and r k , r l are two paired-end reads with the same barcode (yellow circle), so they are sequenced from the same fragment. Because r i , r j map on the head of c i , and r k , r l map on the tail of c j , it is possible to orient and order c i and c j .

  As proposed for read assembly OLC approach,[START_REF] Huson | The greedy path-merging algorithm for contig scaffolding[END_REF],[START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF] and[START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF] compute a transitive reduction of the links based on their distance distribution.Aggressive removing[START_REF] Roy | SLIQ: Simple Linear Inequalities for Efficient Contig Scaffolding[END_REF] predict the orientations based on the links for which the two contigs have the same orientation. A common aggressive link removing is to consider only one link between two contigs of the four cases.[START_REF] Dayarian | SOPRA: Scaffolding algorithm for paired reads via statistical optimization[END_REF],[START_REF] Luo | BOSS: A novel scaffolding algorithm based on an optimized scaffold graph[END_REF] and[START_REF] Qin | LRScaf: Improving draft genomes using long noisy reads[END_REF] keep the maximum-weighted link, in addition,[START_REF] Koren | Bambus 2: Scaffolding metagenomes[END_REF] remove the smallest number of unconfident alignments in the bundle according the distance distribution.[START_REF] Coombe | ARKS: Chromosome-scale scaffolding of human genome drafts with linked read kmers[END_REF] keep only one link if it passes a binomial test of being the real link among the others.Ghurye et al. (

Definition 2. 1 :

 1 Orientation-valid layout Denote by C ⊂ = {c 0 , . . . , c n-1 } a subset of the contig set C and by L ⊂ a subset of the link set L. Without considering an ordering of the contigs in C ⊂ , L ⊂ is an orientation-valid layout if and only if there are exactly two orientation assignments o = (o 0 , . . . , o n-1 ) ∈ {f, r} |C ⊂ | and o = o 0 , . . . , o n-1 and there exists a balanced partition of L

Figure 24 illustrates

 24 Figure24illustrates Definition 2.1 for the bigraph (BG) and for the ungraph (U G) as defined in Sections 1.3 and 1.4, and illustrated in Figure22.

  Figure 24 -Bidirected and undirected graph cycles. (a) and (c): there are exactly two vertex orientation assignments (o = {u f , v r , w f }, and the reverse orientations o) such all the valid paths assign vertex orientations respecting o, otherwise o. (b) and (d): there are at least two valid paths that imply two different orientation assignments with a non-null intersection, e.g. o = {u f , v r , w f } and o = {v r , w f , u r }. It is sufficient to remove one of the red or the yellow link to get a valid-orientation layout.

  For the positioning of the contigs,[START_REF] Roy | SLIQ: Simple Linear Inequalities for Efficient Contig Scaffolding[END_REF] employ a majority voting on the ordering predictions based on inequalities.Randomized Greedy Algorithms for the maximum matching problem Finally,[START_REF] Ghurye | Integrating Hi-C links with assembly graphs for chromosome-scale assembly[END_REF] benifiting from the algorithm proposed in[START_REF] Poloczek | Randomized Greedy Algorithms for the Maximum Matching Problem with New Analysis[END_REF] running in O(|L|) with U G. The total complexity of their approach raises to O(|C||L|).

Proposition 1 .

 1 1 gives the amount of octet DGA consumes. Proposition 1.1: DGA memory consumption The memory size M em(DGA) of DGA (in octets) is equal to: M em(DGA) = 2P(|F| + 1) + 2|L| oct(|F|) + oct 1 2 |L|

2

  Algorithms for DGS, DGF and BGU

1

 1 For the sake of clarity, rev(ind) = ind.Algorithm 2: Find place of edge index in the neighbour list Require: List of tuple list, edge index e, index t of the place of the edge indices in the tuples contained in list. The edge index must be in a tuple. Ensure: Returns the index ind of the tuple containing e in list.1: function get_index(list, e, t) List list, an entry index ind ∈ 0, |list| . Ensure: Only the element in list[ind] is removed, and |list| decreases by one.1: function delete_list_entry(list, ind) ] ← list.pop( )

Algorithm 4 :

 4 Iterate over the predecessors for DGS Require: A vertex label v ∈ Σ V . Ensure: Returns a list containing tuples in Σ V × Σ E representing the predecessors of the vertex labelled by v. 1: function dgs_preds(v) vertex label v ∈ Σ V associated with an orientation vor ∈ {f, r}. Ensure: Returns a list containing tuples in Σ V × {f, r} × Σ E representing the predecessors of the vertex labelled by v with fragment orientation vor. 1: function dgf_preds(v, vor) vertex label v ∈ Σ V associated with an orientation vor ∈ {f, r} representing the predecessors of the vertex labelled by v and with fragment orientation vor. Ensure: Returns a list containing tuples in Σ V × {f, r} × Σ E representing the predecessors of the vertex labelled by v with the fragment orientation vor. 1: function bgu_preds(v, vor) 2: preds ← empty_list( ) 3: for (u, e) ∈ N v do 4: if u < v ∧ Eattr[e][1] = vor then 5: preds.append(u, Eattr[e][0], e)

Algorithm 7 : 1 :

 71 Iterate over the successors for DGS Require: A vertex label v ∈ Σ V . Ensure: Returns a list containing tuples in Σ V × Σ E representing the successors of the vertex labelled by v. 1: function dgs_succs(v) vertex label v ∈ Σ V associated with an orientation vor ∈ {f, r}. Ensure: Returns a list containing tuples in Σ V × {f, r} × Σ E representing the predecessors of the vertex labelled by v with fragment orientation vor. 1: function dgf_succs(v, vor) vertex label v ∈ Σ V associated with an orientation vor ∈ {f, r} representing the successors of the vertex labelled by v and with fragment orientation vor. Ensure: Returns a list containing tuples in Σ V × {f, r} × Σ E representing the predecessors of the vertex labelled by v with the fragment orientation vor. function bgu_succs(v, vor) e) ∈ N v do 4: if v < w ∧ Eattr[e][0] = vor then 5: succs.append(w, Eattr[e][1], e) 6: else if 1 -Eattr[e][1] = vor then 7: succs.append(w, 1 -Eattr[e][0], e)

Algorithm 10 :

 10 Add a new vertex in the graph for DGS Ensure: The length of N + V is increased by two and returns the last vertex label in fragment orientation forward. 1: function dgs_add_vertex( ) The lengths of N - V and N + V are increased by one and returns the last vertex label. 1: function dgf_add_vertex( ) The length of N V is increased by one and returns the last vertex label. 1: function bgu_add_vertex( ) 2:

Algorithm 13 :

 13 Add a new edge and its reverse in the graph for DGS Require: Two vertex labels u and v. Ensure: Returns the new edge index e such that (u, e) ∈ preds(v) and (v, e) ∈ succs(u). 1: function dgs_add_edge(u, v) 2:

  Two vertex labels u and v, associated with their respective fragment orientation uor and vor. Ensure: Returns the new edge index e such that (u, uor, e) ∈ preds(v, vor) and (v, vor, e) ∈ succs(u, uor). 1: function dgf_add_edge(u, uor, v, vor) 2:

1 Algorithm 15 :

 115 Add a new edge and its reverse in the graph for BGU Require: Two vertex labels u and v, associated with their respective fragment orientation uor and vor. Ensure: Returns the new edge index e such that (u, uor, e) ∈ preds(v, vor) and (v, vor, e) ∈ succs(u, uor). 1: function bgu_add_edge(u, uor, v, vor)

  vertex label v. Ensure: The edges containing the vertex or its reverse are deleted from the neighbour lists. The length of N V is reduced by one. 1: function bgu_delete_vertex(v)

  Two vertex labels u and v associated with their fragment orientation uor and vor, respectively, and a edge label. Ensure: The edges are deleted for the two vertices.1: function dgf_delete_edge(u, uor, v, vor, e)

  Two vertex labels u and v and a edge label. Ensure: The edges are deleted for the two vertices.1: function bgu_delete_edge(u, v, e)2:

Figure 31 -

 31 Figure 31 -Repeat degeneration and region orders. 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 their degeneration on their structures. We give the bellow 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.

Figure 32 -

 32 Figure 32 -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.

Figure 33 -

 33 Figure 33 -M DCG example.Each vertex is associated with an occurrence of an oriented contig, and each contig is represented by an even number of vertices. For example, vertex labelled c r,1 means that it comes from contig c = contig(c r,1 ), in its reverse orientation (vor(c r,1 ) = r), and in its second occurrence (vocc(c r,1 ) = 1). The colours are the same as the ones in Figure6ato 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.

( a )Figure 34 -

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

Figure 35 -

 35 Figure 35 -Adjacent repeated fragment sets examples.The two subfigures represent the multiplied link (and its reverse) ((c, f ), (d, f )) ∈ 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∀ (u, v) ∈ E, v, u ∈ 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 ).

Figure 36 -

 36 Figure 36 -Non-exhaustive illustrations for authorised and forbidden order cases for two repeated fragments ((i, j), (k, l)) ∈ P RepF . (a) and (b) Authorised and forbidden orders for P DirF ; (c) and (d) Authorised and forbidden orders for P InvF .

Figure 37 -

 37 Figure 37 -Extracting the genome architecture in M DCG 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 25 returns four regions: m = 1 inverted repeat (the two green arrows pointing downwards) and n = 2 single-copies (the purple and the red arrows). COR = ((a f ), (c r , b f ), (d f )) and SOR = (0 f , 1 f , 2 f , 1 r ), where x y = (x, y) for clarity.

1 :

 1 the initial vertex of the single-copy region containing the starting vertex Require: The starting vertex s, the set of fixed variable values x * e , e ∈ 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. function initial_vertex( ) Two vertices u, v ∈ V , the set of fixed variable values x * e , e ∈ 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_contiguous(u, v, prev_region_code, region_code) Graph M DCG, the starting vertex s, the set of fixed variable values x * e , e ∈ E, the chosen direct fragments DirF * and the chosen inverted fragments InvF * . Ensure: Returns the oriented region sequence SOR, and the oriented contig sequence for each forward oriented region COR. 1: function ilp_solution_to_regions( ) 2: prev_region_code ← SC 3: SOR ← [0 f ] 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

  Figure38-Region graph for the toy example.The toy example's solution given in Figure37results in the graph visualised in Figure33will give a region graph that shows an LSC-IR-SSC-IR architecture as the one given in Figure6a. Two oriented region sequences (genome forms) are obtained by finding the eulerian circuits:0 f → 1 f → 2 f → 1 r and 0 f → 1 f → 2 r → 1 r .They correspond to two multimeric forms that are commonly found in the chloroplast cells, and described in Figure4. (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.

  = (s, f ) then The circuit must start with the starting contig in

  or) ∈ SOC do 10:if((prev_c, prev_or), (c, or)) / N P-hardness of IRP decision problem DIRP is N P-hard.

Figure 39 -

 39 Figure 39 -From a digraph G for LP ST P to a digraph G for IRP. Bold red edges in both sub-figures correspond to the solution path for LP ST P and IRP problems, respectively. (a) G V \{s,t} is the subgraph induced by the vertex set V \ {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 r1 visualise the inverted fragments.

  Figure 39 illustrates the transformation.

Figure 40 -

 40 Figure 40 -Solver running time distributions for perfect artificial data.Points are measured times, the red curves correspond to the best ae bx functions applied on IR length axis.

2

  such that occ i = 2n + or i = 2n, and∃ n ∈ 0, k mult 2 such that occ k = 2n + or k = 2n + 1; -(k, i) / ∈ AInvF as occ k > occ i ; -(l, j) / ∈ AInvF as occ l > occ j ;-(j, l) / ∈ AInvF as or j = 0 and occ j 2 so n ∈ 0,

FragmentC

  Set of contigs 45

alignment

  Nucleotide comparison between at least two sequences assembly graph One of the output of genome assembly method see also read assembly B bidirected graph Graph that contains 3 edge types 32, 47, 87, 88, 136 bidirected multigraph A graph where the same bioriented edge between the same two vertices can occur several times 32, 33, 88, 136, see also multigraph & bidirected graph C chloroplast Organelle in the plants'cells specialised in the photosynthesis process 2, 4-8, 22, 23, 53, 92, 135 chromosome A chromosome is a DNA molecule that contains the genetic information for an organism 4, 13, 21, 26, 135, see also genome & DNA contig Result of the assembly of the reads, nucleotide sequence longer than the reads

  Séquençage ADN et Assemblage de génomes . . . . . . . . . . . . iii Graphes de fragments ADN . . . . . . . . . . . . . . . . . . . . . . iv Échafaudage de génomes de chloroplastes. . . . . . . . . . . . . . v Conclusions et perspectives . . . . . . . . . . . . . . . . . . . . . . vii State-of-the-art. . . . . . . . . . . . . . . . . . . . . . . . . . 91 1.3 Our approach. . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Input data and notation . . . . . . . . . . . . . . . . . . . . . . 92 2.1 Set of contigs C . . . . . . . . . . . . . . . . . . . . . . . . . 92 2.2 Set of links L . . . . . . . . . . . . . . . . . . . . . . . . . . 93 2.3 Mathematically defining genomic regions . . . . . . . . . . . . . 93 Chloroplast scaffolding problem formulations . . . . . . . . . . 94 Graph and repeated fragment sets . . . . . . . . . . . . . . . . 98 4.1 Graph structure . . . . . . . . . . . . . . . . . . . . . . . . . Unknown fragment orientations . . . . . . . . . . . . . . . . . . 4.3 Sequence similarity: single-copy or repeat? . . . . . . . . . . . . 4.4 True sequence divergences or sequencing errors?. . . . . . . . . . The Shortest Common Superstring . . . . . . . . . . . . . . . . 5.2 Overlap-Layout-Consensus . . . . . . . . . . . . . . . . . . . . 5.3 De Bruijn Graph approach . . . . . . . . . . . . . . . . . . . . 5.4 Breaking down the fragment assembly problem . . . . . . . . . . Graph structure for read assembly and scaffolding stages . . . . . 22 6.2 Scaffolding of multimeric chloroplast genome forms . . . . . . . . 22

	Acknowledgement Résumé en français Contents List of Figures List of Tables I Introduction 17 Sequence and Quast metrics for the modified synthetic data i iii ix xv xvii 1 version. Appendix 18 Benchmark 3 v1 contig Quast 19 Benchmark 3 v1 ILP stats 20 Benchmark 3 v2 ILP stats 5 Fragment assembly approaches 5.1 6 Addressed research topics 6.1	132 A17 A17 A18 A20 22
	98
	xi

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 1.1 Chloroplast genome specificities . . . . . . . . . . . . . . . . . . 90 1.2 xviii I INTRODUCTION Joe Hisaishi. (1999). Inner Voyage [Song]. On The Universe Within, Vol.1 & 2. PONY CANYON In this chapter 1 DNA overview 1.1 A recipe for cellular mechanisms . . . . . . . . . . . . . . . . . 1.2 On the motivations for knowing the whole DNA sequence . . . . . 1.3 Various molecular conformations structure the genome . . . . . . 2 Chloroplast genome 2.1 Genome division . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Genome multimeric forms. . . . . . . . . . . . . . . . . . . . . 2.3 Repeats and single-copy regions. . . . . . . . . . . . . . . . . . 2.4 Genome evolution . . . . . . . . . . . . . . . . . . . . . . . . 3 DNA sequencing 3.1 The Shotgun sequencing approach . . . . . . . . . . . . . . . . 3.2 First generation: Sanger and BAC technologies . . . . . . . . . . 3.3 Second generation: high-throughput sequencing . . . . . . . . . . 3.4 Third generation: single molecule run time . . . . . . . . . . . . 3.5 Supplementary sequencing data . . . . . . . . . . . . . . . . . . 4 The challenges of fragment assembly 4.1 Read sequence alignment . . . . . . . . . . . . . . . . . . . . . 4.2

  Notations and fundamental definitions . . . . . . . . . . . . . .

		2.5.3	The maximum (weighted) matching and spanning tree . . . .	50
	2.6 Solving approaches . . . . . . . . . . . . . . . . . . . . . . . .	51
		2.6.1	Greedy approaches . . . . . . . . . . . . . . . . . . . . .	51
	2.6.2 Richard Galliano & Miran Vaupotic & Czech Phil-Fix parameters dynamic programming . . . . . . . . . . . . 2.6.3 Mathematical programming. . . . . . . . . . . . . . . . . harmonic Quartet. (2023). Petite Suite Française: 3 Chloroplast genome assembly III. Espiègle [Song]. On Madreperla. Navona Re-cords 3.1 Chloroplast sequence extraction . . . . . . . . . . . . . . . . . .	52 52 53 53
		3.1.1	Filtering the reads . . . . . . . . . . . . . . . . . . . . .	54
	3.1.2 In this chapter Filtering the contigs . . . . . . . . . . . . . . . . . . . . 3.2 Chloroplast reads assembly . . . . . . . . . . . . . . . . . . . .	55 55
		3.2.1	De Bruijn graph approach . . . . . . . . . . . . . . . . .	55
		3.2.2	Seed-and-extend . . . . . . . . . . . . . . . . . . . . . .	56
	1 Graph structure for fragment assembly 3.3 Chloroplast scaffolding . . . . . . . . . . . . . . . . . . . . . .	56
	Fragment set . . . . . . . . . . . . . . . . . . . . . . . 3.4 Chloroplast assembly validation . . . . . . . . . . . . . . . . . . 1.1 1.1.1	57
	1.1.2 1.2 Directed graph (DG): oriented fragments based . . . . . . . . . . Link set. . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Bidirected graph (BG): oriented walk based . . . . . . . . . . . . 1 Graph structure for fragment assembly
	1.4 Undirected graph (UG): tail-head fragments based. . . . . . . . .
	2 Scaffolding the contigs
	2.1 Scaffolding input data . . . . . . . . . . . . . . . . . . . . . .
	2.2 Subsampling the input data . . . . . . . . . . . . . . . . . . . .
	2.2.1	Bundling the links . . . . . . . . . . . . . . . . . . . . .
	2.2.2	Removing paired-end reads . . . . . . . . . . . . . . . . .
	2.2.3	Removing contigs . . . . . . . . . . . . . . . . . . . . .
	2.2.4	Removing links . . . . . . . . . . . . . . . . . . . . . .
	2.2.5	Partitioning the instances. . . . . . . . . . . . . . . . . .
	2.3 Orienting the contigs . . . . . . . . . . . . . . . . . . . . . . .
	2.3.1	Maximising the sum of used link weights. . . . . . . . . . .
	2.3.2	Remove the minimum number of contigs and links to avoid odd
		reversal cycles by minimising the weight sum . . . . . . . . .
	2.3.3	Maximum log-likelihood function . . . . . . . . . . . . . .
	2.4 Ordering the oriented contigs . . . . . . . . . . . . . . . . . . .
	2.4.1	Maximising the sum of used distances bunches' weights . . . .
	2.4.2	Minimum spanning tree . . . . . . . . . . . . . . . . . .
	2.4.3	Contig positioning only. . . . . . . . . . . . . . . . . . .
	2.5 Orienting and ordering the contigs simultaneously . . . . . . . . .
	2.5.1	Maximising the sum of distance bundles' weights to order linearly
		the contigs . . . . . . . . . . . . . . . . . . . . . . . .
	2.5.2	The heaviest matching, paths and cycles . . . . . . . . . . .

  [START_REF] Donmez | SCARPA: Scaffolding reads with practical algorithms[END_REF] remark that the contig can be misassembled and conclude that the orientation problem should consider removing the contig and the links. Consequently, they associate the orientation problem with finding a

minimum odd cycle traversal in a transformed U G. Indeed, to generalise with the removing of contigs, they transform each link-edge (u, v) into a path uu v v where u and v are two new artificial vertices. Although finding a minimum odd cycle traversal is a N P-hard problem,

[START_REF] Lokshtanov | Simpler Parameterized Algorithm for OCT[END_REF]

propose a fixed-parameter tractable algorithm in O(3 k k|E||V |), where k is the number of vertices to suppress, and V and E are respectively the vertex and the edge sets.

  David Krakauer & Kathleen Tagg. (2020). The Geyser [Song]. On Breath & Hammer.Table Pounding Records Iterating over the predecessors . . . . . . . . . . . . . . . . . . 2.3 Iterating over the successors . . . . . . . . . . . . . . . . . . . 2.4 Adding a vertex . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Adding an edge. . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Deleting a vertex . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Deleting an edge . . . . . . . . . . . . . . . . . . . . . . . . .

	5 Conclusions and perspectives	86
	Section 1 describes several graph implementations for each structure from Chapter II,
	Section 1 (see the fragment graph structures overview in Figure
	In this chapter
	1 Implementations
	3 Time costs
	4 Memory and time cost comparisons

1.1 Directed graph (DG): oriented fragments based . . . . . . . . . .

1.1.1 All oriented fragments directed graph (DGA) . . . . . . . . . 1.1.2 Oriented fragments' successors directed graph (DGS) . . . . . 1.1.3 Forward fragments directed graph (DGF) . . . . . . . . . . 1.2 Bidirected graph (BG): oriented walk based . . . . . . . . . . . . 1.2.1 Unoriented fragments bidirected graph (BGU) . . . . . . . . 1.2.2 Transformation to DGF . . . . . . . . . . . . . . . . . . 1.3 Undirected graph (UG): tail-head fragments based. . . . . . . . . 1.3.1 All oriented fragments undirected graph (UGA). . . . . . . . 1.3.2 Transformation to DGS . . . . . . . . . . . . . . . . . . 1.4 Fragment graph map . . . . . . . . . . . . . . . . . . . . . . . 2 Algorithms for DGS, DGF and BGU 2.1 Subfunctions . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2

  pop( )

	26:	Update the vertex index
	27:	for (u_rev, e) ∈ N + v_rev do
	28:	u ← u_rev
	29:	adj_ind ← get_index(N + u , e, 1)
	30:	N + u [adj_ind] ← (v, e)
	31:	for (w_rev, e) ∈ N + v do
	32:	w ← w_rev
	33:	adj_ind ← get_index(N + w , e, 1)
	34:	N + w [adj_ind] ← (v_rev, e)
	Algorithm 17: Delete a vertex for DGF
	2:	Delete it from its predecessors
	3:	for (u, uor, e) ∈ N -v do
	4:	if uor = f then
	5:	adj_ind ← get_index(N + u , e, 1)
	6:	delete_list_entry(N + u , adj_ind)
	7:	else
	8:	adj_ind ← get_index(N -u , e, 1)
	9:	delete_list_entry(N -u , adj_ind)
	10:	card_edges ← card_edges -2
	11:	delete(N -v )
	12:	Delete it from its successors
	13:	for (w, wor, e) ∈ N + v do
	14:	if wor = f then
	15:	adj_ind ← get_index(N -w , e, 1)
	16:	delete_list_entry(N -w , adj_ind)
	17:	else
	18:	adj_ind ← get_index(N + w , e, 1)
	19:	delete_list_entry(N + w , adj_ind)
	20:	card_edges ← card_edges -2
	21:	delete(N + v )
	22:

Require: A vertex label v. Ensure: The edges containing the vertex or its reverse are deleted from the neighbour lists. The lengths of N - V and N + V are reduced by one. 1: function dgf_delete_vertex(v)

Table 4 -

 4 Calculus detail of basic operations costs. c(x) is the cost of operation x. Add x to the end of the list list list.pop() 1 Delete the last element of the list list delete(list) |list| Delete all the list for i ∈ a, b do x (b -a + 1) × c(x) For loop algorithmic costs of Algorithms 10 to 15.Table 8 gives the algorithmic costs of Algorithms 16 to 18. Table 9 gives the algorithmic costs of Algorithms 19 to 21.

	Operation	Cost Description
	x	0 Memory access
	x ± y	1 + c(x) + c(y) Basic operation
	x ← y	1 + c(y) Affectation
	|x|	1 + c(x) Absolute value or element's length
	if x y	1 + c(x) + c(y) Conditional
	empty_list()	1 Empty list constructor
	list.append(x)	1 + c(x)

Table 5 -

 5 Algorithmic costs for subfunctions. |list| is the length of the input list.

		Best	Worst	Average
	rev		2
	get_index	4 1 + 5|list| 1 + 5	|list| 2
	delete_list_entry	3	4	4 -	1 |list|

Table 6 -Algorithmic costs of iterating over the neighbours.

 6

Table 8 -

 8 Algorithmic costs of deleting a vertex.

	Best	Worst	Average
	DGS 28k + 6 20k 2 + 36k + 8	10

Table 9 -Algorithmic costs of deleting an edge.

 9

		Best	Worst	Average
	DGS	21 10k + 17 5k + 17 -2 k
	DGF	22 10k + 18 5k + 18 -2 k
	BGU	18 20k + 14 10k + 14 -1

k 4 Memory and time cost comparisons

Table 10

 10

provides memory consumption and elementary operation time cost comparisons between DGS, DGF , BGU .

  Complexity validation on artificial data . . . . . . . . . . . . . . 124 10.1.1 Perfect artificial data. . . . . . . . . . . . . . . . . . . . 124 10.1.2 Noisy artificial data . . . . . . . . . . . . . . . . . . . . 124 10.2 Synthetic chloroplast input data. . . . . . . . . . . . . . . . . . 125 10.2.1 Input data generation . . . . . . . . . . . . . . . . . . . 125 10.2.2 The evaluation's metrics . . . . . . . . . . . . . . . . . . 127 10.2.3 Initial version . . . . . . . . . . . . . . . . . . . . . . . 128 10.2.4 Modified version . . . . . . . . . . . . . . . . . . . . . . 131

	6 Hierarchical problem succession	111
	7 From an ILP solution to a genome structure	112
	8 Multiple genome forms	116
	9 N P-completeness	119
	10 Numerical results	123
	10.1 11 Conclusion	132
	12 Discussion and perspectives	133

Table 11 -Toy example of input data. (a)

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

Table 12 -ILP sets and functions corresponding table.

 12

	ILP	RepF PRepF	ARepF	Forbid	Alpha	repfrag	repadj
	DRP DirF P DirF ADirF F orbidDR AlphaDR dirf rag diradj
	IRP InvF P InvF AInvF F orbidIR	AlphaIR	invf rag invadj

  and IRP models 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|).

	max	m p +	isadj e
	p∈P RepF	e∈ARepF
	s.t. CCircuit
	CRepeat
	(C26)
	(C27)	if no repeats to fix
	CF ixRegions otherwise
	Traditionally, SCP finds the maximum weighted circuit.
	SCP model
	max	wex v i v
	v∈V \{s,s}
	s.t. CCircuit
	(C26)
	CF ixRegions if repeats to fix

Verify the validity of the sequence of oriented regions Require:

  It first retrieves the total sequence of oriented contig SOC from the sequence of oriented regions SOR and the sequence of each forward regions COR. It traverses the oriented contigs in SOC, verify if each two consecutive oriented contigs are linked in the link set, and count the number of times a contig occur. At the end, it checks if each contig does not appear more than its multiplicity number of times.It is straightforward to see that Algorithm 26 and Algorithm 27 are linear according the size of SOR and COR. Remember we assume that the belonging test "is x ∈ X?" for an object x in a set X is in Θ[START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF]. 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.

	Algorithm 27:
	Algorithm 26 verifies if the sequence of oriented regions is valid. It requires two traversals of SOR: (i) to identify which regions in the sequence form IRs; (ii) to verify if the order between the IRs is valid (thanks to the use of a LIFO). It also checks if the associated accumulative length of the longest repeats equals at least k. A LIFO ir_lif o is associated with four methods: ir_lif o.put(x) append x to the LIFO; ir_lif o.is_empty returns true if the LIFO is empty; ir_lif o.peek returns the last value in the LIFO; ir_lif o.pop deletes the last value in the LIFO and returns it. Algorithm 27 verifies if the corresponding sequence of oriented contigs is 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 The repeats are pairs of regions 12: IR.add(rid) 13: count_reg.append(ror) 14: Verify parameter k 15: if ir∈IR (2|COR[ir]| -1) < k then 16: return False 17: Verify the order between IRs 18: ir_lif o ← lifo( ) 19: for (rid, ror) ∈ SOR do 20: if rid ∈ IR then 21: if ir_lif o.is_empty( ) ∨ rid = ir_lif o.peek( ) then 22: ir_lif o.put(rid) valid. Algorithm 26: 2: 23: else
	24:	ir_lif o.pop( )
	25:	return ir_lif o.is_empty( )

1: function is_sor_valid( )

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:

Table 15 -Gurobi solver metrics on noisy artificial growing data.

 15 The column descriptions are the same as the one in Table14. 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 L corresponds to the number of noisy links.

	|V|	|E| |SC1| |SC2| |IR| |L|	Time	Opt. % Gap Nodes	Iter.
	186	372	20	20	20 152	1.01	39.00	0.00	1
			26	24	43	30	5.16	39.00
	280	688	20	20	40 212	6.89	79.00	0.00	1
			28	23	89	50	12.26	79.00		19
	366	946	20	20	60 262	42.52 123.50	3.64	1	16
			25	26 132	60	79.28 119.00		50
	452 1208	20	20	80 320	90.06 161.50	1.55	1	27
			22	23 181	78	295.22 159.00		129
	556 1366	20	20 100 322	196.04 199.00	0.00	1	42
			23	24 231	40	244.41 199.00		92
	662 1804	20	20 120 412 1007.59 242.50	1.44	1	84
			29	26 276	90 1434.16 239.00		210
	736 1946	20	20 140 454 1108.09 283.00	1.41	1 228
			24	26 318	92 3540.79 279.00		691
	822 2212	20	20 160 514 1118.76 323.00	1.24	1	86
			26	27 358 112 2449.55 319.00		287
	902 2362	20	20 180 542 1591.18 363.00	1.10	1	91
			26	27 398 100 2576.60 359.00		269
	996 2656	20	20 200 602 2294.85 404.00	1.24	1 116
			26	26 446 120 3747.02 399.00		351

Table 16 -Sequence and Quast metrics for the initial synthetic data version.

 16

Table 17 -Sequence and Quast metrics for the modified synthetic data version.

 17 The caption is the same as in Table16.

	Successful	Quast

Table 18 -Benchmark 3 v1 contig Quast

 18

	Instance	ILPs	%gnm	#mis
	Abies_alba	dr-sc	99.32	0	0
		ir-sc	99.32	0	0
	Acorus_americanus	ir-sc	99.52	0	0
		(the table continues on the next page)
		A17

Table 19 -Benchmark 3 v1 ILP stats

 19

	Table 18, continued
	Instance	ILPs	%gnm	#mis
	Agathis_dammara	dr-ir-sc	80.19	0	0
		ir-dr-sc	80.19	0	0
	Azima_tetracantha	ir-sc	99.93	0	0
	Begonia_pulchrifolia	-		---
	Carpodetus_serratus	ir-sc	98.61	0	0
	Circaeaster_agrestis	ir-sc	99.56	0	0
	Clematis_repens	ir-sc	99.93	0	0
	Commiphora_foliacea	ir-sc	98.34	0	0
	Cucumis_hystrix	ir-sc	99.49	0	0
	Eucommia_ulmoides	ir-sc	99.01	0	0
	Jasminum_tortuosum	ir-sc	99.77	0	0
	Juniperus_scopulorum	sc	98.82	0	0
	Lamprocapnos_spectabilis dr-sc	35.11	0	0
	Lathyrus_pubescens	dr-sc	98.20	0	0
		ir-sc	98.20	0	0
	Lophocereus_schottii	sc	98.63	0	0
	Musa_ornata	ir-sc	99.08	0	0
	Oenothera_glazioviana	ir-sc	98.50	0	0
	Pelargonium_nanum	ir-sc	96.04	0	0
	Podocarpus_totara	sc	98.76	0	0
	Porphyra_purpura	dr-sc	99.95	0	0
	Sagittaria_trifolia	ir-sc	98.55	0	0
	Sciadopitys_verticillata	dr-sc	96.24	0	0
		ir-sc	96.24	0	0
	Sciaphila_densiflora	sc	99.87	0	0
	Selaginella_kraussiana	dr-sc	99.84	0	0
	Selaginella_vardei	dr-sc	99.98	0	0
	Taxus_baccata	sc	98.07	0	0
	Triosteum_pinnatifidum	ir-dr-sc	98.30	0	0
	Uvaria_macrophylla	ir-sc	98.35	0	0
	Welwitschia_mirabilis	ir-sc	99.17	0	0
	Wolffia_australiana	ir-sc	99.83	0	0
	Instance	|C| |L| |V|	|E| Time
	Abies_alba	9 28 20	36		0.02
						0.02
	Acorus_americanus	5 16 16	40		0.04
						0.05
	Agathis_dammara	20 54 50	96		0.20
						0.64
	Azima_tetracantha	7 20 16	28		0.00
	(the table continues on the next page)
		A18

Table 20 -Benchmark 3 v2 ILP stats

 20

	Instance	|C| |L| |V|	|E| Time
	Agathis_dammara	20 54 52 108	0.22
				0.83
	Begonia_pulchrifolia	6 14 14	22	0.00
				0.00
	Carpodetus_serratus	13 42 36	74	0.13
				0.20
	Jasminum_tortuosum	10 28 30	70	0.08
				0.13
	Lamprocapnos_spectabilis 12 36 48 180	0.68
				2.67
	Lathyrus_pubescens	17 44 36	52	0.02
				0.02
	Lophocereus_schottii	12 34 24	34	0.00
				0.00
	Pelargonium_nanum	11 32 36 128	0.33
				0.84
	Podocarpus_totara	13 40 26	40	0.00
				0.00
	Triosteum_pinnatifidum	13 40 32	66	0.07
				0.18

  84, see also forward & {f, r} F f Set of forward fragments 27, 61-63, 65, 67, 71 F r Set of reverse fragments 27, 61, 63, 65, 67, 71 F Set of fragments 27

  , A16, see also reverse Σ F Fragment label set 27, 61, 62, 64, 70 Σ Fun Unoriented fragment label set 27, 61, 66, 68 Σ L Link label set 28 Σ Lcan Canonical link label set 61, 62, 64, 66, 68, 70 G3 SYMBOLS F un Set of unoriented fragments 27, 33 Function diradj Adjacent direct fragment function 103, 104, 107, 110, 115, A11, see also repadj & invadj dirf rag Direct fragment function 101, 107, 114-116, A1, see also repf rag & invf rag invadj Adjacent inverted fragment function 103, 104, 107, 110, 115, A2, see also repadj & diradj invf rag Inverted fragment function 101, 107, 114-116, A1, see also repf rag & dirf rag repadj Adjacent repeated fragment function 104, 107, 109, see also dirf rag & invf rag repf rag Repeated fragment function 101, 107, 109, see also dirf rag & invf rag Graph G Generic graph see also undirected graph, directed graph, bidirected graph, vertex & edge DG Directed fragment multigraph, DG = (V , E, Φ) 29-31, 40, 41, 50, 57, 61, 99, see also directed graph BG Bidirected fragment multigraph, BG = (V , E, Φ, attre) 32-35, 40, 41, 46, 51, 52, 68, see also bidirected graph U G Undirected fragment multigraph, U G = (V , E F , E L , Φ L ) 35-

  see also vertex V f Forward vertex set 66, see also vertex & forward E Edge set 30-35, 48, 61-68, G4, see also edge E F Fragment-edge set 36-38, 70, 71, G4, see also edge, fragment & U G E L Link-edge set 36-38, 70, G4, see also edge, link & U GΦ Incidence function30-35, 61, 63, 65, 67, 68, 116-118, G4, see also edge G4 SYMBOLS Φ L Incidence function for link-edges 36-38, 70, G4, see also edge, link & U G N Neighbour set 61, 68, 74, 76-79, 82-84, see also edge N -Predecesor set 62, 64, 66, 70, 74, 75, 77-79, 81, 82, 84, see also edge, N & N + N + Successor set 62, 64, 66, 70, 74-84, see also edge, N & N - N F Fragment extremity neighbour set see also edge, N & U GA N L Link neighbour set 70, 71, see also edge, N & U GA Σ V Vertex label set 62, 64, 66, 68, 70, 72, 75-77, see also vertex Σ V f Forward vertex label set 66, see also vertex Σ E Edge label set 62, 64, 66, 68, 70, 75-77, see also edge DGA All oriented fragments multidigraph 61, 63, 65, see also DG & directed multigraph DGS Oriented fragments' successors multidigraph 63, 65, 70, 72, 74, 86, 88, 136, 137, see also DG & directed multigraph DGF Forward fragments multidigraph 65, 67, 68, 70, 74, 86, 88, 136, 137, see also DG & directed multigraph BGU Unoriented fragments multibigraph 68, 70, 74, 86, 88, 136, 137, see also BG & bidirected multigraph U GA All fragments' extremities undirected graph multiungraph 70, 72, see also U G & undirected multigraph Set {0, 1} Set of binaries 106, 109 N The set of natural (positive) integer 28, 31, 34, 37, 61, 92, 99, 112, 119-121 Σ nuc Nucleotide alphabet set Σ nuc = {A, C, G, T } 2, 27, 92 R The set of real numbers 92, 99, 105, 106, 108, 111

  Graph structure where the vertices are word of length k, and edges are overlaps of length k -1. Can be vertex-centric or edge-centric. 18-21, 26, 55, see also k-mer directed graph Graph where the edges are directed 16, 17, 29, 87, 91, 92, 99, 136 directed multigraph A graph where the same oriented edge between the same two vertices can occur several times 30, 31, 61, 88, 99, 116, 136, see also multigraph & directed graph

	GLOSSARY
	edge Component of a graph that connects two vertices 17
	E
	G9

, see also read assembly & scaffolding D de Bruijn graph

  , see also read assembly & scaffolding fragment-edge Component of a graph that connects two vertices 35-37, 70-72, 88, 136, see also graph & vertex fragment-neighbour Vertex connected to a given one forming a fragment-edge 70, see also neighbour, fragment-edge & U G G genome Entire set of DNA instructions found in a cell 2-10, 13-15, 17, 18, 22, 23, 92, 135, see also DNA graph Object composed of vertices connected by edges 22, 30, 60, 72, 74, 86-88, 135-137, see also vertex & edge K k-mer Word of lenght k 18-21, 54-56, 126, see also read assembly, SBH, DBG & de Bruijn graph

L link Ordered pair of oriented fragment 22

  see also graph & vertex link-neighbour Vertex connected to a given one forming a link-edge 70-72, see also neighbour, link-edge & U G Organelle found in the cells of most eukaryotes specialised in the aerobic respiration 4, 5, 8 multigraph A graph where there can be several edges between the same two vertices, see also graph multimeric form Copies of genome that structurally differ 6, 53 N neighbour Vertex connected to a given one 60, 63, 65, 67-70, 72, 74, 80-82, 88, 136 nucleotide Basic building block of nucleic acids (RNA and DNA) 2, 3, 7-9, 11, 12, 17, 40, see also DNA O organelle Specialised subunit, usually within a cell, that has a specific function 4, 5, 8 overlap Suffix-prefix alignment type between two sequences 9, 13-18, 50, 87, 88, 135, 136, see also alignment & link overlap graph Graph structure that stores overlaps between reads 17, 22, see also overlap P path A sequence of vertices such that two consecutive vertices are connected by an edge in the graph 17, 31, 34, 35, 37, 38, see also walk, graph, vertex & edge predecessor In-vertex of one edge of one given vertex 61-64, 66, 67, 70, 74-77, see also successor R read DNA fragment from one DNA's strand, output by a sequencer 9

	GLOSSARY
	M
	mitochondrion
	G10

  see also forward S scaffold Sequence of oriented contigs separated by gap 21, 39, 135, see also scaffolding scaffolding Orienting and ordering the contigs22, 23, 30, 35, 39, 40, 87, 135, 136, see also fragment assembly, read assembly, contig & scaffold sequencer Sequencing technology machine 9 sequencing Method to generate nucleotide fragments 2, 8-11, 13-15, 26, 53 strand The DNA molecule is made up of two strands, each of which has a complementary sequence to the other 2, 3, 8, 10-12 string graph Graph structure where vertices are genomic regions, and edges overlaps between them 17, 18, see also overlap graph & overlap successor Out-vertex of one edge of one given vertex 61-67, 70, 74, 76, 77, see also predecessor U undirected graph Graph where the edges are undirected 35, 37, 47, 88, 136 undirected multigraph A unoriented graph where there can be several edges between the same two vertices 33, 36, 68, 70, see also multigraph & directed graph unitig Special case of a contig, it represents a non-ambiguous assembled sequence 18, 19, 21, 30, 50, 56, see also contig, read assembly & scaffolding

V vertex Component of a graph 17

  , see also graph & edge W walk A sequence of vertices such that two consecutive vertices are connected by an edge in the graph 37, see also path, graph, vertex & edge G12 Titre : Assemblage de fragments ADN : structures de graphes et échaffaudage de génomes de chloroplastes Mot clés : Assemblage de génomes, Programmation linéaire en nombre entiers, Répétitions génomiques, Formes génomiques multimériques Résumé : L'obtention de la séquence nucléotidiques d'une molécule ADN necessite sa fragmentation par des technologie de séquençage, et l'assemblage des fragments. Ces fragments sont appelés lectures. Elles souffrent d'erreurs de séquençage et sont considérées sous deux orientations : celle de leur brin ADN d'origine, ou l'inversecomplémentaire pour l'autre brin. L'assemblage se base sur des chevauchements deux à deux entre des lectures orientées, et est composé de trois phases : l'assemblage des lectures pour obtenir des contigs (des séquences plus longues que les lectures), l'échaffaudage des contigs, pour obtenir des échaffaudages (des ordres de contigs orientés), et la complétion des échaffaudages (trouver les séquences de nucléotides séparant les contigs orientés dans les échaffaudages).

k + 2

https://pypi.org/project/revsymg/

In Chapter II, Section 1, L is a multiset while in this chapter, L is a set.

https://khloraa-scaffolding.readthedocs.io/en/latest

To reproduce the results, please refer to https://khloraascaf-results.readthedocs.i o/en/latest/benchmark_4/.

To reproduce the results, please refer to https://khloraascaf-results.readthedocs.i o/en/latest/benchmark_5/.

For more details, please refer to https://khloraascaf-results.readthedocs.io/en/lat est/benchmark_3/

CpGDB: A Comprehensive Database of Chloroplast Genomes http://www.gndu.ac.in/CpG DB/index.aspx

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/

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_a nalysis/.

For more details, please refer to https://khloraascaf-results.readthedocs.io/en/lat est/benchmark_3/synthetic_data_v2/.

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_analy sis/ and https://khloraascaf-results.readthedocs.io/en/latest/benchmark_3/v2_qua st_analysis/.132

https://pypi.org/project/khloraascaf/

MichelPetrucciani. (1994). Estate[Song]. On Live. Blue Note

https://pypi.org/project/revsymg/

https://github.com/GFA-spec/GFA-spec

https://github.com/lh3/miniasm/blob/master/PAF.md

https://pypi.org/project/khloraascaf/

IV GLOBAL EXACT OPTIMISATIONS FOR CHLOROPLAST GENOME MULTIMERIC FORMS SCAFFOLDING

AInvF be the set of adjacent inverted fragments:

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

Integer Linear Programming (ILP) formulation

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

Definition 5.1: Valid circuit in M DCG

Given a graph M DCG = (V, E) and a starting vertex s, where ctg s is the starting contig, or s = f and occ s = 0. A circuit cp in M DCG is valid if:

-it starts and ends with s;

-consecutive vertices u and v in cp are connected 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 = c∈C mult(c) be a constant, N - v and N + v be the sets of predecessors and successors of vertex v ∈ V , respectively.

return SOR, COR

Multiple genome forms

The sequence of oriented regions SOR represents one chloroplast genome form.

Recall that especially with the LSC-IR-SSC-IR architecture (Figure 6a), the SSC can be reversed during the DNA replication phase. In the following, our goal is to retrieve other forms from the one obtained by the hierarchical problem succession. Towards this goal, we introduce a specific assembly graph: the region graph. The discovery of multiple genome forms is associated with the search of Eulerian circuits in this graph. Figure 38 illustrates the procedure for the toy example.

Definition 8.1: Region graph RegGraph

Given the n-uplet COR of forward regions and the sequence SOR of oriented regions, RegGraph = (V reg, Ereg, Φ) denotes a multidigraph named the region graph, such that: We can easily verify that the number of Eulerian circuits is bounded by O(2 m ), where m ≤ m is the number of inverted repeats. Now, given a region graph RegGraph, finding all the Eulerian circuits is equivalent to retrieve all the possible chloroplast genome forms. Each Eulerian circuit traverses exactly the same regions, but not necessary in the same orientations. Figure 38 gives the resulting region graph obtained from the input data given in Table 11.

We accept all the Eulerian circuits, although they may contradict the repeated region interval constraints given in Definitions 3.2 and 3.3. For example, the oriented region sequence

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

In this chapter

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

APPENDIX

In this chapter

APPENDIX

Two combinatorial reductions are necessary -commutative and pairing reductions: i. commutative reduction means that ∀ (i, j) ∈ RepF , (j, i) / ∈ 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:

Proposition 2.2: 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 distance between two direct fragments such that they can both be chosen to form DRs. Let (i, j) ∈ DirF be a direct fragment and suppose that it is chosen to participate into a DR.

-the occurrences of k and l must differ from the ones of i and j, as a vertex can occupy at most one position, so occ j < occ k ;

-according to Proposition 2.1, occ k < occ l and since occurrences are integers, occ k = occ l -1.

Proposition 2.3

Direct fragment set DirF is a minimal set for DRP.

Proof

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

Reminder of the InvF definition:

A4 Proposition 2.5: Pairing reduction for InvF

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

Proof

The idea is to find the minimum occurrence distance between two inverted fragments such that they can both be chosen to form IRs. Let (i, j) ∈ InvF be an inverted fragment and suppose that it is chosen to participate into an

or l be a candidate for the next direct fragment:

-to respect the property given by Proposition 2.4, or k = 1 -or l = 0

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

-according to Proposition 2.4, occ k < occ l and since occurrences are integers, occ k = occ l -1.

Proposition 2.6

Inverted fragment set InvF is a minimal set for IRP.

Proof

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

Pairs of repeated fragment set reductions

Here we consider the reduction operations for P DirF and P InvF , and conclude on their minimality. Just the commutative reduction is necessary.

APPENDIX

Reminder of pairs of inverted fragments set definition:

Proposition 2.9: Commutative reduction for P InvF ∀ ((i, j), (k, l)) ∈ P InvF , ((k, l), (i, j)) / ∈ P InvF

Proof

Let ((i, j), (k, l)) ∈ P InvF :

-if ctg j < ctg k thus ctg l > ctg i so ((k, l), (i, j)) / ∈ P InvF -else ctg j = ctg k ∧ occ j < occ k thus occ l > occ i so ((k, l), (i, j)) / ∈ P InvF Proposition 2.10 Pair of inverted fragment set P InvF is a minimal set for IRP.

Proof

By reductio ad absurdum: if P InvF can be minimised to a P InvF set, then there exists a pair of inverted fragments (p, q) ∈ P InvF | (p, q) / ∈ P InvF . Let be p = (i, j) and q = (k, l), and let q 2 = (m, n) | (p, q 2 ) ∈ P InvF ∧ (q, q 2 ) ∈ P InvF . i, j, k, l, m, n vertices can simultaneously participate in the solution.

Let build a counter-example: suppose that we have the order i k j l m n. 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) / ∈ P InvF . 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.

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:

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

Proposition 2.11: Canonical reductions for ADirF for different identifiers

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

A8

APPENDIX

Proof

The reciprocal are trivial as ADirF ⊂ E. All the edges using only i, i, k or k belong to ADirF as all the vertex' occurrences are even. While all the edges using j, j, l or l cannot belong to ADirF as one of the vertex' occurrence is odd.

Proposition 2.12: 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.

Proof

Let u 1 , u 2 , . . . , u n , n = 2 u mult 2 be same identifiers vertices participating in a solution path for IRP, in this order (i.e. ∀ k ∈ 1, n , u k is before u k+1 in the solution path).

It can be proven a that: regardless the equation set of k -1 order relation R k,k+1 ∈ {<, >} between two vertices u kid and u k+1id (i.e. u kid R k,k+1 u k+1id , ∀ k ∈ 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 defined the following equation set (used in ADirF ):

Have we got enough edges to build a path with n 2 adjacent direct fragments? On the one hand, the maximum number of edges we would need is: 2 u mult 2 -1 (the 2 is caused by direct fragments). On the other hand, ADirF provides:

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

a

With Lucas Robidou we are currently writing a paper and an associated algorithm Proposition 2.14: Canonical reductions for ADirF for equal identifiers

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

Proof

The reciprocals are trivial as ADirF ⊂ E. Only (i, k) ∈ ADirF and (i, k) ∈ ADirF as the occurrences of i, i, k or k are even, and occ i = occ i < occ k = occ 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:

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

Proposition 2.15: Canonical reductions for AInvF for different identifiers

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

Proof

The reciprocals are trivial as AInvF ⊂ E. Let (i, j) ∈ InvF , (k, l) ∈ InvF be two adjacent inverted fragments such that ctg i < ctg k . We will focus on demonstrating the case (i, k) ∈ E. The others follow the same logics.

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 ≤ occ v and occ u ≥ 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 ≤ 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 d mult = 2, and (c f , d f ) ∈ L and (d f , c f ) ∈ 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 q p 2 can participate into a solution path P ath. In such a case, P ath is sub-defined by edges i → k → i and j → l → j.

According to AInvF set definition, (i, k) ∈ AInvF and (j , l) ∈ AInvF . If we constrain the occurrences by the relation occ u ≤ occ v , thus on the one hand, occ i ≤ occ k (= occ l -1), and on the other hand occ l -1 ≥ occ j -1. As p 1 = p 2 , occ i = occ j -1:

Proposition 2.17: Occurences reduction for equal identifiers in AInvF

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

Proof

Let u 1 , u 2 , . . . , u n , n = 2 u mult 2 be same identifiers vertices participating in a solution path for IRP, in this order (i.e. ∀ k ∈ 1, n , u k is before u k+1 in the solution path).

It can be proven a that: regardless the equation set of k -1 order relation R k,k+1 ∈ {<, >} between two vertices u kid and u k+1id (i.e. u kid R k,k+1 u k+1id , ∀ k ∈ 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 ):

Have we got enough edges to build a path with n 2 adjacent inverted fragments? On the one hand, the maximum number of edges we would need is: 2 u mult 2 -1 (the 2 is caused by inverted fragments). On the other hand, AInvF provides:

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

a

With Lucas Robidou we are currently writing a paper and an associated algorithm Proposition 2.18: Canonical reductions for AInvF for equal identifiers

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

APPENDIX