# DNA fragment assembly, graph structures and chloroplast genome scaffolding: comparative analyses, formulations and implementations

Victor Epain

► To cite this version:

Victor Epain. DNA fragment assembly, graph structures and chloroplast genome scaffolding: comparative analyses, formulations and implementations. Other [cs.OH]. Université de Rennes, 2023. English. NNT: 2023URENS083. tel-04516398

HAL Id: tel-04516398

https://theses.hal.science/tel-04516398

Submitted on 22 Mar 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# THÈSE DE DOCTORAT DE

L'UNIVERSITÉ DE RENNES

ÉCOLE DOCTORALE N° 601

Mathématiques, Télécommunications, Informatique, Signal, Systèmes,

Électronique

Spécialité : Informatique

Par

Victor EPAIN

## Assemblage de fragments ADN : structures de graphes et écha-faudage de génomes de chloroplastes

Analyses comparatives, formulations et implémentations

Thèse présentée et soutenue à Centre Inria de l'Université de Rennes, le 27 novembre 2023

Unité de recherche : Inria / IRISA UMR 6074

## Rapporteurs avant soutenance :

Éric ANGEL Professeur des universités, IBISC, Université Paris-Saclay

Annie CHÂTEAU Maitresse de conférence - HDR, LIRMM, Université de Montpellier

## Composition du Jury :

Président : Élisa FROMONT

Examinateurs :
* Éric ANGEL Professeur des universités, IBISC, Université Paris-Saclay
* Annie CHÂTEAU Maitresse de conférence - HDR, LIRMM, Université de Montpellier
* Élisa FROMONT Professeure des universités, IRISA, Université de Rennes
* Camille MARCHET Chargée de recherche CNRS, CRISTAL, Lille
* Mathias WELLER Professor, Institut für Softwaretechnik und Theoretische Informatik, Berlin
* Dominique LAVENIER Directeur de recherche CNRS, IRISA, Rennes

Dir. de thèse : Rumen ANDONOV Professeur des universités, IRISA, Université de Rennes

Co-dir. de thèse : Jean-François GIBRAT Directeur de recherche, INRAe, Université Paris-Saclay
<figure><img src="image_2.png" /></figure>
# REMERCIEMENTS

En fonction des personnes que je remercie, ces mots sont écrits en français ou en anglais / Depending on whom I am addressing, I am writing these acknowledgements in French or in English.

Merci à Annie Château et Éric Angel qui ont évalué et validé ce manuscrit pour
me permettre de soutenir. Vos remarques positives et de fond m'ont encouragé pour
l'oral et m'ont permis de bien préparer les questions. Merci à Elisa Fromont d'avoir
présidé le jury, composé des rapporteurices, de Camille Marchet et de Mathias
Weller. Je les remercie tou·te·s de s'être déplacé·es à Rennes pour assister en
présence à ma soutenance, en particulier Annie qui venait de Montpellier, Mathias
de Berlin, Camille de Lille et Éric de Paris. Nos échanges ont été plaisants, et je
peux affirmer que ma soutenance a été un très beau moment de ma vie.

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

Lors de ma dernière année de prépa math-
physique au Lycée Dupuy de Lôme de Lorient (2017),
j'ai assisté à une conférence donnée par Pierre-Henri
Gouyon pour introduire son livre Le Fil de la vie, La
face immatérielle du vivant, coécrit avec Jean-Louis
Dessalles et Cédric Gaucherel. Leur livre illustrait
entre autres comment les mathématiques et les mo-
dèles de la physique du signal permettaient d'appré-
hender des logiques et des tendances biologiques. À
l'issue de la conférence, P-H Gouyon m'encourageait à revenir vers la biologie (voir
son autographe en figure ci-contre).

J'ai ainsi souhaité continuer en master bioinformatique à Rennes, et la transition
était toute trouvée : une licence d'informatique à l'Université de Rennes 1. Rumen
Andonov m'a d'abord accueilli pour un stage de deux mois. À ce moment-là,
Sébastien Letort me poussait à la rigueur même si je programmais en **Python**,
chose que j'ai gardé, merci à toi. Rumen et Dominique Lavenier m'ont de nouveau
intégré dans l'équipe GenScale pour le stage de M1. Cette fois-ci, nous nous sommes
envolés pour Los Alamos, première fois que je prenais l'avion, première fois que je
traversais l'Atlantique. Merci Rumen pour ce mois de cohabitation, de travail et de
voyages. Thank you Hristo Djidjev for having welcomed us into your laboratory
(LANL) and for the valuable discussions about graph partitioning. Jamais deux
sans trois, mon dernier stage d'étude s'est fait... à GenScale!
Rumen et Dominique m'ont alors présenté Jean-François Gibrat avec qui nous avons préparé le sujet d'une thèse que nous souhaitions soumettre à concours. À l'origine, le titre était : *Développement de méthodes efficaces, précises et conviviales pour corriger, assembler et aligner des lectures issues des technologies de séquençage troisième génération*.

Merci à Olivier Dameron, Annabelle Monnier et Sophie Schbat d'avoir soutenu ma candidature. Merci aux membres du CSID, Guillaume Fertin et Christophe Klopp, pour les discussions scientifiques et leur vérification en faveur du bon déroulement de la thèse. Rumen, tu m'as toujours dit : « l'important, c'est que tu te fasses plaisir dans ta recherche ». Merci infiniment à toi, Dominique et Jean-François pour votre confiance, les très nombreuses discussions scientifiques, et de m'avoir permis d'exercer mes premiers pas de chercheur librement.

Durant mes deux premières années de thèse, j'ai eu l'occasion et la chance d'enseigner dans les matières qui me tenaient à cœur : méthodes algorithmiques sur les graphes et programmation linéaire. Merci Rumen, Kerian Thuillier et Arthur Gontier d'avoir été une super équipe pédagogique. Merci aux étudiant·es pour leurs retours alors que l'on enseignait en hybride, puis à distance confiné·es : leurs messages m'ont confirmé ma volonté d'enseigner. Merci à Pauline Hamon-Giraud pour son travail à nos côtés en stage de M1.

Thank you, Gunnar Klau and Sven Schrinner from the HHU ALBI team in Düsseldorf, for welcoming me during two months while we were threatened by the lockdown due to COVID-19, for your valuable discussions about integer linear programming and problem complexity.

Merci beaucoup Nicolas Buton, Garance Gourdel, Matthieu Bouguéon et Lucas Robidou d'avoir formé l'équipe de doctorant·es que nous avons été : nous garderons de très bons souvenirs de la réalisation de notre court métrage pour Sciences en Cour[t]s, de l'organisation du même festival l'année suivante, de la gestion de l'association Nicomaque. Merci pour votre soutien au quotidien.

Merci à Clara Delahaye, Olivier Dennler, Grégoire Siekaniec et Téo Lemane, nos ainé·es, de nous avoir accompagné·es dans la recherche, dans Sciences en Cour[t]s et Nicomaque.

Elisa Fromont, merci pour ta confiance et de m'avoir laissé la liberté de ton pour l'organisation des séminaires annuels DKM 2023. Un grand merci également à Marie Farge pour nos discussions concernant l'éthique de la publication scientifique
et tes encouragements pour la suite.

Merci Ambre Ayat pour ton aide dans l'organisation et l'animation de l'atelier éthique dans l'ESR. Merci Bernard Friot et Enka Blanchard pour les échanges autour du productivisme dans la recherche et les moyens d'en sortir.

To my (the best) office colleagues, Nicolas Guillaudeux, Khodor Annoush and Arya Kaul: thank you for the discussions and your cheerfulness.

Un merci chaleureux aux personnes des équipes Symbioses, pour les discussions sérieuses aux rires béats – parfois les deux en même temps, à table comme aux pauses : merci à Samuel, François C., Konogan, Meven, Pierre, Stéphanie, Catherine, Émeline, Claire, Emmanuelle, Jacques, Jeanne, Roland, Sandra, Baptiste, Victor, Camille, Riccardo, Olivier B., Sarah, Thomas R., Thomas C., Gildas, Gaëtan, Grégoire P., François M., Karel, Julien, Florian. Merci Marie et Gaëlle pour votre travail de support quotidien à la recherche et les rigolades. Laurence, Sawako, Charles, Anne-Claire et Maeva de la cafétéria de l'Inria, et à l'équipe du restaurant de CentraleSupélec : merci pour votre bonne humeur de tous les jours.

Merci à mes amiès de Lorient et à ma famille pour leur soutien, leurs encouragements et leur patience quand je m'exerçais à la vulgarisation de sujet de thèse.

Merci aux amiès du master bioinformatique pour cette amitié que nous avons su garder malgré les éloignements, et pour la plupart, malgré les thèses que vous prépariez aussi. Je souhaite que l'on puisse exercer nos recherches en assumant nos responsabilités éthiques.
<figure><img src="image_4.png" /></figure>
# RÉSUMÉ EN FRANÇAIS

## Séquençage ADN et Assemblage de génomes

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 21e 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 :
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.

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

## 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). Kececioglu (1991) l'employa pour la première fois. Les
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 *haplotypes structuraux* réversibles coexistent dans un même chloroplaste.

Plusieurs approches ont déjà été proposées pour l'assemblage spécifique d'orga-
nites 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 (Ankenbrand et al.,
2018). 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 (Coissac et al., 2016; Dierckxsens et al., 2017). Jin et al.
(2020) estiment les multiplicités des contigs afin qu'elles soient 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
Andonov et al. (2019) (un échafaudeur appliqué sur des données chloroplastiques)
tiennent compte de solutions multiples correspondant à des haplotypes structuraux.
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 *NPN*-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/).
<figure><img src="image_6.png" /></figure>

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. Le cout de stockage en mémoire et le coût temporel des algorithmes d'itération sur les données, d'ajout ou de suppression d'informations sont calculés pour chacune des implémentations. 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 dont la taille permet 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 d'haplotypes structuraux. 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
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 sou-
haitons 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 présence dans les cellules de plantes. 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.
<figure><img src="image_7.png" /></figure>

# CONTENTS

Remerciements i

Résumé en français v

Séquençage ADN et Assemblage de génomes. v

Graphes de fragments ADN. vi

Échafaudage de génomes de chloroplastes. vii

Conclusions et perspectives ix

Contents xi

List of Figures xvii

List of Tables xix

I Introduction 1

1 DNA overview. 2

1.1 A recipe for cellular mechanisms 3

1.2 On the motivations for knowing the whole DNA sequence. 3

1.3 Various molecular conformations structure the genome 4

2 Chloroplast genome. 4

2.1 Genome division 5

2.2 Repeats and single-copy regions. 5

2.3 Structural haplotypes. 7

2.4 Genome evolution 8
CONTENTS

<table><tr><td>3</td><td>DNA sequencing</td><td>8</td></tr><tr><td>3.1</td><td>The Shotgun sequencing approach</td><td>9</td></tr><tr><td>3.2</td><td>First generation: Sanger and BAC technologies</td><td>9</td></tr><tr><td>3.3</td><td>Second generation: high-throughput sequencing</td><td>10</td></tr><tr><td>3.4</td><td>Third generation: single molecule run time</td><td>11</td></tr><tr><td>3.5</td><td>Supplementary sequencing data</td><td>12</td></tr></table>

<table><tr><td>4</td><td>The challenges of fragment assembly</td><td>13</td></tr><tr><td>4.1</td><td>Read sequence alignment</td><td>13</td></tr><tr><td>4.2</td><td>Unknown fragment orientations</td><td>14</td></tr><tr><td>4.3</td><td>Sequence similarity: single-copy or repeat?</td><td>14</td></tr><tr><td>4.4</td><td>True sequence divergences or sequencing errors?</td><td>15</td></tr></table>

<table><tr><td>5</td><td>Fragment assembly approaches</td><td>15</td></tr><tr><td>5.1</td><td>The Shortest Common Superstring</td><td>15</td></tr><tr><td>5.2</td><td>Overlap-Layout-Consensus</td><td>17</td></tr><tr><td>5.3</td><td>De Bruijn Graph approach</td><td>18</td></tr><tr><td>5.4</td><td>Breaking down the fragment assembly problem</td><td>20</td></tr></table>

<table><tr><td>6</td><td>Addressed research topics</td><td>21</td></tr><tr><td>6.1</td><td>Graph structure for read assembly and scaffolding stages</td><td>22</td></tr><tr><td>6.2</td><td>Scaffolding of chloroplast structural haplotypes</td><td>22</td></tr></table>

II State-of-the-art 23

<table><tr><td>1</td><td>Graph structure for fragment assembly</td><td>24</td></tr><tr><td>1.1</td><td>Notations and fundamental definitions</td><td>24</td></tr><tr><td>1.2</td><td>Directed graph (DG): oriented fragments based</td><td>27</td></tr><tr><td>1.3</td><td>Bidirected graph (BG): oriented walk based</td><td>30</td></tr><tr><td>1.4</td><td>Undirected graph (UG): tail-head fragments based.</td><td>33</td></tr></table>

<table><tr><td>2</td><td>Scaffolding the contigs.</td><td>37</td></tr><tr><td>2.1</td><td>Scaffolding input data</td><td>38</td></tr><tr><td>2.2</td><td>Subsampling the input data.</td><td>40</td></tr><tr><td>2.3</td><td>Orienting the contigs</td><td>44</td></tr><tr><td>2.4</td><td>Ordering the oriented contigs</td><td>46</td></tr><tr><td>2.5</td><td>Orienting and ordering the contigs simultaneously</td><td>47</td></tr><tr><td>2.6</td><td>Solving approaches</td><td>49</td></tr></table>

<table><tr><td>3</td><td>Chloroplast genome assembly</td><td>51</td></tr><tr><td>3.1</td><td>Chloroplast sequence extraction.</td><td>51</td></tr><tr><td>3.2</td><td>Chloroplast reads assembly</td><td>53</td></tr></table>
CONTENTS

<table><tr><td>3.3 Chloroplast scaffolding</td><td>54</td></tr><tr><td>3.4 Chloroplast assembly validation</td><td>55</td></tr></table>

III Fragment graph implementations and comparison 57

1 Implementations 58

1.1 Directed graph (DG): oriented fragments based 59

1.2 Bidirected graph (BG): oriented walk based 66

1.3 Undirected graph (UG): tail-head fragments based. 68

1.4 Fragment graph map 70

2 Algorithms for DGS, DGF and BGU 72

2.1 Subfunctions 72

2.2 Iterating over the predecessors 73

2.3 Iterating over the successors 74

2.4 Adding a vertex 75

2.5 Adding an edge. 76

2.6 Deleting a vertex 78

2.7 Deleting an edge 81

3 Time costs. 82

4 Memory and time cost comparisons 84

5 Conclusions and perspectives 85

IV Global exact optimisations for chloroplast structural haplotype scaffolding 87

1 Introduction 88

1.1 Chloroplast genome specificities. 88

1.2 State-of-the-art. 89

1.3 Our approach. 89

2 Input data and notation 90

2.1 Set of contigs C 90

2.2 Set of links L 91

2.3 Mathematically defining genomic regions 91

3 Chloroplast scaffolding problem formulations 92

4 Graph and repeated fragment sets 96

4.1 Graph structure 97
CONTENTS

<table><tbody><tr><td>4.2</td><td>Repeated fragment sets.</td><td>98</td></tr><tr><td>5</td><td>Integer Linear Programming (ILP) formulation.</td><td>102</td></tr><tr><td>5.1</td><td>Circuit constraints</td><td>103</td></tr><tr><td>5.2</td><td>Repeated regions constraints</td><td>104</td></tr><tr><td>5.3</td><td>Fixing regions constraints.</td><td>107</td></tr><tr><td>5.4</td><td>Speed-up constraints</td><td>108</td></tr><tr><td>5.5</td><td>Scaffolding problems ILP</td><td>108</td></tr><tr><td>6</td><td>Hierarchical problem succession.</td><td>109</td></tr><tr><td>7</td><td>From an ILP solution to a genome structure</td><td>110</td></tr><tr><td>8</td><td>Multiple genome forms.</td><td>114</td></tr><tr><td>9</td><td>NP-completeness</td><td>117</td></tr><tr><td>10</td><td>Numerical results</td><td>121</td></tr><tr><td>10.1</td><td>Complexity validation on artificial data</td><td>122</td></tr><tr><td>10.2</td><td>Synthetic chloroplast input data.</td><td>123</td></tr><tr><td>11</td><td>Conclusion.</td><td>130</td></tr><tr><td>12</td><td>Discussion and perspectives.</td><td>131</td></tr><tr><td>V</td><td>Conclusions and perspectives</td><td>133</td></tr><tr><td>Fragment graph.</td><td></td><td>134</td></tr><tr><td></td><td>Thesis contribution.</td><td>134</td></tr><tr><td></td><td>Short-term future work.</td><td>135</td></tr><tr><td></td><td>Long-term future work</td><td>135</td></tr><tr><td>Chloroplast genome scaffolding.</td><td></td><td>136</td></tr><tr><td></td><td>Thesis contribution.</td><td>136</td></tr><tr><td></td><td>Short-term future work.</td><td>137</td></tr><tr><td></td><td>Long-term future work</td><td>138</td></tr><tr><td>Bibliography</td><td></td><td>139</td></tr><tr><td>Appendix</td><td></td><td>A1</td></tr><tr><td>1</td><td>Repeated fragment set functions</td><td>A1</td></tr></tbody></table>
CONTENTS

<table><tr><td>2</td><td>Reduction of the repeated fragment sets</td><td>A2</td></tr><tr><td>2.1</td><td>Repeated fragment set reductions.</td><td>A2</td></tr><tr><td>2.2</td><td>Pairs of repeated fragment set reductions.</td><td>A5</td></tr><tr><td>2.3</td><td>Adjacent repeated fragment set reductions</td><td>A8</td></tr><tr><td>3</td><td>Metrics.</td><td>A17</td></tr><tr><td>3.1</td><td>Quast metrics</td><td>A17</td></tr><tr><td>4</td><td>Supplementary results</td><td>A18</td></tr><tr><td>4.1</td><td>v1 scaffolding benchmark.</td><td>A18</td></tr><tr><td>4.2</td><td>v2 scaffolding benchmark.</td><td>A20</td></tr><tr><td>Acronyms</td><td></td><td>G1</td></tr><tr><td>Symbols</td><td></td><td>G3</td></tr><tr><td>Computational terms</td><td></td><td>G7</td></tr><tr><td>Glossary</td><td></td><td>G11</td></tr></table>
<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

# LIST OF FIGURES

## I Introduction 2

1 DNA double-strand molecule 3

2 Different DNA conformations 4

3 Simplified DNA distribution in a plant's cell. 5

4 DNA repeats. 6

5 Common chloroplast genome's architectures. 7

6 Inverted repeats causing structural haplotypes: flip-flop inversion during the DNA replication. 8

7 The shotgun sequencing approach. 10

8 Sequencing error types. 11

9 Paired-end sequencing. 12

10 Sequence alignments. 14

11 Unknown fragment's orientation. 15

12 Impact of the genome repeats on the fragment sequencing. 16

13 OLC versusDBG fragment assembly approaches. 19

## II State-of-the-art 26

14 Examples of links 26

15 Fragment and link sets with functions 27

16 Link cases in DG. 28

17 A path in DG. 29

18 Link cases in BG. 31

19 A path in BG. 33

20 Link cases in UG. 35

21 A path in UG. 36

22 Overview of all the fragment graph structures. 37
LIST OF FIGURES

<table><tr><td>23</td><td>Scaffolding input data.</td><td>41</td></tr><tr><td>24</td><td>Bidirected and undirected graph cycles.</td><td>45</td></tr></table>

III Fragment graph implementations and comparison 59

<table><tr><td>25</td><td>DGA implementation.</td><td>60</td></tr><tr><td>26</td><td>DGS implementation.</td><td>62</td></tr><tr><td>27</td><td>DGF implementation.</td><td>64</td></tr><tr><td>28</td><td>BGU implementation.</td><td>67</td></tr><tr><td>29</td><td>UGA implementation.</td><td>69</td></tr><tr><td>30</td><td>Graph structures and their implementations.</td><td>70</td></tr></table>

IV Global exact optimisations for chloroplast structural haplotype scaffolding 93

<table><tr><td>31</td><td>Repeat degeneration and region orders.</td><td>93</td></tr><tr><td>32</td><td>Chloroplast repeat scaffolding.</td><td>96</td></tr><tr><td>33</td><td><i>MDCG</i> example.</td><td>97</td></tr><tr><td>34</td><td>Repeated fragment sets illustration for two contigs <i>c</i> and <i>d</i>.</td><td>99</td></tr><tr><td>35</td><td>Adjacent repeated fragment sets examples.</td><td>101</td></tr><tr><td>36</td><td>Non-exhaustive illustrations for authorised and forbidden order cases for two repeated fragments ((<i>i</i>, <i>j</i>), (<i>k</i>, <i>l</i>)) ∈ <i>PRepF</i>.</td><td>105</td></tr><tr><td>37</td><td>Extracting the genome architecture in <i>MDCG</i> from a <i>CHSP</i> solution.</td><td>111</td></tr><tr><td>38</td><td>Region graph for the toy example.</td><td>115</td></tr><tr><td>39</td><td>From a digraph <i>G</i> for <i>LPSTP</i> to a digraph <i>G'</i> for <i>IRP</i>.</td><td>120</td></tr><tr><td>40</td><td>Solver running time distributions for perfect artificial data.</td><td>123</td></tr></table>
<figure><img src="image_10.png" /></figure>

# LIST OF TABLES

## II State-of-the-art 38

1 Non-exhaustive categorised list of methods and approaches for the scaffolding stage. 39

2 Scaffolding input data properties. 40

3 Chloroplast genome assembly approaches. 52

## III Fragment graph implementations and comparison 82

4 Calculus detail of basic operations costs. 83

5 Algorithmic costs for subfunctions. 83

6 Algorithmic costs of iterating over the neighbours. 83

7 Algorithmic costs of adding a vertex or an edge. 84

8 Algorithmic costs of deleting a vertex. 84

9 Algorithmic costs of deleting an edge. 84

10 Comparison between DGS, DGF and BGU. 85

## IV Global exact optimisations for chloroplast structural haplotype scaffolding 90

11 Toy example of input data. 90

12 ILP sets and functions corresponding table. 105

13 Problem code combinations 110

14 Gurobi solver metrics on perfect artificial growing data. 124

15 Gurobi solver metrics on noisy artificial growing data. 125

16 Sequence and Quast metrics for the initial synthetic data version. 128
LIST OF TABLES

<table><tr><td>17 Sequence and Quast metrics for the modified synthetic data version.</td><td>130</td></tr><tr><td>Appendix</td><td>A18</td></tr><tr><td>18 Benchmark 3 v1 contig Quast</td><td>A18</td></tr><tr><td>19 Benchmark 3 v1 ILP stats</td><td>A19</td></tr><tr><td>20 Benchmark 3 v2 ILP stats</td><td>A20</td></tr></table>
# I INTRODUCTION

Joe Hisaishi. (1999). *Inner Voyage [Song]. On *The Universe Within*, Vol.1 & 2*. PONY CANYON

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

## In this chapter

1 DNA overview 2

1.1 A recipe for cellular mechanisms 3

1.2 On the motivations for knowing the whole DNA sequence 3

1.3 Various molecular conformations structure the genome 4

2 Chloroplast genome 4

2.1 Genome division 5

2.2 Repeats and single-copy regions. 5

2.3 Structural haplotypes. 7

2.4 Genome evolution 8

3 DNA sequencing 8

3.1 The Shotgun sequencing approach 9

3.2 First generation: Sanger and BAC technologies 9

3.3 Second generation: high-throughput sequencing. 10

3.4 Third generation: single molecule run time 11

3.5 Supplementary sequencing data. 12

4 The challenges of fragment assembly 13

4.1 Read sequence alignment. 13

4.2 Unknown fragment orientations. 14

4.3 Sequence similarity: single-copy or repeat? 14

4.4 True sequence divergences or sequencing errors?. 15

5 Fragment assembly approaches 15

5.1 The Shortest Common Superstring 15

5.2 Overlap-Layout-Consensus 17

5.3 De Bruijn Graph approach 18

5.4 Breaking down the fragment assembly problem 20
<table><tr><td>6</td><td>Addressed research topics</td><td>21</td></tr><tr><td>6.1</td><td>Graph structure for read assembly and scaffolding stages</td><td>22</td></tr><tr><td>6.2</td><td>Scaffolding of chloroplast structural haplotypes</td><td>22</td></tr></table>

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 20th century. They impact many areas, from ecological understanding to economics, through agronomics and medicine.

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.

# 1 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 (Franklin and Gosling, 1953) lead to the double-strand model proposed by Watson and Crick in 1953. 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 $\Sigma_{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 *AATGCCA* is a strand (or a DNA sequence), then *TGGCATT* 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 (Trevors, 1996) to billions for eukaryotes (Kidwell, 2002). 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.
<figure><img src="image_12.png" /></figure>

■ **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,
*AATGCCA* is the reverse of *TGGCATT*, and vice-versa.
This DNA molecule's length is 7 bp.

1.1 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.

1.2 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 Rice and Green (2019) 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.

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

<figure><img src="image_13.png" /><figcaption>■ Figure 2 – Different DNA conformations</figcaption></figure>

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

## 2 Chloroplast genome

Chloroplasts are plant and algae cell organelles. They are the result of endosymbiosis between bacteria and eukaryotic cells (Gould, 2012). 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
1986, the first chloroplast genomes were sequenced and assembled. Those are chloroplasts of a liverwort, *Marchantia polymorpha* (Ohyama et al., 1986), and of a tobacco plant, *Nicotiana tabacum* (Shinozaki et al., 1986). Knowing the sequence of organelle genomes allows performing ecological, phylogenetic and evolutionary studies (Daniell et al., 2016; Sun et al., 2020; Long et al., 2023).

## 2.1 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 (Bendich, 1987; Kobayashi et al., 2002). The number of copies have been shown to decrease with the leaf age (Kumar et al., 2014). Figure 3 schematises the genome distribution in a plant cell.

<figure><img src="image_14.png" /><figcaption>■ Figure 3 – Simplified DNA distribution in a plant's cell.</figcaption></figure>

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.

## 2.2 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);

- 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 4 gives two examples of perfect repeats.

<figure><img src="image_15.png" /><figcaption>■ **Figure 4 – DNA repeats.**</figcaption></figure>

A blue subsequence is the reverse-complement of a green subsequence, i.e. *GCAA* is the reverse of *TTGC*. (a) *GCAA* is direct repeated in the top strand, thus *TTGC* is direct repeated in the bottom strand. (b) *GCAA* is inverted repeated in the top strand, thus *TTGC* is inverted repeated in the bottom strand.

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 (Inverted Repeat (IR)) separated by a Long Single-Copy (LSC) and a Short Single-Copy (SSC) sequences (Bock and Knoop, 2012). Figure 5a illustrates IR-quadripartite chloroplast genome. For example, Thode et al. (2021) 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 Kim and Lee (2005). 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 Direct Repeat (DR), generally shorter than the IRs, can also be found (Palmer, 1985), as well as more complex structure involving both type of repeats (Tsai and Strauss, 1989). Figure 5 summarises the main chloroplast genome structures found in the literature.
<figure><img src="image_16.png" /><figcaption>■ **Figure 5** – Common chloroplast genome's architectures.</figcaption></figure>

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 that illustrates one is the reverse-complement sequence of the other; (b) the two blue DR arrows are in the same direction that illustrates both have the same nucleotide sequence; (c) the two types of repeat can simultaneously exist in the chloroplast genome, and DRs are shorter than IRs.

## 2.3 Structural haplotypes

The molecular form of the copies of chloroplast genomes differ (Bendich, 2004), but the circular form has been mainly studied. For example Seyer et al. (1981) 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 (Bendich, 2004).

For this reason, in this manuscript we focus on the chloroplast circular genome copies. During the DNA molecular replication of one cpDNA molecule, flip-flop inversions can occur so that the resulting replicated molecules' form differs. This transformation is due to the presence of IRs, 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 *structural haplotypes* (Palmer, 1983; Deng et al., 1989; Wang and Lanfear, 2019). Figure 6 schematises a flip-flop
inversion for a chloroplast genome with IRs.

<figure><img src="image_17.png" /><figcaption>■ **Figure 6** – Inverted repeats causing structural haplotypes: flip-flop inversion during the DNA replication.</figcaption></figure>

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 (*heteroplasmy*).

## 2.4 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 (Xiao-Ming et al., 2017). Horizontal transfers (a transfer of genes not by heredity) were also studied from the plant to its organelles.

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

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.

## 3.1 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 Staden (1979) 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 Mehdi et al. (2017) and Pervez et al. (2022).

## 3.2 First generation: Sanger and BAC technologies

In the late 1970s Sanger's research succeeded in sequencing much longer fragments
(Sanger et al., 1977). 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 \times 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.,
1989). Venter et al. (1996) 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.
<figure><img src="image_18.png" /><figcaption>■ Figure 7 – The shotgun sequencing approach.</figcaption></figure>

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.

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.

3.3 Second generation: high-throughput sequencing

As mentioned above, the first decade of the 21st century witnessed the emergence
of new sequencing technologies. Refered to as Next Generation Sequencing tech-
nologies (NGS), they produce millions of short reads at low prices allowing the
sequencing of numerous organisms.

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
ACGTCAT
ACGTCGT
(a)

ACGTCAT
ACGTCAAT
(b)

ACGTCAT
ACGTC_T
(c)

■ **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 latter are categorised as indels.

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.

**454 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 and 8c).

**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, Korbel et al. 2007), they reach the kilobases at the expense of accuracy.

## 3.4 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.
<figure><img src="image_19.png" /><figcaption>■ **Figure 9 – Paired-end sequencing.**</figcaption></figure>

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).

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, Hon et al. 2020).

**Oxford nanopore sequencing (2014)** In Oxford Nanopore Technology (ONT), the DNA fragments pass through 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, Delahaye and Nicolas (2021) measure it to 5% – 8% but remark that the base-caller (the software interpreting ionic signals) is upgrading fast.

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

# 4 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. Shapiro (1967) 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 Gingeras et al. (1979) and Staden (1980). They enable the first shotgun assembly project with the assembly of the 50 kbp virus λ using reads of length 200 bp (Sanger et al., 1982).

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

<figure><img src="image_20.png" /><figcaption>(a) Overlap</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Containment</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) Share</figcaption></figure>

■ **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).

4.2 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.

4.3 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 maize genome or in bacteria genomes Kidwell 2002). 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
<figure><img src="image_23.png" /><figcaption>■ **Figure 11 – Unknown fragment’s orientation.**</figcaption></figure>

$\mathcal{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$.

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.

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

Peltola et al. (1984) tackle the problem of fragment assembly with the Shortest Common Superstring (SCS) algorithm. Given an error ratio $\delta$, each read must participate in the superstring $F$, in only one orientation, modified in sequence with an edit distance less than $\delta$ 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
<figure><img src="image_24.png" /><figcaption>■ **Figure 12 – Impact of the genome repeats on the fragment sequencing.**</figcaption></figure>

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 *abcdefgh*, or *abfdecgh* (if the repeats are not exactly the same, switching *c* and *f* can misassemble the genome), or even *abcgh* that shortens the genome.

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. Räihä and Ukkonen (1981) have previously modelled SCS as the search for a Hamiltonian path in a *directed graph* (*digraph*) that represents overlaps. In that case, the fragment assembly problem is *NP*-hard.

Li (1990) 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, Myers (1995) found that the SCS can wrongly make the genome shorter by merging the repeats.
## 5.2 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*. Kececioglu (1991) 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 Kececioglu (1991) and Kececioglu and Myers (1995). 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, Myers (1995) 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* implemented this approach and enabled to assemble the genome of *Drosophila melanogaster* (Myers et al., 2000).

**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 (Batzoglou et al., 2002) 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 (Miller et al., 2008), the BOG strategy is part of hifiasm (Cheng et al., 2021) which is dedicated to haplotype resolving thanks to HiFi long reads.

**OLC complexity** In Myers (2005), 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, Medvedev et al. (2007) proved the problem to be $NP$-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.

## 5.3 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. Idury and Waterman (1995) associated the genome as a double-Eulerian superwalk in *G*. If the Eulerian superwalk is unique, it should correspond to the genome. The first motivation was the analogy with the search for an Eulerian path, which is simpler 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
<figure><img src="image_25.png" /><figcaption>(a) Read sequencing</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) OLC approach</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(c)DBG approach</figcaption></figure>

■ Figure 13 – OLC versusDBG 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 theDBG 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.

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, Pevzner et al. 2001).
**Local and variant approaches** ManyDBG variant algorithms have implemented their own heuristic for specific targets. Velvet and SOAPdenovo merge bubbles thanks to a Dijkstra-like Breadth-First Search (BFS) algorithm (Zerbino and Birney, 2008; Li et al., 2010). SPAdes applies theDBG approach with different values of k and iteratively merges the results (Bankevich et al., 2012). While the latter software are suitable for short paired-end reads, theDBG 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 (Kolmogorov et al., 2019). 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 (Ruan and Li, 2020). 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 (Ekim et al., 2021).

**DBG complexity** Pevzner (1989) 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 *NP*-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 (Hierholzer and Wiener, 1873). However,DBG corresponds to finding a super-walk. Medvedev and Brudno (2009) associateDBG with the Chinese postman problem (that echoes the SCS approach) on bidirected graphs, and propose an exact polynomial time algorithm. They next linkDBG 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 proveDBG to be in *P*.

5.4 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 align-
ment 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 theDBG 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.

# 6 Addressed research topics

In this thesis we focus mainly on two aspects of genome assembly.
## 6.1 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.

## 6.2 Scaffolding of chloroplast structural haplotypes

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 *NP*-complete. Furthermore, our formulation enables us to tackle the presence of chloroplast structural haplotypes. 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.
<figure><img src="image_28.png" /></figure>

# STATE-OF-THE-ART

Richard Galliano & Miran Vaupotic & Czech Phil-harmonic Quartet. (2023). Petite Suite Française: III. Espiègle [Song]. On *Madreperla*. Navona Records

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

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

## In this chapter

1 Graph structure for fragment assembly 24

1.1 Notations and fundamental definitions 24

1.1.1 Fragment set 24

1.1.2 Link set. 26

1.2 Directed graph (DG): oriented fragments based 27

1.3 Bidirected graph (BG): oriented walk based 30

1.4 Undirected graph (UG): tail-head fragments based. 33

2 Scaffolding the contigs 37

2.1 Scaffolding input data 38

2.2 Subsampling the input data. 40

2.2.1 Bundling the links 42

2.2.2 Removing paired-end reads. 42

2.2.3 Removing contigs 42

2.2.4 Removing links 43

2.2.5 Partitioning the instances. 44

2.3 Orienting the contigs. 44

2.3.1 Maximising the sum of used link weights. 45

2.3.2 Remove the minimum number of contigs and links to avoid odd reversal cycles. 46

2.3.3 Maximum log-likelihood function 46

2.4 Ordering the oriented contigs. 46

2.4.1 Maximising the sum of used distances bunches' weights. 46

2.4.2 Minimum spanning tree 47

2.4.3 Contig positioning only. 47

2.5 Orienting and ordering the contigs simultaneously. 47

2.5.1 Maximising the sum of distance bundles' weights to order linearly the contigs 47

2.5.2 The heaviest matching, paths and cycles. 48
<table><tbody><tr><td>2.5.3</td><td>The maximum (weighted) matching and spanning tree</td><td>48</td></tr><tr><td>2.6</td><td>Solving approaches</td><td>49</td></tr><tr><td>2.6.1</td><td>Greedy approaches</td><td>49</td></tr><tr><td>2.6.2</td><td>Fix parameters dynamic programming</td><td>50</td></tr><tr><td>2.6.3</td><td>Mathematical programming</td><td>50</td></tr><tr><td>3</td><td>Chloroplast genome assembly</td><td>51</td></tr><tr><td>3.1</td><td>Chloroplast sequence extraction</td><td>51</td></tr><tr><td>3.1.1</td><td>Filtering the reads</td><td>52</td></tr><tr><td>3.1.2</td><td>Filtering the contigs</td><td>53</td></tr><tr><td>3.2</td><td>Chloroplast reads assembly</td><td>53</td></tr><tr><td>3.2.1</td><td>De Bruijn graph approach</td><td>53</td></tr><tr><td>3.2.2</td><td>Seed-and-extend</td><td>54</td></tr><tr><td>3.3</td><td>Chloroplast scaffolding</td><td>54</td></tr><tr><td>3.4</td><td>Chloroplast assembly validation</td><td>55</td></tr></tbody></table>

# 1 Graph structure for fragment assembly

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

## 1.1 Notations and fundamental definitions

### 1.1.1 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 $\mathcal{F}_{un}$ the set of unoriented fragments, in bijection with the set of integer label $\Sigma_{\mathcal{F}_{un}} = \llbracket 0, |\mathcal{F}_{un}| \rrbracket$, and denote by $unfid: \mathcal{F}_{un} \hookrightarrow \Sigma_{\mathcal{F}_{un}}$ the labelling function. Let $\Sigma_{nuc} = \{A, T, G, C\}$ be the nucleotide alphabet and $unseq: \mathcal{F}_{un} \to \Sigma_{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 $\bar{\cdot}$ gives the complementary orientation (its reverse), i.e. $\bar{f} = r$ and $\bar{r} = f$.

The set of oriented fragments results from the Cartesian product of $\mathcal{F}_{un}$ and $\{f, r\}$:

► **Definition 1.3: Oriented fragment sets**

Let $\mathcal{F} = \mathcal{F}_{un} \times \{f, r\} = \mathcal{F}_f \sqcup \mathcal{F}_r$ be the set of oriented fragments, associated with the projection functions $\pi_{\mathcal{F}_{un}}: \mathcal{F} \to \mathcal{F}_{un}$ and $for: \mathcal{F} \to \{f, r\}$, $\Sigma_{\mathcal{F}} = \llbracket 0, 2|\mathcal{F}_{un}| \rrbracket$ be the set of fragment integer identifiers and $fid: \mathcal{F} \hookrightarrow \Sigma_{\mathcal{F}}$ be the associated labelling function. Let denote by $seq: \mathcal{F} \to \Sigma_{nuc}^+$ the function that gives the nucleotide sequence of an oriented fragment. It holds that:

- $\mathcal{F}_f$ is the set of forward fragments, such that for each $a_f \in \mathcal{F}_f$:

- $for(a_f) = f;$

- $seq(a_f) = unseq(\pi_{\mathcal{F}_{un}}(a_f)).$

- $\mathcal{F}_r$ is the set of reverse fragments, such that for each $a_r \in \mathcal{F}_r$:

- $for(a_r) = r;$

- $seq(a_r) = \overline{unseq(\pi_{\mathcal{F}_{un}}(a_r))}.$

Where $\bar{\cdot}$ 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 \in \mathcal{F}$, $\bar{a}$ denotes its reverse, where:

$$ \text{— } |fid(\bar{a}) - fid(a)| = 1; $$

$$ \text{— } for(\bar{a}) = \overline{for(a)} \text{ (a and } \bar{a} \text{ have complementary orientations in } \{f, r\}\text{)}; $$

$$ \text{— } seq(\bar{a}) = \overline{seq(a)} \text{ (the nucleotide sequence of } \bar{a} \text{ is the reverse of this of } a\text{).} $$

Figure 15 illustrates the fragment sets and their associated functions.

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

► **Definition 1.5: Link set**

Denote by $\mathcal{L} = (\mathcal{L}', m_{\mathcal{L}})$ the multiset of ordered pairs of fragments, where $\mathcal{L}' \subset \mathcal{F}^2$ is a finite set, and $m_{\mathcal{L}}: \mathcal{L}' \to \mathbb{N}_{>0}$ gives the multiplicity of each element of $\mathcal{L}'$ in $\mathcal{L}$. Also let $linkp: \mathcal{L} \twoheadrightarrow \mathcal{L}'$ be the function that gives the fragment couple associated with a link. Denote by $\Sigma_{\mathcal{L}} = \llbracket 0, |\mathcal{L}| \rrbracket$ the set of link integer labels and its associated function $lid: \mathcal{L} \hookrightarrow \Sigma_{\mathcal{L}}$.

<figure><img src="image_31.png" /><figcaption>■ **Figure 14 – Examples of links**</figcaption></figure>

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.
## ► Property 1.1: Link reverse symmetry

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

$$ \forall (u, v) \in \mathcal{F}^2, m_{\mathcal{L}}((u, v)) = m_{\mathcal{L}}((\bar{v}, \bar{u})) $$

$$ \text{i.e. } (u, v) \in \mathcal{L}' \iff (\bar{v}, \bar{u}) \in \mathcal{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 \in \mathcal{L}$, where $linkp(l) = (u, v)$, $\bar{l}$ denotes its reverse such that $linkp(\bar{l}) = (\bar{v}, \bar{u}) = (\overline{u, v}) = \overline{linkp(l)}$.

<figure><img src="image_32.png" /><figcaption>■ **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.</figcaption></figure>

Figure 15 illustrates 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.

## 1.2 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 Kececioglu (1991).
For the read assembly stage, DG is described in **FALCON** (Chin et al., 2016), **HINGE** (Kamath et al., 2017), **Shasta** (Shafin et al., 2020) and **hifiasm** (Cheng et al., 2021). Concerning the scaffolding stage, GAT refers to an unitig graph (Andonov et al., 2019).

► **Definition 1.7: Oriented fragments based directed multigraph (DG)**

Let $DG = (V, E, \Phi)$ be a directed multigraph (multidigraph) such that:

- $V$ is the set of vertices, where each $v \in V$ represents one of the two orientations for a fragment in $\mathcal{F}$;

- $E$ is the multiset of edges, where each $e \in E$ represents one link in $\mathcal{L}$;

- $\Phi: E \to \{(u, v) \mid u, v \in V, u \neq v\}$ is the incidence function that associates an edge with an ordered pair of vertices.

Denote by *frag*: $V \hookrightarrow \mathcal{F}$ the bijective function that associates each vertex with one and only one fragment, and by *link*: $E \hookrightarrow \mathcal{L}$ the bijective function that associates each edge with one and only one link.

Figure 16 illustrates the link case representations in DG.

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

(a) $(a_f, b_f) \in \mathcal{L}'$

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

(b) $(b_f, a_f) \in \mathcal{L}'$

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

(c) $(a_f, b_r) \in \mathcal{L}'$

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

(d) $(b_r, a_f) \in \mathcal{L}'$

<caption>■ **Figure 16 – Link cases in DG.**</caption>

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.

► **Property 1.2: Sizes of DG**

For $DG = (V, E, \Phi)$, $|V| = |\mathcal{F}|$ and $|E| = |\mathcal{L}|$.

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.
► **Definition 1.8: The reverse operation in DG**

For each vertex $v \in V$ and edge $e \in E$, $\bar{v} \in V$ and $\bar{e} \in E$ denote their reverse. From the bijective properties of functions $frag$ and $link$, it holds that:

$$— frag(\bar{v}) = \overline{frag(v)};$$

$$— link(\bar{e}) = \overline{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, \dots, v_{n-1}) \in V^n, n \in \mathbb{N}$ in $DG = (V, E, \Phi)$ is valid if the following two properties hold:

**Contiguity**

$$\forall i \in [0, n-1], \exists e \in E \mid \Phi(e) = (v_i, v_{i+1})$$

**Exclusive orientation**

$$\forall v \in p, \bar{v} \notin 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><img src="image_37.png" /><figcaption>■ **Figure 17 – A path in DG.**</figcaption></figure>

$p = v_f u_f x_r$ is a valid path coloured in red. $\bar{p} = x_f u_r v_r$ is its reverse and is also valid. However, $p$ and $\bar{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**

$$ \neg \forall a \in \mathcal{F}, \exists! v \in V \mid \text{frag}(v) = a. $$

$$ \neg \forall l \in \mathcal{L}, \exists! e \in E \mid \text{link}(e) = l. $$

a∃! $x \in X$ means "there is a unique $x$ in $X$".

▷ Proof

From Definition 1.7, frag: $V \rightsquigarrow \mathcal{F}$ and link: $E \rightsquigarrow \mathcal{L}$ are bijective.

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

## 1.3 Bidirected graph (BG): oriented walk based

The bidirected graph (bigraph) was introduced by Edmonds and Johnson (1970). Myers (1995) is the first to suggest a bigraph structure (BG) to store the links exploiting the reverse symmetry of the fragment and the link sets $\mathcal{F}$ and $\mathcal{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 Gritsenko et al. (2012) and suggested by Myers (1995).

For the read assembly stage, BG is described in Bambus (Pop et al., 2004), Minimus (Sommer et al., 2007) and in Edena (Hernandez et al., 2008). Concerning the scaffolding stage, BG is found in SOPRA (Dayarian et al., 2010), Bambus 2 (Koren et al., 2011), MIP (Salmela et al., 2011), Opera (Gao et al., 2011), SSPACE (Boetzer et al., 2011), GRASS (Gritsenko et al., 2012), BOSS (Luo et al., 2017), LACHESIS (Burton et al., 2013), ScaffoldScaffolder (Bodily et al., 2016) and SLR (Luo et al., 2019).

► **Definition 1.10: Oriented walk based bidirected multigraph (BG)**

Let $BG = (V, E, \Phi, \text{attre})$ be a bidirected multigraph (multibigraph) such that:

- $V$ is the set of vertices, where each $v \in V$ represents the two orientations for a fragment in $\mathcal{F}$;

- $E$ is the multiset of edges, where each $e \in E$ represents one link and its reverse in $\mathcal{L}$;

- $\Phi: E \to \{\{u, v\} | u, v \in V, u \neq v\}$ is the incidence function that associates an edge with a pair of vertices;
— *attre*: $E \to \{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 *unfrag*: $V \hookrightarrow \mathcal{F}_{un}$ the bijective function that associates each vertex with one and only one unoriented fragment, and by:

$$ \textit{link}: \{(u, v, e) \mid e \in E, \Phi(e) = \{u, v\}\} \hookrightarrow \mathcal{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.

<figure><img src="image_39.png" /><figcaption>■ **Figure 18 – Link cases in BG.**</figcaption></figure>

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 $unfid(a) < unfid(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, \Phi, \textit{attre})$, $|V| = \frac{1}{2}|\mathcal{F}|$ and $|E| = \frac{1}{2}|\mathcal{L}|$.
Less explicitly than in Definition 1.8, Definition 1.11 extends the reverse operations given in Definitions 1.4 and 1.6 for BG.

### ► Definition 1.11: The reverse operation in BG

Let $e \in E$ be an edge, $\Phi(e) = \{u, v\}$, $a = \text{unfrag}(u)$ and $b = \text{unfrag}(v)$. Assume that $\text{unfid}(a) < \text{unfid}(b)$. As function $\text{link}$ is bijective, it holds that:

$$ \text{link}(v, u, e) = \overline{\text{link}(u, v, e)}; $$

$$ \text{attre}(e) = (\text{for}(a), \text{for}(b)) \text{ and } \text{attre}(e) \begin{bmatrix} 0 & - \\ - & 0 \end{bmatrix} = (\overline{\text{for}(b)}, \overline{\text{for}(a)}). $$

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, \dots, v_{n-1}) \in V^n$, $n \in \mathbb{N}$ in BG = $(V, E, \Phi, \text{attre})$ is valid if the following two properties hold:

**Contiguity** $\forall i \in [0, n-1], \exists e_i, e_{i+1} \in E$ s.t.:

$$ \Phi(e_i) = \{v_{i-1}, v_i\} \wedge \Phi(e_{i+1}) = \{v_i, v_{i+1}\} $$

$$ \text{link}(v_{i-1}, v_i, e_i)[1]^a = \text{link}(v_i, v_{i+1}, e_{i+1})[0] $$

**Exclusive orientation** $\forall v_i, v_j \in p, v_i = v_j \implies i = j$

Finally, Property 1.5 defines the reverse symmetry in BG, that inherits from Property 1.1.

### ► Property 1.5: Reverse symmetry in BG

$$ \forall l \in \mathcal{L}, \exists! e \in E, \Phi(e) = \{u, v\} \mid \text{link}(u, v, e) = l \supset \text{link}(v, u, e) = l. $$

### ▷ Proof

From Definition 1.10, $\text{link}: \{(u, v, e) | e \in E, \Phi(e) = \{u, v\}\} \twoheadrightarrow \mathcal{L}$ is bijective.

<footnote>aThe notation $x[i]$ denotes the $i$th element in $x$.</footnote>

<footnote>a$\supset$ is the exclusive disjunction notation.</footnote>
<figure><img src="image_40.png" /><figcaption>(a) bidirected graph view</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(b) undirected graph view</figcaption></figure>

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

⊂

## 1.4 Undirected graph (UG): tail-head fragments based

Huson et al. (2002) first describe the *undirected graph (ungraph)* structure (UG) 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.

The structure mainly appears in method for the scaffolding stage, as in **SCARPA** (Donmez and Brudno, 2013), **Scaftools** (Briot et al., 2014), **BESST** (Sahlin et al., 2014), **fragScaff** (Adey et al., 2014), in **Chateau and Giroudeau** (2015) and **Weller et al.** (2015), **ScaffMatch** (Mandric and Zelikovsky, 2015), **SALSA2** (Ghurye et al., 2019), **LRScaf** (Qin et al., 2019), in **Davot et al.** (2022) and **Aganezov et al.** (2022). Concerning the read assembly stage, *UG* can be found in **Myers** (2005) and **Li** (2016).
► **Definition 1.13: Tail-head fragments based undirected multigraph**
(UG)

Let $UG = (V, E_F, E_L, \Phi_L)$ be an multiungraph such that:

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

- $E_F = \{\{v_t, v_h\} | v_t, v_h \in V, v_t \neq v_h\}$ is the set of fragment-edges ($v_t$ for the tail, $v_h$ for the head), where each $e \in E_F$ represents both the forward and the reverse orientation of one fragment;

- $E_L$ is the multiset of link-edges, where each $e \in E_L$ represents both one link and its reverse;

- $\Phi_L: E_L \to \{\{u, v\} | u, v \in V, u \neq v\}$ is the incidence function that associates a link-edge with two vertices that are not connected by a fragment-edge.

Denote by:

$$ \text{frag}: \{(u, v) | \{u, v\} \in E_F\} \twoheadrightarrow \mathcal{F} $$

the bijective function that associates each ordered fragment-edge with one and only one (oriented) fragment, and by:

$$ \text{link}: \{(u, v, e) | e \in E_L, \Phi_L(e) = \{u, v\}\} \twoheadrightarrow \mathcal{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 Chateau and Giroudeau (2015), we can define the set $E_F$ of fragment-edges as a perfect matching of the vertices. Figure 20 illustrates each link case in $UG$.

► **Property 1.6: Sizes of UG**

For $UG = (V, E_F, E_L, \Phi_L)$, $|V| = |\mathcal{F}|$, $E_F = \frac{1}{2}|\mathcal{F}|$ and $|E_L| = \frac{1}{2}|\mathcal{L}|$.

Analogously to Definitions 1.8 and 1.11, Definition 1.14 formalises the reverse operation in $UG$.
<figure><img src="image_42.png" /><figcaption>(a) $(a_f, b_f) \in \mathcal{L}'$</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>(b) $(b_f, a_f) \in \mathcal{L}'$</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>(c) $(a_f, b_r) \in \mathcal{L}'$</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>(d) $(b_r, a_f) \in \mathcal{L}'$</figcaption></figure>

■ **Figure 20 – Link cases in UG.**

$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 UG.

► **Definition 1.14: The reverse operation in UG**

— Let $\{u, v\} \in E_F$. Hence $\text{frag}(v, u) = \overline{\text{frag}(u, v)}$.

— Let $e \in E_L, \Phi_L(e) = \{u, v\}$. Hence $\text{link}(v, u, e) = \overline{\text{link}(u, v, e)}$.

Definition 1.15 specifies the conditions for a path in UG to be valid, illustrated in Figure 21.

► **Definition 1.15: Path in UG**

A path $p = (v_0, v_1, \dots, v_{2k-1}) \in V^{2k}, k \in \mathbb{N}$ in UG = $(V, E_F, E_L, \Phi_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:

$$ \forall i \in [0, k], \{v_{2i}, v_{2i+1}\} \in E_F $$

— even-numbered edges in the path are link-edges:

$$ \forall i \in [1, k], \exists e \in E_L \mid \Phi_L = \{v_{2i-1}, v_{2i}\} $$

**Exclusive orientation** $\forall v_i, v_j \in p, v_i = v_j \implies i = j$

Note that walking in UG 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 UG.

<figure><img src="image_46.png" /><figcaption>■ **Figure 21 – A path in UG.**</figcaption></figure>

A valid path is coloured in red. As UG 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 UG, that inherits from Property 1.1.

► **Property 1.7: Reverse symmetry in UG**

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

— $\forall a \in \mathcal{F}, \exists! \{u, v\} \in E_{\mathcal{F}} \text{ s.t.:}$

$$ \text{frag}(u, v) = a \sqcup \text{frag}(v, u) = a $$

— $\forall l \in \mathcal{L}, \exists! e \in E_{\mathcal{L}}, \Phi_{\mathcal{L}}(e) = \{u, v\} \text{ s.t.:}$

$$ \text{link}(u, v, e) = l \sqcup \text{link}(v, u, e) = l $$

▷ Proof

From Definition 1.13, $\text{frag}: \{(u, v) | \{u, v\} \in E_{\mathcal{F}}\} \iff \mathcal{F}$ and $\text{link}: \{(u, v, e) | e \in E_{\mathcal{L}}, \Phi_{\mathcal{L}}(e) = \{u, v\}\} \iff \mathcal{L}$ are bijective.

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

Quick reference

For a quick overview of the different fragment graph structures, the readers can refer to Figure 22.
<figure><img src="image_49.png" /><figcaption>■ **Figure 22 – Overview of all the fragment graph structures.**</figcaption></figure>

Let *a* and *b* two fragments. Each vertical dashed line separates two graph structures. Each colour is associated with a link case (and its reverse): (*a*<sub>*f*</sub>, *b*<sub>*f*</sub>) ∈ 𝒰' in blue, (*b*<sub>*f*</sub>, *a*<sub>*f*</sub>) ∈ 𝒰' in green, (*a*<sub>*f*</sub>, *b*<sub>*r*</sub>) ∈ 𝒰' in yellow and (*b*<sub>*r*</sub>, *a*<sub>*f*</sub>) ∈ 𝒰' 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.

## 2 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 paired-end 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, UG) can be compared in Figure 22.

## 2.1 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 theDBG 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 *UG* respectively stand for directed, bidirected and undirected graphs.

<table><thead><tr><th>Dedicated input data</th><th>Approach</th><th>Software</th><th>Graph</th></tr></thead><tbody><tr><td rowspan="12">Paired-end/Mate-pair</td><td>Huson et al. (2002)</td><td>–</td><td>UG</td></tr><tr><td>Pop et al. (2004)</td><td>Bambus</td><td>BG</td></tr><tr><td>Dayarian et al. (2010)</td><td>SOPRA</td><td>BG</td></tr><tr><td>Koren et al. (2011)</td><td>Bambus 2</td><td>BG</td></tr><tr><td>Salmela et al. (2011)</td><td>MIP</td><td>BG</td></tr><tr><td>Gao et al. (2011)</td><td>Opera</td><td>BG</td></tr><tr><td>Boetzer et al. (2011)</td><td>SSPACE</td><td>BG**</td></tr><tr><td>Roy et al. (2012)*</td><td>SLIQ</td><td>–</td></tr><tr><td>Gritsenko et al. (2012)</td><td>GRASS</td><td>BG</td></tr><tr><td>Donmez and Brudno (2013)</td><td>SCARPA</td><td>UG</td></tr><tr><td>Sahlin et al. (2014)</td><td>BESST</td><td>UG</td></tr><tr><td>Mandric and Zelikovsky (2015)</td><td>ScaffMatch</td><td>UG</td></tr><tr><td>Luo et al. (2017)</td><td>BOSS</td><td>BG</td></tr><tr><td>Andonov et al. (2019)</td><td>GAT</td><td>DG</td></tr><tr><td rowspan="2">Hi-C</td><td>Burton et al. (2013)</td><td>LACHESIS</td><td>BG</td></tr><tr><td>Ghurye et al. (2019)</td><td>SALSA2</td><td>UG</td></tr><tr><td rowspan="2">Long reads</td><td>Qin et al. (2019)</td><td>LRScaf</td><td>UG</td></tr><tr><td>Luo et al. (2019)</td><td>SLR</td><td>BG</td></tr><tr><td rowspan="2">Linked reads</td><td>Adey et al. (2014)</td><td>fragScaff</td><td>UG</td></tr><tr><td>Coombe et al. (2018)*</td><td>ARKS</td><td>–</td></tr><tr><td rowspan="4">Generic link</td><td>Briot et al. (2014)</td><td>Scaftools</td><td>UG</td></tr><tr><td>Chateau and Giroudeau (2015)</td><td>–</td><td>UG</td></tr><tr><td>Weller et al. (2015)</td><td>–</td><td>UG</td></tr><tr><td>Davot et al. (2022)</td><td>–</td><td>UG</td></tr></tbody></table>

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).

<footnote>*SLIQ and ARKS are pre-scaffolder methods.**Suggested in supplementary method.</footnote>
**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.

<table><thead><tr><th>Sequencing data</th><th>Link</th><th>Distance</th></tr></thead><tbody><tr><td>Paired-end*</td><td>✓</td><td>✓</td></tr><tr><td>Mate pair*</td><td>✓</td><td>✓</td></tr><tr><td>Hi-C</td><td>✓</td><td></td></tr><tr><td>Long reads</td><td>✓</td><td>✓</td></tr><tr><td>Reference</td><td>✓</td><td>✓</td></tr><tr><td>Optical mapping</td><td>✓</td><td>✓</td></tr><tr><td>Linked reads**</td><td>✓</td><td></td></tr></tbody></table>

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

<footnote>*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.</footnote>
<figure><img src="image_50.png" /></figure>

■ **Figure 23 – 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) OPM 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$.
### 2.2.1 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 Huson et al. (2002), Dayarian et al. (2010), Koren et al. (2011), Gao et al. (2011), Gritsenko et al. (2012), Salmela et al. (2011), Donmez and Brudno (2013), Briot et al. (2014), Mandric and Zelikovsky (2015) and Luo et al. (2017) adopt with paired-end data, Burton et al. (2013) and Ghurye et al. (2019) for Hi-C data, and Adey et al. (2014) and Coombe et al. (2018) 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.

### 2.2.2 Removing paired-end reads

**Detecting artefact data** Low mapping score are ignored, as Luo et al. (2017) does. If only one read of the pair is aligned, the whole pair is ignored (Burton et al., 2013; Coombe et al., 2018). Mandric and Zelikovsky (2015) and Luo et al. (2017) remove paired-end alignments from the bundle if the distance they provide between the contigs is over a confident interval. Luo et al. (2017) 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, Ghurye et al. (2019) trim the reads that are aligned after a ligation.

**Detecting repeats** When one of the two paired reads maps on distinct contigs, Mandric and Zelikovsky (2015), Burton et al. (2013) and Ghurye et al. (2019) do not consider the pair. Luo et al. (2017) ignore the paired-end reads when they map to a highly covered contig.

### 2.2.3 Removing contigs

**Detecting artefact data** Dayarian et al. (2010) 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, Burton et al. (2013) remove a contig if it does not contain enough restriction sites (only for their clustering step). In the long read scaffolding method of Luo et al. (2019), short contigs are ignored. Finally, Qin et al. (2019) 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 (Gao et al., 2011; Sahlin et al., 2014;
Mandric and Zelikovsky, 2015; Qin et al., 2019). For Hi-C data, Burton et al. (2013) identify a contig as a repeat if the number of restriction sites is high. Luo et al. (2019) remove a contig if it maps in the middle of different long reads, and if its neighbour contigs differ. Koren et al. (2011) computes all the shortest path between two contigs: a contig appearing in too many paths is dismissed.

### 2.2.4 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 (Gao et al., 2011; Gritsenko et al., 2012; Donmez and Brudno, 2013; Luo et al., 2017; Qin et al., 2019). Note that in Gao et al. (2011) 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 (Sahlin et al., 2014; Qin et al., 2019). In Ghurye et al. (2019), the confidence of neighbours links is compared, and only the “best buddy weighted” link is kept (a $O(|\mathcal{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. Qin et al. (2019) dismiss a contig alignment if the long read is contained in the contig.

**Reducing redundant links** As proposed for read assembly OLC approach, Huson et al. (2002), Koren et al. (2011) and Qin et al. (2019) compute a transitive reduction of the links based on their distance distribution.

**Aggressive removing** Roy et al. (2012) 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. Dayarian et al. (2010), Luo et al. (2017) and Qin et al. (2019) keep the maximum-weighted link, in addition, Koren et al. (2011) remove the smallest number of unconfident alignments in the bundle according the distance distribution. Coombe et al. (2018) keep only one link if it passes a binomial test of being the real link among the others. Ghurye et al. (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(|\mathcal{L}|(|\mathcal{C}| + |\mathcal{L}|))$ procedure). Finally, if two contigs are aligned against the same end coordinate of a long read, the best mapping is kept in Luo et al. (2019).
### 2.2.5 Partitioning the instances

Biconnected components can be solved independently For *BG* and *UG* fragment graphs (Sections 1.3 and 1.4), finding articulation edges and removing them enables to partition the graph into biconnected components. Dayarian et al. (2010) show that each component can be solved independently without a loss of optimality. Salmela et al. (2011) and Donmez and Brudno (2013) also apply this divide and conquer strategy. Finding all the articulation points can be done in $O(|C|^2)$

Clustering according to the number of chromosomes In the Burton et al. (2013) Hi-C scaffolding method, the number of chromosomes *N* is an input. Based on the links' weight, the contig set is partitioned in *N* parts with a hierarchical agglomerative clustering approach. Note that the clustering does not take into consideration the relative orientations of the contigs.

## 2.3 Orienting the contigs

Here we present the methods that focus on the contig orientation problem independently of their ordering. Kececioglu and Myers (1995) 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.

### ► Definition 2.1: Orientation-valid layout

Denote by $\mathcal{C}^{\mathcal{C}} = \{c_0, \dots, c_{n-1}\}$ a subset of the contig set $\mathcal{C}$ and by $\mathcal{L}^{\mathcal{C}}$ a subset of the link set $\mathcal{L}$. Without considering an ordering of the contigs in $\mathcal{C}^{\mathcal{C}}, \mathcal{L}^{\mathcal{C}}$ is an orientation-valid layout if and only if there are exactly two orientation assignments $o = (o_0, \dots, o_{n-1}) \in \{f, r\}^{|\mathcal{C}|}$ and $\bar{o} = (\overline{o_0}, \dots, \overline{o_{n-1}})$ and there exists a balanced partition of $\mathcal{L}^{\mathcal{C}} = \mathcal{L}_f^{\mathcal{C}} \cup \mathcal{L}_r^{\mathcal{C}}$ such that:

$$- |\mathcal{L}_f^{\mathcal{C}}| = |\mathcal{L}_r^{\mathcal{C}}| \text{ and } \forall l \in \mathcal{L}_f^{\mathcal{C}}, \bar{l} \in \mathcal{L}_r^{\mathcal{C}};$$

- $o$ and $\bar{o}$ are respectively consistent in $\mathcal{L}_f^{\mathcal{C}}$ and $\mathcal{L}_r^{\mathcal{C}}$, i.e. $\forall l \in \mathcal{L}_f^{\mathcal{C}}$ where $linkp(l) = (a, b)$, and so $\bar{l} \in \mathcal{L}_r^{\mathcal{C}}$ where $linkp(\bar{l}) = (\bar{b}, \bar{a})$, $\exists! i, j \in [0, n[, i \neq j$ such that:

$$- a = (c_i, o_i) \iff \bar{a} = (c_i, \overline{o_i});$$

$$- b = (c_j, o_j) \iff \bar{b} = (c_j, \overline{o_j}).$$
Figure 24 illustrates Definition 2.1 for the bigraph (BG) and for the ungraph (UG) as defined in Sections 1.3 and 1.4, and illustrated in Figure 22.

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

■ **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 $\bar{o}$) such that all the valid paths assign vertex orientations respecting $o$, otherwise $\bar{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.

2.3.1 Maximising the sum of used link weights

As originally proposed by Kececioglu and Myers (1995), 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. Kececioglu (1991) reduced this formulation to the maximum-weight cut
problem whose decision version is $\mathcal{NP}$-complete. Pop et al. (2004) 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. Dayarian et al.
(2010) aim to minimise a cost function (simply, the sign of the weight is inverse
comparing to the other methods). In Bambus 2 (Koren et al., 2011) (as for its
predecessor Bambus), 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.

Luo et al. (2017) 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. Luo et al. (2019) 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

Donmez and Brudno (2013) 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 *UG*. 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 *NP*-hard problem, Lokshtanov et al. (2009) propose a fixed-parameter tractable algorithm in $O(3^k|E||V|)$, where *k* is the number of vertices to suppress, and *V* and *E* are respectively the vertex and the edge sets.

### 2.3.3 Maximum log-likelihood function

At the opposite of the previous methods, Burton et al. (2013) 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.

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

Pop et al. (2004) reduce the ordering problem to the optimal linear arrangement problem. It consists in assigning a coordinate on an axe to each contig, a *NP*-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 Luo et al. (2019) for long read based scaffolding.

Donmez and Brudno (2013) refer to the feedback arc set problem (*NP*-hard, Karp 1972) and define the contig ordering problem as finding a minimal set of
edges whose removal leads to a directed acyclic graph. Koren et al. (2011) position
the contigs in confidence intervals of three times the standard deviation.

2.4.2 Minimum spanning tree

For Hi-C data, Burton et al. (2013) 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.

2.4.3 Contig positioning only

Donmez and Brudno (2013) 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.

Gao et al. (2011) 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 (Goldfarb and Idnani,
1983).

2.5 Orienting and ordering the contigs simultaneously

Already for the read assembly stage, Kececioglu (1991) 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

Huson et al. (2002) aim to maximise the sum of weights of happy mate-edges in
*UG*. 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 *NPN*-complete
by reducing the bandwidth problem to it.
Gao et al. (2011) focus on a bounded version of the graph bandwidth problem, for which the complexity is in $O(|\mathcal{C}|^w|\mathcal{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, Salmela et al. (2011) 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.

Gritsenko et al. (2012) 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

Briot et al. (2014) maximise the weight of the remaining links such that at most one link connects each contig's end in *UG*. This approach can result in several paths. Building on this, Weller et al. (2015) are looking for at most $\sigma_p$ paths and $\sigma_c$ cycles (where $\sigma_p, \sigma_c$ are parameters) that cover all the contig-edges in *UG*. They follow the scaffolding definition given in Chateau and Giroudeau (2014), 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 $\sigma_p$ paths and $\sigma_c$ cycles. Except in the case where $n = \sigma_p + 2\sigma_c$, for which the problem is in $\mathcal{P}$, the scaffolding problem (decision version) remains $\mathcal{NP}$-complete.

In the global Integer Linear Programming (ILP) formulation of Andonov et al. (2019), 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

Mandric and Zelikovsky (2015) reduce the scaffolding problem to the Maximum-Weight Acyclic 2-Matching (MWA2M) problem on the link-edges in *UG*. The latter one's decision version is $\mathcal{NP}$-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(|\mathcal{L}||\mathcal{C}|^3)$.

Concerning the scaffolding with linked read, Adey et al. (2014) adopt a similar strategy as in Burton et al. (2013) 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, Qin et al. (2019) only extract valid simple paths in *UG* and return them as scaffolds.

## 2.6 Solving approaches

### 2.6.1 Greedy approaches

As in the majority of the cases the scaffolding problem (orienting and ordering separately or in the same time) is *NP*-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. Huson et al. (2002) are the first to present the greedy path merging algorithm. Koren et al. (2011) adopt this strategy for the orienting and the ordering of the contigs separately. They use the same algorithm as in Kececioglu and Myers (1995): a $O(|\mathcal{C}| + |\mathcal{L}|)$ complexity algorithm on *BG*. Mandric and Zelikovsky (2015) propose a greedy heuristic for the MWA2M problem that chooses the heaviest feasible link-edge in *UG*. 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(|\mathcal{C}|\log(|\mathcal{C}|))$ with a max heap implementation. In the continuity of Chateau and Giroudeau (2015), Davot et al. (2022) propose a $O(|\mathcal{C}||\mathcal{L}|\sigma_c^2)$ 3-approximation algorithm for connected cluster graphs (*UG*), where $\sigma_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. Boetzer et al. (2011) 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 Pop et al. (2004), Luo et al. (2017) first order the longest contigs
and connect them with a linear BFS procedure.

**Simulated annealing** Dayarian et al. (2010) 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.

**Spanning tree** As described in the previous sections, Roy et al. (2012) implement a spanning tree algorithm for the orienting of the contigs while Adey et al. (2014) implement it for both the orienting and the ordering of the contigs.

**Majority voting on predictions** For the positioning of the contigs, Roy et al. (2012) employ a majority voting on the ordering predictions based on inequalities.

**Randomised Greedy Algorithms for the maximum matching problem** Finally, Ghurye et al. (2019) benefiting from the algorithm proposed in Poloczek and Szegedy (2012) running in $O(|\mathcal{L}|)$ with UG. The total complexity of their approach raises to $O(|\mathcal{C}|\cdot|\mathcal{L}|)$.

## 2.6.2 Fix parameters dynamic programming

Gao et al. (2011) propose a $O(|\mathcal{C}|^w|\mathcal{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. Weller et al. (2015) give an exact algorithm in $O(|\mathcal{C}|\sigma_p^2)$ for trees (UG graph), where $\sigma_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}\sigma_p\sigma_c|\mathcal{C}|)$, where $tw$ is the tree width for the decomposition, and $\sigma_c$ is equivalent to $\sigma_p$ for cycles.

## 2.6.3 Mathematical programming

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

**Mixed Integer Linear Programming (MILP)** Salmela et al. (2011) 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. Briot et al. (2014) however model the link choice variables as binaries and connect a vertex with at most one link in *UG*. Luo et al. (2017) maximise the sum of weights of used links for orienting the contigs on iteratively increasing subgraphs. Andonov et al. (2019) model in one objective function the longest path and the maximum number of satisfied constraints on distances. Luo et al. (2019) apply the same strategy as in Luo et al. (2017) on both the orienting and the ordering problems.

**Mixed Integer Quadratic Programming (MIQP)** Gritsenko et al. (2012) 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)** Donmez and Brudno (2013) use linear programming to position the oriented and ordered contigs to minimise the positioning approximation according to the distance distributions on the links.

# 3 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, Freudenthal et al. (2020) 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 deal with the structural haplotypes 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).

## 3.1 Chloroplast sequence extraction

Although there are routine data filters, as in Bakker et al. (2016) and McKain and afinit (2017), here we discuss specific approaches to identifying chloroplast data (original and transformed) and extracting them.
■ **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. **Structural haplotypes:** if the method is chloroplast structural haplotypes aware, i.e. if it addresses the issue of several chloroplast forms in the sequencing data.

<table><thead><tr><th>Paper and software</th><th>Approach</th><th>Filter</th><th>Structural haplotypes</th></tr></thead><tbody><tr><td>Coissac et al. (2016)*<br>ORG.Asm</td><td>Seed-and-Extend</td><td>?</td><td>?</td></tr><tr><td>Bakker et al. (2016)<br>IOGA</td><td>Pipeline</td><td>••</td><td></td></tr><tr><td>Dierckxsens et al. (2017)<br>NOVOPlasty</td><td>Seed-and-Extend</td><td>•</td><td></td></tr><tr><td>McKain and afinit (2017)<br>Fast-Plast</td><td>Pipeline</td><td>••</td><td></td></tr><tr><td>Ankenbrand et al. (2018)<br>chloroExtractor</td><td>Pipeline</td><td>••</td><td></td></tr><tr><td>Sancho et al. (2018)<br>Chloroplast assembly protocol</td><td>Pipeline</td><td>•</td><td></td></tr><tr><td>Andonov et al. (2019)<br>GAT</td><td>Scaffolding</td><td></td><td>✓</td></tr><tr><td>Jin et al. (2020)<br>GetOrganelle</td><td>Pipeline<br>+ Scaffolding</td><td>•••</td><td>✓</td></tr></tbody></table>

### 3.1.1 Filtering the reads

**De novo filtering** Before any assembly of the raw data, Ankenbrand et al. (2018) 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 Sancho et al. (2018) align the reads against one reference chloroplast genome, Ankenbrand et al. (2018) against several closely related species, Bakker et al. (2016) against several but not

<footnote>*Cannot find a description of the method, software manual not found.</footnote>
necessarily closely related. Jin et al. (2020) in addition align the reads against
organelle subsequences. In the latter, they denote the reads they extract as seed
reads. Note that in Sancho et al. (2018) 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. Bakker et al. (2016) 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. Jin et al. (2020) extend the seed contigs (from the assembly of the seed reads) with the remaining reads.

3.1.2 Filtering the contigs

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

*De novo filtering* McKain and afinit (2017) and Jin et al. (2020) 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, Jin et al. (2020) 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.

3.2 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.

3.2.1 De Bruijn graph approach

Recall that **DBG** is a **dbgraph** read assembly approach.

**Baiting assemblies** As mentioned before, Bakker et al. (2016) iteratively produce several read assemblies with several sizes of *k*-mer (recall Chapter I, Section 5.3) thanks to SOAPdenovo2 (Luo et al., 2012). 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** Sancho et al. (2018) use Velvet (Zerbino and Birney, 2008), whereas Bakker et al. (2016), Ankenbrand et al. (2018) and Jin et al. (2020) prefer SPAdes that internally tests different size of *k*-mer (Bankevich et al., 2012). Unlike 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.

### 3.2.2 Seed-and-extend

The S&E approach is a local assembling approach adapted for “simple” genomes, particularly with few repeats. Coissac et al. (2016) implement a S&E whose description could not be found at the time of writing this thesis. The S&E strategy in Dierckxsens et al. (2017) 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.

## 3.3 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 Sancho et al. (2018) use the paired-end based scaffolder SSPACE (Boetzer et al., 2011). Finally, the latter fill the gap in the scaffolds with GapFiller (Boetzer and Pirovano, 2012). Andonov et al. (2019) 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** Jin et al. (2020) 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.
## 3.4 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 corresponding to structural haplotypes.

**De novo evaluation** Bakker et al. (2016) select the best assemblies according to the Assembly Likelihood Estimation (ALE) (Clark et al., 2013). The assembly candidates come from both the final assemblies during the iterative read baiting step and the final step.

**Reference genome evaluation** Ankenbrand et al. (2018) keep the contigs that align sufficiently against reference genomes.

**Multiple solutions** Andonov et al. (2019) detect in the *DG* path, representing the scaffold, pairs of subpaths $p = (v_0, \dots, v_{n-1})$, $p' = (v'_0, \dots, v'_{n-1})$ such that: $\forall i \in [0, n[, frag(v_i) = \overline{frag(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 structural haplotypes. 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 Andonov et al. (2019).
<figure><img src="image_52.png" /></figure>
<figure><img src="image_53.png" /></figure>

# III FRAGMENT GRAPH IMPLEMENTATIONS AND COMPARISON

David Krakauer & Kathleen Tagg. (2020). The Geyser [Song]. On *Breath & Hammer*. Table Pounding Records

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

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

## In this chapter

<table><tbody><tr><td>1</td><td>Implementations</td><td>58</td></tr><tr><td>1.1</td><td>Directed graph (DG): oriented fragments based</td><td>59</td></tr><tr><td>1.1.1</td><td>All oriented fragments directed graph (DGA)</td><td>59</td></tr><tr><td>1.1.2</td><td>Oriented fragments' successors directed graph (DGS)</td><td>61</td></tr><tr><td>1.1.3</td><td>Forward fragments directed graph (DGF)</td><td>63</td></tr><tr><td>1.2</td><td>Bidirected graph (BG): oriented walk based</td><td>66</td></tr><tr><td>1.2.1</td><td>Unoriented fragments bidirected graph (BGU)</td><td>66</td></tr><tr><td>1.2.2</td><td>Transformation to DGF</td><td>68</td></tr><tr><td>1.3</td><td>Undirected graph (UG): tail-head fragments based</td><td>68</td></tr><tr><td>1.3.1</td><td>All oriented fragments undirected graph (UGA)</td><td>68</td></tr><tr><td>1.3.2</td><td>Transformation to DGS</td><td>70</td></tr><tr><td>1.4</td><td>Fragment graph map</td><td>70</td></tr><tr><td>2</td><td>Algorithms for DGS, DGF and BGU</td><td>72</td></tr><tr><td>2.1</td><td>Subfunctions</td><td>72</td></tr><tr><td>2.2</td><td>Iterating over the predecessors</td><td>73</td></tr><tr><td>2.3</td><td>Iterating over the successors</td><td>74</td></tr><tr><td>2.4</td><td>Adding a vertex</td><td>75</td></tr><tr><td>2.5</td><td>Adding an edge</td><td>76</td></tr><tr><td>2.6</td><td>Deleting a vertex</td><td>78</td></tr><tr><td>2.7</td><td>Deleting an edge</td><td>81</td></tr><tr><td>3</td><td>Time costs</td><td>82</td></tr><tr><td>4</td><td>Memory and time cost comparisons</td><td>84</td></tr></tbody></table>
# 5 Conclusions and perspectives

85

Section 1 describes several graph implementations for each structure from Chapter II, Section 1 (see the fragment graph structures overview in Figure 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.

# 1 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 $\mathcal{F}$ and to the link set $\mathcal{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 ($\mathcal{F} = \mathcal{F}_f \sqcup \mathcal{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, $\mathcal{L} = \mathcal{L}_f \sqcup \mathcal{L}_r$. $\mathcal{L}_f$ contains the links for which in the pair, the fragment with the smallest $\Sigma_{\mathcal{F}}$ label is in forward orientation, i.e. $\mathcal{L}_f = \{l \in \mathcal{L} \mid \text{for}(a) = f, \text{fid}(a) = \min_{x \in \text{link}_p(l)} \text{fid}(x)\}$. Similarly to the unoriented fragment label set $\Sigma_{\mathcal{F}_{\text{un}}}$, let $\Sigma_{\mathcal{L}_{\text{can}}} = [0, \frac{1}{2}|\mathcal{L}|[}$ be the label set of canonical links and *canlid*: $\mathcal{L} \twoheadrightarrow \Sigma_{\mathcal{L}_{\text{can}}}$ its associated function, such that $\forall l \in \mathcal{L}, \text{canlid}(l) = \text{canlid}(\bar{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 $\mathcal{L}_f$ or in $\mathcal{L}_r$.

**Adjacency list notation** The adjacency lists are represented by functions $f: \mathbb{N} \rightarrow 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 \in 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 \in \mathbb{N}$, $oct(n) = \lfloor \frac{\log_2 n}{8} \rfloor$ gives the number of octets necessary to memorise an integer between 0 and $n-1$.

## 1.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*, *Φ*).

### 1.1.1 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.
<figure><img src="image_56.png" /><figcaption>LEGEND</figcaption></figure>

■ **Figure 25 – DGA 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 \in V$, $vind(v) = v_f$ if $frag(v) \in \mathcal{F}_f$, $v_r$ otherwise. For each link and its reverse, $e$ is their canonical label.

► **Definition 1.1: DGA labels and functions**

Let $vind: V \hookrightarrow \Sigma_V$ and $eind: E \hookrightarrow \Sigma_E$, where $\Sigma_V = \Sigma_F$ and $\Sigma_E = \Sigma_{L_{can}}$ are the label sets for the vertices and the edges, respectively. The predecessor and the successor lists are defined by two functions:

$$ \text{predl}: \Sigma_V \to N_V^- = \{N_v^-, \forall v \in V \mid N_v^-: \llbracket 0, d_v^- \rrbracket \to \Sigma_V \times \Sigma_E \} $$

$$ \text{succl}: \Sigma_V \to N_V^+ = \{N_v^+, \forall v \in V \mid N_v^+: \llbracket 0, d_v^+ \rrbracket \to \Sigma_V \times \Sigma_E \} $$

For each vertex $v \in V$, $N_v^-: \llbracket 0, d_v^- \rrbracket \to \Sigma_V \times \Sigma_E$ and $N_v^+: \llbracket 0, d_v^+ \rrbracket \to \Sigma_V \times \Sigma_E$ give its predecessors and its successors, respectively. The following property
holds:

$$
\begin{align*}
\forall e \in E, \Phi(e) &= (u, v), \exists! n \in [0, d_u^+], \exists! m \in [0, d_v^-] \text{ s.t.} \\
\text{succl}(vind(u))(n) &= (vind(v), eind(e)) \\
\text{predl}(vind(v))(m) &= (vind(u), eind(e))
\end{align*}
$$

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 \in \mathcal{F}_f$ and $a_r = \bar{a}_f \in \mathcal{F}_r$ and their associated vertices $v_f$ and $v_r$ ($frag(v_f) = a_f$, $frag(v_r) = a_r$), the labels of $v_f$ and $v_r$ respect $vind(v_r) - vind(v_f) = 1$. For each edge $e \in E$, its label respects $eind(e) = eind(\bar{e})$.

Proposition 1.1 gives the amount of octet *DGA* consumes.

► **Proposition 1.1: DGA memory consumption**

The memory size *Mem(DGA)* of *DGA* (in octets) is equal to:

$$
Mem(DGA) = 2P(|\mathcal{F}| + 1) + 2|\mathcal{L}| \left( \text{oct}(|\mathcal{F}|) + \text{oct}\left(\frac{1}{2}|\mathcal{L}|\right) \right)
$$

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|)).

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

1.1.2 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.
<figure><img src="image_58.png" /><figcaption>■ Figure 26 – DGS implementation.</figcaption></figure>

Each squared mathematical formula provides the link case and its reverse (in parentheses).
The list below 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 \hookrightarrow \Sigma_V$ and $eind: E \hookrightarrow \Sigma_E$, where $\Sigma_V = \Sigma_F$ and $\Sigma_E = \Sigma_{L_{can}}$ are the label sets for the vertices and the edges, respectively. The successor list is defined by:

$$
\mathrm{succl}: \Sigma_V \to N_V^+ = \{N_v^+, \forall v \in V \mid N_v^+ : [0, d_v^+ [ \to \Sigma_V \times \Sigma_E]\}
$$

For each vertex $v \in V$, $N_v^- : [0, d_v^-] \to \Sigma_V \times \Sigma_E$ and $N_v^+ : [0, d_v^+] \to \Sigma_V \times \Sigma_E$
give its predecessors and its successors, respectively. The following property
holds:

$$ \forall e \in E, \Phi(e) = (u, v), \exists! n \in [0, d_u^+], \exists! m \in [0, d_{\bar{v}}^+] \text{ s.t.} $$

$$ \mathrm{succl}(\mathrm{vind}(u))(n) = (\mathrm{vind}(v), \mathrm{eind}(e)) $$

$$ \mathrm{succl}(\mathrm{vind}(\bar{v}))(m) = (\mathrm{vind}(\bar{u}), \mathrm{eind}(\bar{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 \in \mathcal{F}_f$ and $a_r = \overline{a_f} \in \mathcal{F}_r$ and their associated vertices $v_f$ and $v_r$ ($\mathrm{frag}(v_f) = a_f, \mathrm{frag}(v_r) = a_r$), the labels of $v_f$ and $v_r$ respect $\mathrm{vind}(v_r) - \mathrm{vind}(v_f) = 1$. For each edge $e \in E$, its label respects $\mathrm{eind}(e) = \mathrm{eind}(\bar{e})$.

Proposition 1.2 gives the amount of octet DGS consumes.

### ► Proposition 1.2: DGS memory consumption

The memory size $Mem(DGS)$ of DGS (in octets) is equal to:

$$ Mem(DGS) = P(|\mathcal{F}| + 1) + |\mathcal{L}| \left( \mathrm{oct}(|\mathcal{F}|) + \mathrm{oct}\left(\frac{1}{2}|\mathcal{L}|\right) \right) $$

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(|\mathcal{F}| + 1)$. Then, for each link, for each successor, the neighbour's label and the edge's label are provided, i.e. $|\mathcal{L}|(\mathrm{oct}(|\mathcal{F}|) + \mathrm{oct}(\frac{1}{2}|\mathcal{L}|))$. △

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

### 1.1.3 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.
<figure><img src="image_60.png" /><figcaption>■ **Figure 27 – DGF implementation.**</figcaption></figure>

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 \to \Sigma_{V_f}$ and $eind: E \to \Sigma_E$, where $\Sigma_{V_f} = \Sigma_{\mathcal{F}_{un}}$ and $\Sigma_E = \Sigma_{\mathcal{L}_{can}}$
are the label sets for the forward vertices and the edges, respectively. The
predecessor and the successor lists are defined by two functions:

$$
\texttt{predl: } \Sigma_V \to N_{V_f}^- = \{N_v^-, \forall v \in V_f \mid N_v^- : \llbracket 0, d_v^- \rrbracket \to \Sigma_{V_f} \times \{f, r\} \times \Sigma_E \}
$$

$$
\texttt{succl: } \Sigma_V \to N_{V_f}^+ = \{N_v^+, \forall v \in V_f \mid N_v^+ : \llbracket 0, d_v^+ \rrbracket \to \Sigma_{V_f} \times \{f, r\} \times \Sigma_E \}
$$

where $V_f = \{v \in V \mid \text{frag}(v) \in \mathcal{F}_f\}$. For each forward vertex $v \in V_f$, $N_v^- : \llbracket 0, d_v^- \rrbracket \to \Sigma_{V_f} \times \{f, r\} \times \Sigma_E$ and $N_v^+ : \llbracket 0, d_v^+ \rrbracket \to \Sigma_{V_f} \times \{f, r\} \times \Sigma_E$ give its predecessors
and its successors, respectively. The following property holds:

$$
\begin{align*}
& \forall e \in E, \Phi(e) = (u, v), \\
& \mathit{frag}(u) \in \mathcal{F}_f \implies \exists! n \in [0, d_u^+ [ \text{ s.t.} \\
& \phantom{\mathit{frag}(u) \in \mathcal{F}_f \implies \exists! n \in [0, d_u^+ [ \text{ s.t.}} \mathit{succl}(\mathit{vind}(u))(n) = (\mathit{vind}(v), \mathit{for}(\mathit{frag}(v)), \mathit{eind}(e)) \\
& \mathit{frag}(u) \in \mathcal{F}_r \implies \exists! n \in [0, d_u^- [ \text{ s.t.} \\
& \phantom{\mathit{frag}(u) \in \mathcal{F}_r \implies \exists! n \in [0, d_u^- [ \text{ s.t.}} \mathit{predl}(\mathit{vind}(u))(n) = (\mathit{vind}(v), \mathit{for}(\mathit{frag}(\bar{v})), \mathit{eind}(\bar{e})) \\
& \mathit{frag}(v) \in \mathcal{F}_f \implies \exists! m \in [0, d_v^- [ \text{ s.t.} \\
& \phantom{\mathit{frag}(v) \in \mathcal{F}_f \implies \exists! m \in [0, d_v^- [ \text{ s.t.}} \mathit{predl}(\mathit{vind}(v))(m) = (\mathit{vind}(u), \mathit{for}(\mathit{frag}(u)), \mathit{eind}(e)) \\
& \mathit{frag}(v) \in \mathcal{F}_r \implies \exists! m \in [0, d_v^- [ \text{ s.t.} \\
& \phantom{\mathit{frag}(v) \in \mathcal{F}_r \implies \exists! m \in [0, d_v^- [ \text{ s.t.}} \mathit{succl}(\mathit{vind}(v))(m) = (\mathit{vind}(u), \mathit{for}(\mathit{frag}(\bar{u})), \mathit{eind}(\bar{e}))
\end{align*}
$$

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 \in V$, its label respects $vind(v) = vind(\bar{v})$. For each edge $e \in E$, its label respects $eind(e) = eind(\bar{e})$.

Proposition 1.3 gives the amount of octet DGF consumes.

► **Proposition 1.3: DGF memory consumption**

The memory size *Mem(DGF)* of *DGF* (in octets) is equal to:

$$
Mem(DGF) = P(|\mathcal{F}| + 2) + |\mathcal{L}| \left( \mathrm{oct}(|\mathcal{F}|) + \mathrm{oct}\left(\frac{1}{2}|\mathcal{L}|\right) \right)
$$

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).

### 1.2.1 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 ↛ Σ<sub>V</sub> and eind: E ↛ Σ<sub>E</sub>, where Σ<sub>V</sub> = Σ<sub>F<sub>un</sub></sub> and Σ<sub>E</sub> = Σ<sub>L<sub>can</sub></sub> are respectively the label sets for the vertices and the edges. The neighbour and the edge attribute lists are defined by two functions:

$$nborl: \Sigma_V \to N_V = \{N_v, \forall v \in V \mid N_v: \llbracket 0, d_v \rrbracket \to \Sigma_V \times \Sigma_E\}$$

$$attrel: \Sigma_E \to \{f, r\}^2$$

For each vertex $v \in V$, $N_v: \llbracket 0, d_v \rrbracket \to \Sigma_V \times \Sigma_E$ gives its neighbours. The following property holds:

$$\begin{align*}
\forall e \in E, \Phi(e) &= \{u, v\}, \exists! n \in \llbracket 0, d_u \rrbracket, \exists! m \in \llbracket 0, d_v \rrbracket \text{ s.t.} \\
&nborl(vind(u))(n) = (vind(v), eind(e)) \\
&nborl(vind(v))(m) = (vind(u), eind(e)) \\
&attrel(eind(e)) = attre(e)
\end{align*}$$

Proposition 1.4 gives the amount of octet BGU consumes.

#### ► Proposition 1.4: BGU memory consumption

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

$$Mem(BGU) = P\left(\frac{1}{2}|\mathcal{F}| + 1\right) + |\mathcal{L}|\left(oct\left(\frac{1}{2}|\mathcal{F}|\right) + oct\left(\frac{1}{2}|\mathcal{L}|\right)\right) + oct(|\mathcal{L}|)$$
LEGEND

$$l \in \mathcal{L}' \iff (\bar{l} \in \mathcal{L}')$$

$N_v$

Eattr

$$ \begin{array}{ccc} u & \rightarrow & v \end{array} $$

e

uor vor

$$u_f \rightarrow v_f \iff (v_r \rightarrow u_r)$$

$$ \begin{array}{ccc} u & \rightarrow & v \ e \\ v & \rightarrow & u \ e \end{array} $$

e

f

f

$$u_f \rightarrow v_r \iff (v_f \rightarrow u_r)$$

$$ \begin{array}{ccc} u & \rightarrow & v \ e \\ v & \rightarrow & u \ e \end{array} $$

e

f

r

$$u_f \leftarrow v_f \iff (v_r \leftarrow u_r)$$

$$ \begin{array}{ccc} u & \rightarrow & v \ e \\ v & \rightarrow & u \ e \end{array} $$

e

r

r

$$u_f \leftarrow v_r \iff (v_f \leftarrow u_r)$$

$$ \begin{array}{ccc} u & \rightarrow & v \ e \\ v & \rightarrow & u \ e \end{array} $$

e

r

f

■ **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(\frac{1}{2}|\mathcal{F}|+1)$. Each vertex is associated with its neighbours and the edge's label. Each edge is stored twice ($|\mathcal{L}|$). So it gives $|\mathcal{L}|(\text{oct}(\frac{1}{2}|\mathcal{F}|) + \text{oct}(\frac{1}{2}|\mathcal{L}|))$. Finally, *Eattr* provides for each edge two binary values (one bit per orientation), i.e. $\text{oct}(|\mathcal{L}|)$.

<figure><img src="image_61.png" /></figure>
### 1.2.2 Transformation to DGF

It is possible to transform BGU to DGF. For each $u \in \Sigma_V$, for each $(v, e) \in \Sigma_V \times \Sigma_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^+$.

## 1.3 Undirected graph (UG): tail-head fragments based

Section 1.3.1 describes the first implementation for UG and Section 1.3.2 describes how to transform UGA to DGS. Recall from Definition 1.13 that UG = $(V, E_F, E_L, \Phi_L)$.

### 1.3.1 All oriented fragments undirected graph (UGA)

UGA 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 UGA.

**► Definition 1.8: UGA labels and functions**

Let $vind: V \hookrightarrow \Sigma_V$, and $eind: E_L \hookrightarrow \Sigma_E^a$, where $\Sigma_V = \Sigma_F$ and $\Sigma_E = \Sigma_{Lcan}$ are the label sets for the vertices and the link-edges, respectively. The fragment-neighbour function is denoted by $nborlf: \Sigma_V \hookrightarrow \Sigma_V$. For each vertex $v \in V$, $nborlf(v)$ gives the other extremity of the fragment. The link-neighbour list is defined by:

$$nborll: \Sigma_V \rightarrow N_V^L = \{N_v^L, \forall v \in V \mid N_v^L: [0, d_v^L] \rightarrow \Sigma_V \times \Sigma_E\}$$

For each vertex $v \in V$, $N_v^L: [0, d_v^L] \rightarrow \Sigma_V \times \Sigma_E$ provides its link-neighbours. The following property holds:

$$\begin{align*}
& \forall e \in E_L, \Phi_L(e) = \{u, v\}, \exists! n \in [0, d_u^L], \exists! m \in [0, d_v^L \text{ s.t.} \\
& nborll(vind(u))(n) = (vind(v), eind(e)) \\
& nborll(vind(v))(m) = (vind(u), eind(e))
\end{align*}$$
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 \in \mathcal{F}_f$ and its reverse $a_r = \overline{a_f} \in \mathcal{F}_r$ and their associated vertices $\{v_t, v_h\} \in E_{\mathcal{F}}$ ($frag(v_t, v_h) = a_f$, $frag(v_h, v_t) = a_r$), the label of $v_t$

<figure><img src="image_62.png" /><figcaption>■ **Figure 29 – UGA implementation.**</figcaption></figure>

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_V^\mathcal{L}$ 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.

<footnote>ᵃWe do not provide edge labels for fragment-edges.</footnote>
and $v_h$ respect $vind(v_h) - vind(v_t) = 1$.

Proposition 1.5 gives the amount of octet UGA consumes.

### Proposition 1.5: UGA memory consumption

The memory size $Mem(UGA)$ of UGA (in octets) is equal to:

$$Mem(UGA) = P(|\mathcal{F}| + 1) + |\mathcal{L}| \left( \text{oct}(|\mathcal{F}|) + \text{oct}\left(\frac{1}{2}|\mathcal{L}|\right) \right)$$

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(|\mathcal{F}| + 1)$. Then, for each link, for each extremity, the neighbour's label and the edge's label are provided, i.e. $|\mathcal{L}|(\text{oct}(|\mathcal{F}|) + \text{oct}(\frac{1}{2}|\mathcal{L}|))$.

□

## 1.3.2 Transformation to DGS

UGA can easily be adapted to obtain DGS. For each $v \in \Sigma_V$, for each $n \in [0, d_v^\mathcal{L}]$ 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$.

## 1.4 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 theoretical 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).
<figure><img src="image_63.png" /></figure>

LEGEND

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

Graph structures
ordered chronologically

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

Graph implementations
linked by transformations

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

Oriented fragment set

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

Link set

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

Vertex set

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

Edge set

Graph structure

Theoretical graph type

Theoretical number of vertices and edges

Graph implementation

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

Sets of vertices for all the
oriented fragments
and for only the forward

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

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

Adjacency lists

i integer of range $0 \le i < \frac{1}{2}|\mathcal{F}|$

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

Lists of predecessors
and successors

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

ii integer of range $0 \le ii < |\mathcal{F}|$

e integer of range $0 \le e < \frac{1}{2}|\mathcal{L}|$

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

Lists of link-edges
and fragment-edges

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

b boolean value $b \in \{0, 1\}$
## 2 Algorithms for DGS, DGF and BGU

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.

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

**► Algorithm 1: Reverse operation on index**

**Require:** Integer *ind*.

**Ensure:** Returns $2^{\lfloor \frac{ind}{2} \rfloor} + 1 - (ind \bmod 2)$.

1: **function** REV(*ind*)

2: **if** *ind* | 2 **then**

3: **return** *ind* + 1

4: **return** *ind* - 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*)

2: *ind* ← 0

3: **while** list[*ind*][*t*] ≠ *e* and *ind* < |*list*| - 1 **do**

4: *ind* ← *ind* + 1

5: **return** *ind*
### ► Algorithm 3: Delete an entry in a list and arrange it

**Require:** 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*)

2: if *ind* = |*list*| - 1 **then** ► *Just delete the last element*

3: *list*.POP()

4: **else** ► *Replace the element of interest by the last one*

5: *list[*ind*] ← *list*.POP()

## 2.2 Iterating over the predecessors

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

### ► Algorithm 4: Iterate over the predecessors for DGS

**Require:** A vertex label *v* ∈ Σ<sub>V</sub>.

**Ensure:** Returns a list containing tuples in Σ<sub>V</sub> × Σ<sub>E</sub> representing the predecessors of the vertex labelled by *v*.

1: **function** DGS_PREDS(*v*)

2: *preds* ← EMPTY_LIST()

3: **for** (*u*, *e*) ∈ N<sub>v</sub><sup>+</sup> **do**

4: *preds*.APPEND(ū, *e*)

5: **return** *preds*

### ► Algorithm 5: Iterate over the predecessors for DGF

**Require:** A vertex label *v* ∈ Σ<sub>V</sub> associated with an orientation *vor* ∈ {*f*, *r*}.

**Ensure:** Returns a list containing tuples in Σ<sub>V</sub> × {*f*, *r*} × Σ<sub>E</sub> representing the predecessors of the vertex labelled by *v* with fragment orientation *vor*.

1: **function** DGF_PREDS(*v*, *vor*)

2: *preds* ← EMPTY_LIST()

3: if *vor* = *f* **then** ► *Forward fragment*

4: **for** (*u*, *uor*, *e*) ∈ N<sub>v</sub><sup>-</sup> **do**

5: *preds*.APPEND(*u*, *uor*, *e*)

6: **else** ► *Reverse fragment*

7: **for** (*u*, *uor*, *e*) ∈ N<sub>v</sub><sup>+</sup> **do**
8: `preds.append(u, 1 - uor, e)`

9: `return preds`

**► Algorithm 6:** **Iterate over the predecessors for BGU**

**Require:** A vertex label $v \in \Sigma_V$ associated with an orientation $vor \in \{f, r\}$ representing the predecessors of the vertex labelled by $v$ and with fragment orientation $vor$.

**Ensure:** Returns a list containing tuples in $\Sigma_V \times \{f, r\} \times \Sigma_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) \in N_v$ **do**

4: **if** $u < v \land Eattr[e][1] = vor$ **then**

5: `preds.append(u, Eattr[e][0], e)`

6: **else if** $1 - Eattr[e][0] = vor$ **then**

7: `preds.append(u, 1 - Eattr[e][1], e)`

8: `return preds`

## 2.3 Iterating over the successors

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

**► Algorithm 7: Iterate over the successors for DGS**

**Require:** A vertex label $v \in \Sigma_V$.

**Ensure:** Returns a list containing tuples in $\Sigma_V \times \Sigma_E$ representing the successors of the vertex labelled by $v$.

1: **function** DGS succes($v$)

2: `succs = EMPTY_LIST()`

3: **for** $(w, e) \in N_v^+$ **do**

4: `succs.append(w, e)`

5: `return succs`
### ► Algorithm 8: Iterate over the successors for DGF

**Require:** A vertex label $v \in \Sigma_V$ associated with an orientation $vor \in \{f, r\}$.

**Ensure:** Returns a list containing tuples in $\Sigma_V \times \{f, r\} \times \Sigma_E$ representing the predecessors of the vertex labelled by $v$ with fragment orientation $vor$.

1: **function** DGF succes($v, vor$)

2: $succs \leftarrow EMPTY_LIST(\,)$

3: **if** $vor = f$ **then** ▷ Forward fragment

4: **for** ($w, wor, e) \in N_v^+$ **do**

5: $succs$.APPEND($w, wor, e$)

6: **else** ▷ Reverse fragment

7: **for** ($w, wor, e) \in N_v^-$ **do**

8: $succs$.APPEND($w, 1 - wor, e$)

9: **return** $succs$

### ► Algorithm 9: Iterate over the successors for BGU

**Require:** A vertex label $v \in \Sigma_V$ associated with an orientation $vor \in \{f, r\}$ representing the successors of the vertex labelled by $v$ and with fragment orientation $vor$.

**Ensure:** Returns a list containing tuples in $\Sigma_V \times \{f, r\} \times \Sigma_E$ representing the predecessors of the vertex labelled by $v$ with the fragment orientation $vor$.

1: **function** BGU succes($v, vor$)

2: $succs \leftarrow EMPTY_LIST(\,)$

3: **for** ($w, e) \in N_v$ **do**

4: **if** $v < w \land 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$)

8: **return** $succs$

## 2.4 Adding a vertex

The following algorithms allocate a place for new vertices by providing their label.
► Algorithm 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()

2: $N_V^+$.APPEND(EMPTY_LIST())

3: $N_V^+$.APPEND(EMPTY_LIST())

4: **return** $|N_V^+| - 2$

► Algorithm 11: Add a new vertex in the graph for DGF

**Ensure:** The lengths of $N_V^-$ and $N_V^+$ are increased by one and returns the last vertex label.

1: **function** DGF_ADD VERTEX()

2: $N_V^-$.APPEND(EMPTY_LIST())

3: $N_V^-$.APPEND(EMPTY_LIST())

4: **return** $|N_V^-| - 1$

► Algorithm 12: Add a new vertex in the graph for BGU

**Ensure:** The length of $N_V$ is increased by one and returns the last vertex label.

1: **function** BGU_ADD VERTEX()

2: $N_V$.APPEND(EMPTY_LIST())

3: **return** $|N_V| - 1$

## 2.5 Adding an edge

The following algorithms append new edges by providing their label.

► Algorithm 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) \in \text{PREDS}(v)$ and $(v, e) \in \text{SUCCS}(u)$.

1: **function** DGS_ADD EDGE(u, v)

2: $\triangleright$ ind_edges is the number of edges. It is always even. $\triangleleft$

3: $N_u^+$.APPEND(v, ind_edges)

4: $N_v^+$.APPEND($\bar{u}$, ind_edges)
5: $ind_edges \leftarrow ind_edges + 1$

6: $card_edges \leftarrow card_edges + 2$

7: **return** $ind_edges - 1$

► Algorithm 14: Add a new edge and its reverse in the graph for DGF

**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) \in \text{PREDS}(v, vor)$ and $(v, vor, e) \in \text{SUCCS}(u, uor)$.

1: **function** DGF_ADDEdge($u, uor, v, vor$)

2: ▷ $ind_edges$ is the number of edges. It is always even. ◁

3: **if** $uor = f$ **then** ▷ $u_f \to v_f$ or $u_f \to v_r$

4: $N_u^+$.APPEND($v, vor, ind_edges$)

5: **else** ▷ $u_f \leftarrow v_f$ or $u_f \leftarrow v_r$

6: $N_u^-$.APPEND($v, 1 - vor, ind_edges$)

7: **if** $vor = f$ **then** ▷ $u_f \to v_f$ or $u_r \to v_f$

8: $N_v^-$.APPEND($u, uor, ind_edges$)

9: **else** ▷ $u_f \leftarrow v_f$ or $u_r \leftarrow v_f$

10: $N_v^+$.APPEND($u, 1 - uor, ind_edges$)

11: $ind_edges \leftarrow ind_edges + 1$

12: $card_edges \leftarrow card_edges + 2$

13: **return** $ind_edges - 1$

► Algorithm 15: 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) \in \text{PREDS}(v, vor)$ and $(v, vor, e) \in \text{SUCCS}(u, uor)$.

1: **function** BGU_ADDEdge($u, uor, v, vor$)

2: $N_u$.APPEND($v, ind_edges$)

3: $N_v$.APPEND($u, ind_edges$)

4: **if** $u < v$ **then**

5: $Eattr$.APPEND($uor, vor$)

6: **else**

7: $Eattr$.APPEND($1 - vor, 1 - uor$)

8: $ind_edges \leftarrow ind_edges + 1$

9: $card_edges \leftarrow card_edges + 1$
10: **return** ind_edges - 1

## 2.6 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 \mid 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.

1: **function** DGS_DELETE VERTEX($v$)

2: ▷ Delete it from its predecessors

3: $v_rev \leftarrow v + 1$

4: **for** $(u_rev, e) \in N_{v_rev}^+$ **do**

5: $u \leftarrow \overline{u_rev}$

6: $adj\_ind \leftarrow \text{GET\_INDEX}(N_u^+, e, 1)$

7: DELETE_LIST_ENTRY($N_u^+, adj\_ind$)

8: $card_edges \leftarrow card_edges - 2$

9: DELETE($N_{v_rev}^+$)

10: ▷ Delete it from its successors

11: **for** $(w, e) \in N_v^+$ **do**

12: $w_rev \leftarrow \overline{w}$

13: $adj\_ind \leftarrow \text{GET\_INDEX}(N_{w_rev}^+, e, 1)$

14: DELETE_LIST_ENTRY($N_{w_rev}^+, adj\_ind$)

15: $card_edges \leftarrow card_edges - 2$

16: DELETE($N_v^+$)

17: ▷ Delete the whole vertex

18: **if** $v = |N_V^+| - 2$ **then**

19: ▷ It is the last, just pop it and its reverse

20: $N_V^+.\text{POP}($)

21: $N_V^+.\text{POP}($)

22: **else**

23: ▷ Replace it and its reverse by the last and its reverse

24: $N_{v_rev}^+ \leftarrow N_V^+.\text{POP}($)

25: $N_v^+ \leftarrow N_V^+.\text{POP}($)

26: ▷ Update the vertex index
27: **for** `(u_rev, e) ∈ N<sub>v_rev</sub><sup>+</sup> **do**

28: $u \leftarrow \overline{u\_rev}$

29: $adj\_ind \leftarrow \text{GET\_INDEX}(N_u^+, e, 1)$

30: $N_u^+[adj\_ind] \leftarrow (v, e)$

31: **for** `(w_rev, e) ∈ N_v^+ **do**

32: $w \leftarrow \overline{w\_rev}$

33: $adj\_ind \leftarrow \text{GET\_INDEX}(N_w^+, e, 1)$

34: $N_w^+[adj\_ind] \leftarrow (v\_rev, e)$

► Algorithm 17: Delete a vertex for DGF

**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$)

2: ▷ Delete it from its predecessors

3: **for** `(u, uor, e) ∈ N_v^-$ **do**

4: if uor = f then

5: $adj\_ind \leftarrow \text{GET\_INDEX}(N_u^+, e, 1)$

6: DELETE_LIST_ENTRY($N_u^+$, $adj\_ind$)

7: else

8: $adj\_ind \leftarrow \text{GET\_INDEX}(N_u^-, e, 1)$

9: DELETE_LIST_ENTRY($N_u^-$, $adj\_ind$)

10: $card\_edges \leftarrow 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 \leftarrow \text{GET\_INDEX}(N_w^-, e, 1)$

16: DELETE_LIST_ENTRY($N_w^-$, $adj\_ind$)

17: else

18: $adj\_ind \leftarrow \text{GET\_INDEX}(N_w^+, e, 1)$

19: DELETE_LIST_ENTRY($N_w^+$, $adj\_ind$)

20: $card\_edges \leftarrow card\_edges - 2$

21: DELETE($N_v^+$)

22: ▷ Delete the whole vertex

23: **if** $v = |N_v^+| - 1$ **then**

24: ▷ It is the last, just pop it and its reverse

25: $N_v^- \text{.POP()}$
26: $N_V^+.POP()$

27: **else**

28: ▷ Replace it and its reverse by the last and its reverse ▷

29: $N_v^- \leftarrow N_V^+.POP()$

30: $N_v^+ \leftarrow N_V^+.POP()$

31: ▷ Update the vertex index ▷

32: **for** $(u, uor, e) \in N_v^-$ **do**

33: **if** $uor = f$ **then**

34: $adj\_ind \leftarrow GET\_INDEX(N_u^+, e, 1)$

35: $N_u^+[adj\_ind] \leftarrow (v, f, e)$

36: **else**

37: $adj\_ind \leftarrow GET\_INDEX(N_u^-, e, 1)$

38: $N_u^-[adj\_ind] \leftarrow (v, r, e)$

39: **for** $(w, wor, e) \in N_v^+$ **do**

40: **if** $wor = f$ **then**

41: $adj\_ind \leftarrow GET\_INDEX(N_w^-, e, 1)$

42: $N_w^-[adj\_ind] \leftarrow (v, f, e)$

43: **else**

44: $adj\_ind \leftarrow GET\_INDEX(N_w^+, e, 1)$

45: $N_w^+[adj\_ind] \leftarrow (v, r, e)$

**► Algorithm 18: Delete a vertex for BGU**

**Require:** A 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$)

2: **for** $(w, e) \in N_v$ **do**

3: $adj\_ind \leftarrow GET\_INDEX(N_w, e, 1)$

4: DELETE_LIST_ENTRY($N_w$, $adj\_ind$)

5: $card\_edges \leftarrow card\_edges - 1$

6: DELETE($N_v$)

7: ▷ Delete the whole vertex ▷

8: **if** $v = |N_V| - 1$ **then**

9: ▷ It is the last, just pop it ▷

10: $N_V$.POP()

11: **else**

12: ▷ Replace it by the last ▷

13: $N_v \leftarrow N_V$.POP()
14: ▷ Update the vertex index ◁

15: **for** (*w*, *e*) ∈ *N*<sub>*v*</sub> **do**

16:     *adj_ind* ← GET_INDEX(*N*<sub>*w*</sub>, *e*, 1)

17:     *N*<sub>*w*</sub>[*adj_ind*] ← (*v*, *e*)

18:     **if** *w* > *v* **then**

19:         ▷ The last (the greatest) identifier becomes a fewer and breaks
             the identifier order ◁

20:     *Eattr*[*e*] ← (1 - *Eattr*[*e*][1], 1 - *Eattr*[*e*][0])

## 2.7 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 an edge label.

**Ensure:** The edges are deleted for the two vertices.

1: **function** DGS_DELETEEdge(*u*, *v*, *e*) ◁

2:     ▷ Remove *v* from *u*'s successors ◁

3:     *adj_ind* ← GET_INDEX(*N*<sub>*u*</sub><sup>+</sup>, *e*, 1)

4:     DELETE_LIST_ENTRY(*N*<sub>*u*</sub><sup>+</sup>, *adj_ind*) ◁

5:     ▷ Remove *u* from *v*'s predecessors*. ◁

6:     *v_rev* ← $\bar{v}$

7:     *adj_ind* ← GET_INDEX(*N*<sub>*v_rev*</sub><sup>+</sup>, *e*, 1)

8:     DELETE_LIST_ENTRY(*N*<sub>*v_rev*</sub><sup>+</sup>, *adj_ind*) ◁

9:     *card_edges* ← *card_edges* - 2

### ▶ Algorithm 20: Delete an edge for DGF

**Require:** Two vertex labels *u* and *v* associated with their fragment orientation *uor* and *vor*, respectively, and an edge label.

**Ensure:** The edges are deleted for the two vertices.

1: **function** DGF_DELETEEdge(*u*, *uor*, *v*, *vor*, *e*) ◁

2:     ▷ Remove *v* from *u*'s successors ◁

3:     **if** *uor* = *f* **then**

4:         *list* ← *N*<sub>*u*</sub><sup>+</sup>

5:     **else**

6:         *list* ← *N*<sub>*u*</sub><sup>-</sup>
7: $adj\_ind \leftarrow GET\_INDEX(list\_succs, e, 2)$

8: `DELETE_LIST_ENTRY(list\_succs, adj\_ind)`

9: $\triangleright$ *Remove *u* from *v*'s predecessors.* $\triangleleft$

10: **if** *vor* = *f* **then**

11: $list\_preds \leftarrow N_v^-$

12: **else**

13: $list\_preds \leftarrow N_v^+$

14: $adj\_ind \leftarrow GET\_INDEX(list\_preds, e, 2)$

15: `DELETE_LIST_ENTRY(list\_preds, adj\_ind)`

16: $card\_edges \leftarrow card\_edges - 2$

► Algorithm 21: Delete an edge for BGU

**Require:** Two vertex labels *u* and *v* and an edge label.

**Ensure:** The edges are deleted for the two vertices.

1: **function** BGU_DELETEEdge(*u*, *v*, *e*)

2: $\triangleright$ *Remove *v* from *u*'s neighbours.* $\triangleleft$

3: $adj\_ind \leftarrow GET\_INDEX(N_u, e, 1)$

4: `DELETE_LIST_ENTRY(N_u, adj\_ind)`

5: $\triangleright$ *Remove *u* from *v*'s neighbours.* $\triangleleft$

6: $adj\_ind \leftarrow GET\_INDEX(N_v, e, 1)$

7: `DELETE_LIST_ENTRY(N_v, adj\_ind)`

8: $card\_edges \leftarrow card\_edges - 1$

3 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 $|\mathcal{F}| \gg k$ and $|\mathcal{L}| \gg 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 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.
■ Table 4 – Calculus detail of basic operations costs.

$c(x)$ is the cost of operation $x$.

<table><thead><tr><th>Operation</th><th>Cost</th><th>Description</th></tr></thead><tbody><tr><td>$x$</td><td>0</td><td>Memory access</td></tr><tr><td>$x \pm y$</td><td>$1 + c(x) + c(y)$</td><td>Basic operation</td></tr><tr><td>$x \leftarrow y$</td><td>$1 + c(y)$</td><td>Affectation</td></tr><tr><td>$|x|$</td><td>$1 + c(x)$</td><td>Absolute value or element's length</td></tr><tr><td><b>if</b> $x \le y$</td><td>$1 + c(x) + c(y)$</td><td>Conditional</td></tr><tr><td><b>empty_list()</b></td><td>1</td><td>Empty list constructor</td></tr><tr><td><b>list.append(x)</b></td><td>$1 + c(x)$</td><td>Add $x$ to the end of the list</td></tr><tr><td><b>list.pop()</b></td><td>1</td><td>Delete the last element of the list</td></tr><tr><td><b>delete(list)</b></td><td>$|list|$</td><td>Delete all the list</td></tr><tr><td><b>for</b> $i \in [a, b] \text{ do } x \ (b - a + 1) \times c(x)$</td><td></td><td>For loop</td></tr></tbody></table>

■ Table 5 – Algorithmic costs for subfunctions.

$|list|$ is the length of the input list.

<table><thead><tr><th></th><th>Best</th><th>Worst</th><th>Average</th></tr></thead><tbody><tr><td>REV</td><td></td><td>2</td><td></td></tr><tr><td>GET_INDEX</td><td>4</td><td>1 + 5|list|</td><td>1 + 5[<sup>|list|</sup> / 2]</td></tr><tr><td>DELETE_LIST_ENTRY</td><td>3</td><td>4</td><td>4 - <sup>1</sup> / |list|</td></tr></tbody></table>

■ Table 6 – Algorithmic costs of iterating over the neighbours.

<table><caption>(a) Iterating over the predecessors</caption><thead><tr><th></th><th>Best</th><th>Worst</th><th>Average</th></tr></thead><tbody><tr><td>DGS</td><td></td><td>3k + 4</td><td></td></tr><tr><td>DGF</td><td>k + 3</td><td>2k + 3</td><td><sup>3</sup> / 2k + 3</td></tr><tr><td>BGU</td><td>6k + 2</td><td>12k + 2</td><td><sup>17</sup> / 2k + 2</td></tr></tbody></table>

<table><caption>(b) Iterating over the successors</caption><thead><tr><th></th><th>Best</th><th>Worst</th><th>Average</th></tr></thead><tbody><tr><td>DGS</td><td></td><td>k + 2</td><td></td></tr><tr><td>DGF</td><td>k + 3</td><td>2k + 3</td><td><sup>3</sup> / 2k + 3</td></tr><tr><td>BGU</td><td>6k + 2</td><td>12k + 2</td><td><sup>17</sup> / 2k + 2</td></tr></tbody></table>
<table><caption>■ Table 7 – Algorithmic costs of adding a vertex or an edge.</caption><caption>(a) Adding a vertex</caption><thead><tr><th></th><th>Best</th><th>Worst</th><th>Average</th></tr></thead><tbody><tr><th>DGS</th><td>5</td><td></td><td></td></tr><tr><th>DGF</th><td>5</td><td></td><td></td></tr><tr><th>BGU</th><td>3</td><td></td><td></td></tr></tbody></table><table><caption>(b) Adding an edge</caption><thead><tr><th></th><th>Best</th><th>Worst</th><th>Average</th></tr></thead><tbody><tr><th>DGS</th><td></td><td>11</td><td></td></tr><tr><th>DGF</th><td>9</td><td>11</td><td>10</td></tr><tr><th>BGU</th><td>9</td><td>11</td><td>10</td></tr></tbody></table>

<table><caption>■ Table 8 – Algorithmic costs of deleting a vertex.</caption><thead><tr><th></th><th>Best</th><th>Worst</th><th>Average</th></tr></thead><tbody><tr><th>DGS</th><td>28k + 6</td><td>20k<sup>2</sup> + 36k + 8</td><td>10<sup>|F|-1</sup>⁄<sub>|F|</sub>k<sup>2</sup> + 36<sup>|F|-2</sup>⁄<sub>|F|</sub>k + <sup>1-|F|</sup>⁄<sub>|F|</sub></td></tr><tr><th>DGF</th><td>24k + 4</td><td>20k<sup>2</sup> + 28k + 6</td><td>10<sup>|F|-1</sup>⁄<sub>|F|</sub>k<sup>2</sup> + 28<sup>|F|-4</sup>⁄<sub>|F|</sub>k + 4<sup>|F|-1</sup>⁄<sub>|F|</sub></td></tr><tr><th>BGU</th><td>22k + 3</td><td>40k<sup>2</sup> + 38k + 7</td><td>20<sup>|F|-1</sup>⁄<sub>|F|</sub>k<sup>2</sup> + 24<sup>|F|-1</sup>⁄<sub>|F|</sub>k + <sup>11</sup>⁄<sub>2</sub> <sup>|F|-14</sup>⁄<sub>|F|</sub></td></tr></tbody></table>

<table><caption>■ Table 9 – Algorithmic costs of deleting an edge.</caption><thead><tr><th></th><th>Best</th><th>Worst</th><th>Average</th></tr></thead><tbody><tr><th>DGS</th><td>21</td><td>10k + 17</td><td>5k + 17 - <sup>2</sup>⁄<sub>k</sub></td></tr><tr><th>DGF</th><td>22</td><td>10k + 18</td><td>5k + 18 - <sup>2</sup>⁄<sub>k</sub></td></tr><tr><th>BGU</th><td>18</td><td>20k + 14</td><td>10k + 14 - <sup>1</sup>⁄<sub>k</sub></td></tr></tbody></table>

## 4 Memory and time cost comparisons

Table 10 provides memory consumption and elementary operation time cost comparisons between *DGS*, *DGF* and *BGU*. Each row presents a score ranging from one to three stars for the three implementations, simplifying the comparison under specific conditions. In the same row, the cost functions of *DGS*, *DGF* and *BGU* are all polynomial functions with the same degree. The polynomial function with the smallest coefficient on the highest degree is used as the basis for comparison. Let $f_{min}$ be this function, and *g* and *h* be the two others. Let $c = \lim_{x \to \infty} \frac{g(x)}{f_{min}(x)}$ and $d = \lim_{x \to \infty} \frac{h(x)}{f_{min}(x)}$ be the results of the asymptotic comparisons. The scores are determined whether *c* and *d* belong to the following intervals: *** [1, 1.2]; ** [1.2, 2]; ★ [2, +∞[.
■ Table 10 – Comparison between DGS, DGF 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.

<table><thead><tr><th></th><th>DGS</th><th>DGF</th><th>BGU</th></tr></thead><tbody><tr><td colspan="4"><strong>Memory consumption</strong></td></tr><tr><td></td><td>*</td><td>*</td><td>***</td></tr><tr><td colspan="4"><strong>Iterating over the predecessors</strong></td></tr><tr><td>Best</td><td>*</td><td>***</td><td>*</td></tr><tr><td>Worst</td><td>**</td><td>***</td><td>*</td></tr><tr><td>Average</td><td>*</td><td>***</td><td>*</td></tr><tr><td colspan="4"><strong>Iterating over the successors</strong></td></tr><tr><td>Best</td><td>***</td><td>***</td><td>*</td></tr><tr><td>Worst</td><td>***</td><td>*</td><td>*</td></tr><tr><td>Average</td><td>***</td><td>**</td><td>*</td></tr><tr><td colspan="4"><strong>Adding a vertex</strong></td></tr><tr><td>Best</td><td>**</td><td>**</td><td>***</td></tr><tr><td>Worst</td><td>**</td><td>**</td><td>***</td></tr><tr><td>Average</td><td>**</td><td>**</td><td>***</td></tr><tr><td colspan="4"><strong>Adding an edge</strong></td></tr><tr><td>Best</td><td>**</td><td>***</td><td>***</td></tr><tr><td>Worst</td><td>***</td><td>***</td><td>***</td></tr><tr><td>Average</td><td>***</td><td>***</td><td>***</td></tr><tr><td colspan="4"><strong>Deleting a vertex</strong></td></tr><tr><td>Best</td><td>**</td><td>***</td><td>***</td></tr><tr><td>Worst</td><td>***</td><td>***</td><td>*</td></tr><tr><td>Average</td><td>***</td><td>***</td><td>*</td></tr><tr><td colspan="4"><strong>Deleting an edge</strong></td></tr><tr><td>Best</td><td>***</td><td>**</td><td>***</td></tr><tr><td>Worst</td><td>***</td><td>***</td><td>*</td></tr><tr><td>Average</td><td>***</td><td>***</td><td>*</td></tr></tbody></table>

5 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