THÈSE / UNIVERSITÉ DE RENNES 1
sous le sceau de l'Université Bretagne Loire

pour le grade de
DOCTEUR DE L'UNIVERSITÉ DE RENNES 1

Mention : Informatique

Ecole doctorale MATISSE

présentée par

Antoine Limasset

préparée à l'unité de recherche 6074 IRISA
INSTITUT DE RECHERCHE EN INFORMATIQUE ET
SYSTEMES ALEATOIRES
ISTIC

Nouvelles approches
pour l'exploitation
des données de
sequençage haut débit

Thèse soutenue à Rennes
le 12 Juillet 2017

devant le jury composé de :

Hélène TOUZET

Directeur de recherche, CRYSTAL / rapporteur

Thierry LECROQ

Directeur de recherche, LITIS / rapporteur

Christophe KLOPP

Ingenieur de recherche, BIA / examinateur

Vincent LACROIX

Maitre de conférence, LBBE / examinateur

Gregory KUCHEROV

Directeur de recherche, LIGM / examinateur

Rumen ANDONOV

Professeur, IRISA / examinateur

Dominique LAVENIER

Directeur de recherche, IRISA / directeur de thèse

Pierre PETERLONGO

Chargé de recherche, IRISA / co-directeur de thèse
<figure><img src="image_1.png" /></figure>
# Acknowledgments

## Administration

Je voudrais tout d'abord remercier l'école doctorale Matisse, l'Université Rennes 1, IRISA et INRIA ainsi que leur personnels sans qui cette thèse et plus généralement cet environnement de recherche et d'enseignement scientifique n'existeraient tout simplement pas. Je remercie doublement les personnels d'administration à qui j'ai eu affaire et à qui j'ai dû donner du fil à retordre. Une pensée spéciale à Marie Le Roïc qui n'a jamais paru effrayée par mon organisation plus que doûteuse et à l'administration de l'école doctorale pour avoir réussi à me faire remplir un dossier correctement pour la première fois en 26 ans.

**Recherche** Je voudrais bien évidemment remercier tous mes collaborateurs. Merci à Eric Rivals de m'avoir jeté dans le monde de la bioinformatique et de m'y avoir donné goût pendant mon stage de L3. Merci à Dominique Lavenier de m'avoir non seulement conseillé mais guidé vers un stage puis une thèse entre de très bonnes mains. Merci a Rayan Chikhi et Paul Medvedev de m'avoir lancé dans le petit monde du graphe de De Bruijn, je me demande souvent quand est-ce que je vais pouvoir en sortir. Je remercie ici, avec quelques années de retard, Rayan de m'avoir accueilli en ami plus qu'en encadrant, je garde de très bons souvenirs, peut-être un peu trop français, de ce voyage aux Etats-Unis. Je voudrais spécialement remercier ces premiers encadrants pour leur bienveillance qui m'a permis de me conforter dans ma capacité à contribuer à la formidable construction qu'est la recherche, je ne suis plus tout à fait le même depuis.

Je ne vais pas prendre assez de place pour remercier Pierre Peterlongo pour m'avoir encadré et suivi pendant plus de trois ans. Je pense avoir été une thésard difficile et malgré nos disjonctions tu m'as toujours laissé libre dans mon fonctionnement. Je ne pense pas que j'aurais pu trouver une meilleure personne pour cette aventure. Je trouve que nous avons formé une belle paire et je suis impressionné par la vitesse à laquelle les choses avancent.

Merci également à tout le groupe Colib'read, plein de grandes rencontres ont eu lieu pour moi lors de ces réunions. Si je garde la vision d'un groupe de recherche comme un groupe d'amis c'est un peu grâce à vous.

Je voulais aussi remercier Romain Feron, je pourrais me targuer sûrement longtemps que mon premier encadrement se soit si bien passé. J'espère que tu en garderas de bons souvenirs et que tu continueras dans ta lancée.
Petite pensée pour Thomas Rokicki et Dianne Puges qui ont trouvé l'assemblage génomique suffisamment intéressant pour participer à un petit projet avec moi. Je regrette un peu que celui-ci n'ait jamais abouti mais il n'est jamais trop tard.

Merci à Jean-François Flot de m'avoir invité à Bruxelles et d'être aussi enthousiaste. Je nous souhaite beaucoup de succès à l'avenir.

Merci à tous mes co-bureaux de m'avoir supporté, et je souhaite à Marie Chevallier une bonne et heureuse année. Plus généralement je remercie toute l'équipe GenScale de m'avoir acueilli et pour cette bonne ambiance au travail qui n'a tout simplement pas de prix.

### These

Je remercie très solennellement Yannick Zakowski Hugo Bazille Camille Marchet et Pierre Peterlongo d'avoir accompli un travail digne des travaux d'Asterix: corriger ma prose. Je regrette presque de ne pas avoir compté le nombres de fautes corrigées par vos soins mais sachez que votre tâche a été colossale. Une pensée pour Pierre qui a dû passer plusieurs nuits blanches à relire les différentes versions de celle-ci, à se fouler le bras en rayant des centaines de pages. Un très très grand merci à Camille pour l'organisation de mon pot de thèse et des festivités, si j'avais moins procrastiné j'aurais pu davantage mettre la main a la pâte et j'aurais pu mettre Paris en bouteille.

### Enseignement

J'ai eu le plus grand des plaisirs lors des missions d'enseignement qui m'ont été accordées par l'ISTIC. Je ne peux remercier tous mes élèves mais sait-on jamais, merci de votre interêt (parfois) et de votre sympathie (souvent), vous avez été super. Je voudrais également profiter de cette occasion pour remercier quelques-un de mes professeurs. Merci a Mohamed Rebbal de m'avoir fait comprendre l'essence de la résolution de problèmes. Merci à Nicolas De Granrut et Gérard Debeaumarché pour m'avoir fait confiance et permis d'accéder aux études supérieures auxquelles j'aspirais. Ma vie serait bien différente sans votre soutien et je travaille à rester digne de celui-ci.

### Amis

Ces années à Rennes m'ont vu changer d'une manière saisissante et je vais remercier ici essentiellement des Rennais qui m'ont entouré lors de mes péripéties dans ce beau pays, que personne n'en prenne ombrage. Merci à Justine de m'avoir suivi jusque ici et d'avoir été ma première experience de vie commune. J'espère que tu garderas une place spéciale à nos élucubrations et diverses aventures farfelues en pays breton, elles ont pour ma part toute leur place dans la boite des bons souvenirs. Merci a Simon Rokicki et Bastien Pasdeloup d'avoir les premiers accueilli un Cachannais expatrié et d'avoir fait de moi un des vôtres. Je remercie tout les membres des Projet C2Silicium et UberML, pour avoir contribué au lignes les plus hilarantes de mon CV. Encore une fois je ne vais pas prendre assez de temps pour remercier suffisament la Complexiteam. Merci pour tout ces super moments, je m'excuse officiellement pour mon Warwick Top et j'accuse le mauvais matchup. Plus sérieusement quoi que l'on vous dise, sachez-le, dans mon coeur, nous sommes au moins Diamant. Merci à Joris Lamare d'avoir gardé le contact malgré la distance et d'être toujours là pour accueillir des voyageurs déglingués et pour nous proposer les meilleurs plan de la Terre, pas moins. Un des plus grand merci à Gaëtan Benoit pour avoir été mon plus proche collègue pendant toutes ces années. Pour ne citer qu'une de nos nombreuses aventures, merci d'avoir porté le projet de court-métrage sur l'assemblage génomique et sache que, quoi que le président du jury en dise, notre film
reste le meilleur. Je vais également remercier les Doudous, pour tout votre soutien, vous n'avez jamais répondu absents lorsque j'avais besoin d'aide, je serais etonné de vous en avoir rendu le dixième. Je n'oublie pas un remerciement spécial à Simon pour son aide et ses discussions précieuses sur la compilation et les serveurs Minecraft. Merci à Olivier Ruas et Hugo Bazille pour avoir essayé de me maintenir en forme le long de cette thèse. Merci à Gurvan Cabon d'avoir été le coloc le moins chiant du monde. Dur à classer je vais remercier également mes amis à quatre pattes. Merci Belial pour avoir écrit plus que moi sur mon document de thèse (notamment le fameux éedFFFFFFFFFFFFFFFFFFFFFFFFFFFF). Merci à Lux d'avoir été la plus gentille des chattes et mon premier animal. Merci à Jelly d'être un si bon petit tambour. Merci a BBhash de m'avoir laissé des doigts pour travailler.

**Famille** Je voudrais remercier mes parents sans qui rien de tout cela n'aurait été possible. Plus sérieusement, je vous remercie de m'avoir donné les clefs pour en arriver là. Merci pour tout votre amour et tout votre soutien. Tout ce que vous m'avez apporté a fait de moi quelqu'un d'entier et pour cela je vous remercie infiniment. Je voudrais finalement remercier ici Camille Marchet que je considère aujourd'hui comme la meilleure chose qui me soit arrivé. Si un des remerciements n'est pas assez développé c'est celui-ci. Je n'aurais pas été jusqu'ici sans toi, merci pour ces bons moments.

**Divers** A mon grand regret ceux-ci ne pourront jamais lire ces remerciements. Merci a Nadine et Griffin, mes deux ordinateurs de travail, pour leur fidélité. Que Nadine me pardonne mon inattention, j'espère que tu vis heureux avec ton nouveau propriétaire qui a dû avoir bien du mal à réinstaller MacOS sur un Linux. Merci au Genocluster de m'avoir supporté tout ce temps, pardonne-moi si je t'en ai parfois trop demandé. Pour finir sur une note sérieuse je voudrais remercier les projets suivants et leurs collaborateurs : Libreoffice Linux Scihub Github Travis Gnu.
# Abstract

## Novel approaches for the exploitation of high throughput sequencing data

In this thesis we discuss computational methods to deal with DNA sequences provided by high throughput sequencers. We will mostly focus on the reconstruction of genomes from DNA fragments (genome assembly) and closely related problems. These tasks combine huge amounts of data with combinatorial problems. Various graph structures are used to handle this problem, presenting trade-off between scalability and assembly quality.

This thesis introduces several contributions in order to cope with these tasks. First, novel representations of assembly graphs are proposed to allow a better scaling. We also present novel uses of those graphs apart from assembly and we propose tools to use such graphs as references when a fully assembled genome is not available. Finally we show how to use those methods to produce less fragmented assembly while remaining tractable.
# Nouvelles approches pour l'exploitation des données de séquencage haut débit

Cette thèse a pour sujet les méthodes informatiques traitant les séquences ADN provenant des séquenceurs haut débit. Nous nous concentrons essentiellement sur la reconstruction de génomes à partir de fragments ADN (assemblage génomique) et sur des problèmes connexes. Ces tâches combinent de très grandes quantités de données et des problèmes combinatoires. Différentes structures de graphe sont utilisées pour répondre à ces problèmes, présentant des compromis entre passage à l'échelle et qualité d'assemblage.

Ce document introduit plusieurs contributions pour répondre à ces problèmes. De nouvelles représentations de graphes d'assemblage sont proposées pour autoriser un meilleur passage à l'échelle. Nous présentons également de nouveaux usages de ces graphes, différents de l'assemblage, ainsi que des outils pour utiliser ceux-ci comme références dans les cas où un génome de référence n'est pas disponible. Pour finir nous montrons comment utiliser ces méthodes pour produire un meilleur assemblage en utilisant des ressources raisonnables.
<figure><img src="image_2.png" /></figure>
# Contents

<table><tr><td><strong>Acknowledgments</strong></td><td style="text-align: right;"><strong>1</strong></td></tr><tr><td><strong>Abstract</strong></td><td style="text-align: right;"><strong>3</strong></td></tr><tr><td><strong>Table of contents</strong></td><td style="text-align: right;"><strong>4</strong></td></tr><tr><td><strong>1 Introduction</strong></td><td style="text-align: right;"><strong>7</strong></td></tr><tr><td style="padding-left: 2em;">1.1 DNA, genomes and sequencing</td><td style="text-align: right;">8</td></tr><tr><td style="padding-left: 2em;">1.2 Genome assembly</td><td style="text-align: right;">11</td></tr><tr><td style="padding-left: 4em;">1.2.1 Greedy</td><td style="text-align: right;">14</td></tr><tr><td style="padding-left: 4em;">1.2.2 Overlap Layout Consensus</td><td style="text-align: right;">15</td></tr><tr><td style="padding-left: 4em;">1.2.3 De Bruijn graphs</td><td style="text-align: right;">20</td></tr><tr><td style="padding-left: 2em;">1.3 Assembly hardness</td><td style="text-align: right;">31</td></tr><tr><td style="padding-left: 4em;">1.3.1 Repeats</td><td style="text-align: right;">31</td></tr><tr><td style="padding-left: 4em;">1.3.2 Scaffolding</td><td style="text-align: right;">33</td></tr><tr><td style="padding-left: 4em;">1.3.3 Multiple genomes</td><td style="text-align: right;">34</td></tr><tr><td style="padding-left: 2em;">1.4 Outline of the thesis</td><td style="text-align: right;">36</td></tr><tr><td><strong>2 Handling assembly</strong></td><td style="text-align: right;"><strong>37</strong></td></tr><tr><td style="padding-left: 2em;">2.1 The assembly burden</td><td style="text-align: right;">38</td></tr><tr><td style="padding-left: 2em;">2.2 Overlap graph scalability</td><td style="text-align: right;">41</td></tr><tr><td style="padding-left: 2em;">2.3 The scaling story of the de Bruijn graph representation</td><td style="text-align: right;">43</td></tr><tr><td style="padding-left: 4em;">2.3.1 Generic graphs</td><td style="text-align: right;">44</td></tr><tr><td style="padding-left: 4em;">2.3.2 Theoretical limits</td><td style="text-align: right;">45</td></tr><tr><td style="padding-left: 4em;">2.3.3 Kmer Counting</td><td style="text-align: right;">46</td></tr><tr><td style="padding-left: 4em;">2.3.4 Probabilistic de Bruijn graphs</td><td style="text-align: right;">46</td></tr><tr><td style="padding-left: 4em;">2.3.5 Navigational data structures</td><td style="text-align: right;">47</td></tr><tr><td style="padding-left: 4em;">2.3.6 Massively parallel assembly</td><td style="text-align: right;">48</td></tr><tr><td style="padding-left: 2em;">2.4 Efficient de Bruijn graph representation</td><td style="text-align: right;">49</td></tr><tr><td style="padding-left: 4em;">2.4.1 Compacted de Bruijn graph</td><td style="text-align: right;">49</td></tr><tr><td style="padding-left: 4em;">2.4.2 De Bruijn graph construction</td><td style="text-align: right;">49</td></tr><tr><td style="padding-left: 4em;">2.4.3 Unitig enumeration in low memory: BCALM</td><td style="text-align: right;">51</td></tr><tr><td style="padding-left: 4em;">2.4.4 Assembly in low memory using BCALM</td><td style="text-align: right;">54</td></tr><tr><td style="padding-left: 2em;">2.5 Efficient de Bruijn graph construction</td><td style="text-align: right;">56</td></tr><tr><td style="padding-left: 4em;">2.5.1 BCALM2 algorithm</td><td style="text-align: right;">56</td></tr><tr><td style="padding-left: 4em;">2.5.2 Implementation details</td><td style="text-align: right;">57</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>2.5.3 Large genome de Bruijn graphs</td>
      <td>59</td>
    </tr>
    <tr>
      <td>2.6 Indexing large sets</td>
      <td>62</td>
    </tr>
    <tr>
      <td><b>3 The de Bruijn graph as a reference</b></td>
      <td><b>81</b></td>
    </tr>
    <tr>
      <td>3.1 Genome representations</td>
      <td>82</td>
    </tr>
    <tr>
      <td>3.1.1 Reference sequences</td>
      <td>82</td>
    </tr>
    <tr>
      <td>3.1.2 Genome graphs</td>
      <td>84</td>
    </tr>
    <tr>
      <td>3.1.3 De Bruijn graphs as references</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.2 Read mapping on the de Bruijn graph</td>
      <td>91</td>
    </tr>
    <tr>
      <td>3.2.1 An efficient tool for an NP-Complete problem</td>
      <td>91</td>
    </tr>
    <tr>
      <td>3.2.2 Mapping refinements</td>
      <td>92</td>
    </tr>
    <tr>
      <td>3.3 De novo, reference guided, read correction</td>
      <td>93</td>
    </tr>
    <tr>
      <td>3.3.1 Limits of kmer spectrum</td>
      <td>93</td>
    </tr>
    <tr>
      <td>3.3.2 Reads correction by de Bruijn graph mapping</td>
      <td>94</td>
    </tr>
    <tr>
      <td><b>4 Improving an assembly paradigm</b></td>
      <td><b>111</b></td>
    </tr>
    <tr>
      <td>4.1 Assembly challenges</td>
      <td>112</td>
    </tr>
    <tr>
      <td>4.1.1 Two assembly paradigms</td>
      <td>112</td>
    </tr>
    <tr>
      <td>4.1.2 Why read length matters: polyploid assembly</td>
      <td>112</td>
    </tr>
    <tr>
      <td>4.1.3 Taking the best of both worlds</td>
      <td>113</td>
    </tr>
    <tr>
      <td>4.2 Building very high order de Bruijn graph</td>
      <td>116</td>
    </tr>
    <tr>
      <td>4.2.1 Successive mapping strategy</td>
      <td>116</td>
    </tr>
    <tr>
      <td>4.2.2 Beyond read length</td>
      <td>121</td>
    </tr>
    <tr>
      <td>4.3 Assembly results</td>
      <td>123</td>
    </tr>
    <tr>
      <td>4.3.1 Haploid genome</td>
      <td>123</td>
    </tr>
    <tr>
      <td>4.3.2 Diploid genome</td>
      <td>124</td>
    </tr>
    <tr>
      <td>4.3.3 Future works</td>
      <td>125</td>
    </tr>
    <tr>
      <td><b>5 Conclusions and perspectives</b></td>
      <td><b>129</b></td>
    </tr>
    <tr>
      <td>5.1 Conclusion</td>
      <td>130</td>
    </tr>
    <tr>
      <td>5.2 Proposed methods</td>
      <td>132</td>
    </tr>
    <tr>
      <td>5.3 Future of sequencing</td>
      <td>133</td>
    </tr>
    <tr>
      <td>5.3.1 Third generation sequencing</td>
      <td>133</td>
    </tr>
    <tr>
      <td>5.3.2 Long range short reads sequencing</td>
      <td>133</td>
    </tr>
    <tr>
      <td>5.4 Future of this work and perspectives</td>
      <td>134</td>
    </tr>
    <tr>
      <td><b>Bibliography</b></td>
      <td><b>148</b></td>
    </tr>
    <tr>
      <td>Table of contents</td>
      <td></td>
    </tr>
  </tbody>
</table>
# Chapter 1

## Introduction

"People say they want simple things, but they don't. Not always."

John D. Cook
In this chapter we first introduce the global context of DNA sequencing and genome assembly. Then we provide a high level description of the main methods used in this field. Thereafter we describe the limits and challenges faced nowadays in genome assembly. We finish this chapter by an outline of the thesis.

## 1.1 DNA, genomes and sequencing

Deoxyribonucleic acid or DNA is a molecule that stores biological information used in the functioning of all known living organisms. Most DNA molecules consist of two polynucleotides strands coiled around each other to form a double helix, composed of simpler component called nucleotides (Figure 1.1). Each nucleotide, or base, can be either identified as adenine (A), cytosine (C), guanine (G) or thymine (T). The bases of the two separate strands are bound together, according to base pairing rules: A with T and C with G, to make double-stranded DNA. The sequence of these four bases along the backbone encodes biological information. DNA is usually stored directly in the cytoplasm for prokaryotes, and in the nucleus and different organelles in eukaryotic cells. Depending on organisms, DNA can be circular or not, and characterize organisms through structures called chromosomes. The whole information within the DNA molecule of an organism is called its genome. Conceptually a genome can be represented as a word or a set of words on the alphabet {A, C, G, T}. Genomes are usually large sequences, a human genome has more than 3 billions bases arranged into 46 chromosomes. In certain viruses, the genome can be encoded in a closely related molecule, the RNA, instead of the DNA. Importantly, there exists a range of redundancy of the genomic information across the living. Some cells store their genomic information into a single set of unpaired chromosomes (haploidy) while others (like humans) have two copies of each chromosome (diploidy). Some species have even more than two copies (polyploidy).

Since DNA discovery in 1953 [1], genome study has shown itself to have huge implications in both academical and industrial fields like agronomy, medicine and ecology. Indeed the knowledge of the genome sequences gives a tremendous access to living organisms characteristics and properties, and is now commonly at the core of biology studies. Sequencing is the operation that consists in determining the bases sequence of a DNA molecule and to encode it on a numerical support for its analysis. Being able to access and study the Human genome [2] [3] is considered as one of the major milestone in scientific history. But this genomic information is partial and imperfect, as no sequencer is yet able to directly output a completely sequenced genome. Thus the apparition of sequencing technologies has created a whole field in computational biology to handle this incredibly discovery-promising data. Several kind of sequencing technologies exist, but they share common properties:

* They only produce fragments of our genome ("reads")

* The locations of the fragments are unknown

* The fragments may contain errors

We call those genomic fragments/substrings "reads" to reflect the fact that they were "read" from the genome. Because of the DNA structure, both strands are present and sequenced. The strands are bases complemented (A to T and C to G) and read in the
opposite way. A genome containing ACCTGC therefore may present reads as CCTG or CAGG: CAGG being the reverse-complement of CCTG read from the opposite strand. Several methods have been proposed to sequence DNA based on a wide range of technologies that will not be described here. The differences between sequencing technologies are essentially errors and read length distribution. We can define an error rate of a read by the ratio of the number of incorrectly sequenced bases by the size of the read.

The sequencing technologies are therefore categorized in three generations according to these characteristics and order of appearance.

* Sanger sequencing [4] was the first method available. It produces sequences of some kilo-bases length with a high accuracy (error rate around 1%).

* Next generation sequencing (NGS), often referred to as Illumina/Solexa sequencing [5], is the most broadly used technology. It produces shorter reads (hundreds of bases) with high accuracy (< 1%). This technology presents an order of magnitude higher throughput and cheaper sequencing than the former technology.

* Third generation sequencing (TGS) is the last generation of sequencing technology, which includes Single Molecule Real Time sequencing [6] and Nanopore sequencing [7], producing very long sequences, up to hundreds kilo-bases. However, they exhibit a very high error rate (up to 30%).

In addition to those properties, each method may show different biases due to the protocol employed. We can mention that in the short reads from NGS, some regions rich in G/C bases are less covered [8] and read ends present higher error rate [9].

Sanger sequencing is almost no longer used because of its cost and low throughput. As for the third sequencing generation, the different technologies are extremely recent

<figure><img src="image_3.png" /><figcaption>Figure 1.1: DNA representation.</figcaption></figure>
and are still currently being developed in a fast moving environment. Nowadays, NGS
remains the most popular technique and we will focus on this type of sequencing in this
document. However TGS will eventually become widely used and the present work was
partially designed to open onto the usage of such data.

For all existing technologies, the main challenge comes from the lack of information
about the origin of a read. Each read could come from any strand in any position of any
DNA molecule introduced in the sequencer. This absence of context and the small size
of the sequences obtained, relatively to a genome size, make it difficult to use reads as
such. Ideally, we would need the access to the underlying genomes in their entirety.

Since the beginning of sequencing of DNA molecules, genomes are produced by struc-
turing and ordering reads information. Then these reconstructed genomes can be used
as references. Reference genomes are the best insight we have about the one-dimensional
organization of information in living cells. They give access not only to the gene se-
quences that lead to proteins, but also to flanking sequences that altogether impact the
functioning of living beings [10]. They also reveal the inner organization of the genome
such as genes relative positions or chromosomes structure. Helping understanding the
genomes and organisms evolution, as well as how all the living is ruled by the encoded
information. Besides, reference genomes can be seen as an entry point for biologists to
use other kinds of data. For instance, they may add information about the known genes
positions and functions to annotate the genome [11] [12].

Reference genome reconstruction is therefore crucial in various domains where raw,
out of context reads are unusable. The task of reordering the reads to recompose the
sequenced genome is called genome assembly. Tools designed for this task, called assem-
blers, have to make no *a priori* hypothesis over the location or the strand of each read
and try to reconstruct the original sequences by ordering the reads relatively. As it will
be detailed further, genome assembly is especially complex as the bases distribution is far
from being uniform. Genomes present specific patterns such as large repeated sequences
(repeats), regions with very specific distributions of nucleotides or extremely repeated
sequences of nucleotides. Such patterns make genomes different from a uniformly dis-
tributed sequence of nucleotides. [13] shows that a human genome is largely constituted
of repeated sequences of significant lengths.

**Reference genomes used in this document**

In this document we will present various results based on different genomes. For the
sake of consistency, we choose a small number of well known and well studied genomes.
The first one is the genome of the Escherichia coli (*E. coli*). *E. coli* is a bacteria with a
genome of 4.6 megabase pairs constituted of one circular chromosome. The second one
is the genome of Caenorhabditis elegans (*C. elegans*). *C. elegans* is a nematode and was
the first multicellular organism to have its whole genome sequenced. Its genome counts
6 chromosomes and is 100 megabase pairs long. Pictures of the two organisms are shown
on Figure 1.2. The third genome is the human reference genome. The version used
was GRCh38, accessible at https://www.ncbi.nlm.nih.gov/grc/human. The genome
counts 23 chromosomes and is 3,234 megabase pairs long.
<figure><img src="image_4.png" /><figcaption>Figure 1.2: *E. coli* bacteria and *C. elegans* worm from https://upload.wikimedia.org/wikipedia/commons/3/32/EscherichiaColi_NIAID.jpg and http://www.socmucimm.org/wp-content/uploads/2014/06/C.-elegans.jpg</figcaption></figure>

"DNA, genomes and sequencing" core messages:

* Genomes are large sequences of "ACTG", whose knowledge is essential to biological studies

* To access this information, we use machines called sequencers

* We are not able to obtain a whole genome directly out of the sequencers but only "reads" that can be seen as fragments of the genome

* Reads are way smaller than the genome and they contain errors

* The task of recovering the original genome from the "reads" is called genome assembly

## 1.2 Genome assembly

Two kinds of genome assembly may be distinguished: reference guided assembly and *de novo* assembly. The reference guided assembly consists in the assembly of a genome when we already have a reference for the species of the individual sequenced. We expect the new genome to be very close to the reference and we are interested in the differences between the individuals. This type of assembly is much easier because we only have to find the differences between the two genome sequences and we can mimic the reference genome to order the reads. Reference guided assemblers as STAGE [14] and Chip-seq [15] consist in two main steps:

* Reads Alignment on the reference

* Consensus between mapped sequences

If this way of proceeding makes the assembly step easier and much less costly, it can seem unsatisfying for two reasons. First because of the biases that the method present. We make the prior hypothesis that the genome to assemble is very close to the reference.
This may mislead the assembly onto something too similar to the reference. Secondly the method is obviously not self-sufficient since a reference needs a prior reference to be constructed. We do not have access to many references genome as it can be seen on Figure 1.3 even if this number should increase dramatically in the next decades. For these reasons we will focus on *de novo* (without reference) assembly in this document.

Since we do not know the original position of the reads, we will try to position them with respect to each others. For this we mainly rely on the notion of "overlaps" between reads. We say that two sequences A and B overlap if a suffix of A is equal to a prefix of B. In some cases we might consider that A and B overlap if a suffix of A is similar to a prefix of B according to some distance. The main intuition being: "If a read A overlaps a read B then it is likely that A and B are consecutive in the genome" (if the overlap is large enough to be significant). Following this intuition, assembly consists in searching for reads that overlap each others. The larger the overlap, the more significant it is, as overlaps of few bases may be spurious. In Figure 1.4 the first and fifth reads share an overlap of 2 bases (GC) although the two reads do not come from the same location of the genome in gray.

Putting aside specific biases, we can depict a sequencing experiment as a uniform distribution of reads along the genome as a first approximation. Large overlaps then mean redundancy in the reads, as bases within an overlap will be present in reads sharing this overlap. In order to obtain large overlaps between sequences, we rely on a high redundancy in the sequencing dataset. The notion of coverage (or depth) of a genome by a sequencing is often used to quantify this redundancy. For a genome of 4 millions

<figure><img src="image_5.png" /><figcaption>Figure 1.3: Number of finished and unfinished genomes in gold database (https://gold.jgi.doe.gov/statistics). Unfinished (or "draft") genomes may be represented as a set of widely fragmented sequences.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 1.4: Intuition of how assembly is possible. The first set of sequences on the left are the reads to be assembled. Overlaps between those reads are computed before we can compact the reads according to the large overlaps found. Compaction between two reads AO and OB according to the overlap O is obtained by concatenating AO with the suffix B excluding O (AAAATCGAGC and CGAGCGCGGA become AAAATCGAGCGCGGA). In this toy example, in a first step, we compact the yellow overlaps and obtain longer sequences. In the second step we compact the cyan overlaps and get the original sequenced genome. Double strand aspect of the DNA is not considered here.</figcaption></figure>

<table><caption>Table 1.1: Expected missing fraction of the genome according to the sequencing depth. Source: adapted from http://www.genome.ou.edu/poisson_CALC.html.</caption><thead><tr><td>Genome coverage</td><td>% Genome not sequenced</td><td>% Genome sequenced</td></tr></thead><tbody><tr><td>0.25</td><td>78</td><td>22</td></tr><tr><td>0.5</td><td>61</td><td>39</td></tr><tr><td>0.75</td><td>47</td><td>53</td></tr><tr><td>1</td><td>37</td><td>63</td></tr><tr><td>2</td><td>14</td><td>86</td></tr><tr><td>5</td><td>0.6</td><td>99.4</td></tr><tr><td>10</td><td>0.0005</td><td>99.995</td></tr></tbody></table>

bases, a sequencing of 4 millions reads of 100 base pairs will present a coverage of 100X because it contains 100 times more nucleotides than the genome. A high coverage is an important factor for at least 3 reasons. The higher the coverage:

* The larger the overlaps between reads are (on average)

* The less chance we have to miss regions of the genome because they are not sequenced (Table 1.1)

* The easier it is to deal with sequencing errors

With a high degree of redundancy we will be able via statistical methods to detect stochastic errors and remove them. Some assembly strategies try to correct reads before assembling them.
<figure><img src="image_7.png" /><figcaption>Figure 1.5: Example of greedy assembly. The overlaps between the first and last reads (AAAATCGA and ATCGATTT) is the largest and therefore compacted first in a greedy manner. This first compaction produces AAAATCGATTT. Then the largest overlap is between the two other reads (CGAGCGCG and GCGGATCG) that are compacted into CGAGCGCGGATCG. The assembler produced two misassembled sequences that are not the shortest common superstring.</figcaption></figure>

We can distinguish three main families of assemblers: "Greedy", "Overlap Layout Consensus" and "de Bruijn graph" as detailed in the three following sections.

### 1.2.1 Greedy

This family of assemblers is the conceptually simpler. The idea is to find a shortest common superstring (SCS) of a set of sequences. Given a set of reads, a SCS is a string T of minimal size such that every read is a substring of T. Since finding the shortest common superstrings is an NP-complete problem [16], greedy assemblers, as their name suggests, apply greedy heuristics. The heuristic performs the compaction of the largest overlap if a read overlaps with several reads. The algorithm can be outlined by:

* Index reads

* Select two reads with the largest overlap

* Merge the two reads

* Repeat

The result is of course not guaranteed to be optimal as the greedy strategy may induce assembly errors, especially around repeated sequences (Figure 1.5). Popular Sanger assembler such as TIGR [17] or CAP3 [18] were greedy and broadly used because of their efficiency. This strategy was also reintroduced later to handle very short reads (around 25 to 50 bases) and implemented in assemblers like Ssake [19] and Vcake [20]. Those tools produce acceptable results on simple genomes. On more redundant genomes, such approaches produce too much assembly errors and other techniques are now favored. We can also criticize the model of the shortest common superstring, as in the presence
<figure><img src="image_8.png" /><figcaption>Figure 1.6: Toy example of the shortest common superstrings of a sequencing dataset. In this example, the genome contains a repeat in light and dark red. When the shortest common superstrings are computed we observe that none of the two SCS match with the genome because they are both shorter. In this example the SCS computation do not lead to the genome reconstruction.</figcaption></figure>

of repeated sequences, multiple SCS can exist and genome may not be one of the SCS (Figure 1.6).

### 1.2.2 Overlap Layout Consensus

The overlap layout consensus paradigm core notion is the overlap graph. This framework is the most general among the three paradigms as we can argue that all assemblers use an overlap graph implicitly. The objective is to know how all reads can be positioned in relation to each others, to represent those connections in a graph and to consider all overlaps (not only maximal ones) to produce a solution. We know how reads can be ordered by knowing how they overlap. The overlap graph is a graph where reads are nodes, connected if they overlap significantly (Figure 1.7). The algorithm can be outlined by:

* Overlap: calculate pairwise overlaps between reads

* Layout: look for a parsimonious solution (as a generalized Hamiltonian path visiting each node at least once while minimizing the total string length)

* Consensus: merging reads, using redundancy to correct sequencing errors

The first OLC assembler was Celera [21] and was designed to handle Sanger sequences. Celera uses a BLAST-like [22] approach to compare each read to the others and to find significant overlaps. Then it compacts the overlaps presenting no ambiguity (Figure 1.8) and tries to apply heuristics on the complex cases involving repeats. The final sequences
<figure><img src="image_9.png" /><figcaption>Figure 1.7: Overlap graph toy example where only overlaps of size 3 or more are considered.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 1.8: Example of non ambiguous compactions in green and unsafe compactions in red. The green compactions are the only possible choices, thus there is no ambiguity on which compaction should be performed. But the third read could be compacted with two other reads, the two compactions are indistinguishable. Choosing one compaction over the other could lead to assembly error.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 1.9: How repeated sequences can result in unsolvable graph. Two genomes A and B are sequenced, the reads covering their different regions are compacted into a simplified assembly graph of unambiguous sequences represented by blocks. Two repeats are present in red and purple. Several remarks can be made. First, the two genome A and B share the same simplified assembly graph structure despite being different. Secondly, several minimal solution may exist for a given assembly graph. Thirdly, because of repeats, the genome B will not be a minimal solution of the given assembly graph.</figcaption></figure>

are produced via a consensus to remove most sequencing errors. For comparison, in the toy example of Figure 1.5, the OLC approach would achieve to solve the assembly by considering even non maximal overlaps and to produce the correct path, composed of strictly less nucleotides than the greedy solution.

The fact is that, with either paradigm, a perfect assembly is in most cases impossible to obtain. Sometimes the information available is not sufficient to make sound choices. In those cases the parsimonious strategy of no choosing between two indistinguishable possibilities is applied (Figure 1.8). This results into fragmented assemblies constituted of consensus sequences that are supposed to be genome substrings. We call those sequences "contigs" for contiguous consensus sequence [23]. In the example of the Figure 1.9 an assembly graph is created from the reads information. The graph can be simplified by compacting reads that overlap unambiguously into contigs. Assembly can become complex in multiple ways. First, different assemblies can be proposed from this graph even considering only minimal solutions. The green and yellow contigs are interchangeable in the two minimal solutions of Figure 1.9. Secondly, different genomes can share very similar assembly graphs. In Figure 1.9, both genomes A and B would be represented by the same simplified assembly graph. Thirdly, sometimes the solution is not a minimal substring: the genome B is not generated as minimal solution of the assembly graph. The most parsimonious way is therefore to output the proposed contigs represented by the colored blocks. To give orders of magnitude of the fragmentation of a genome into contigs we can look at published assemblies. A *E. coli* genome has almost a hundred contigs https://www.ncbi.nlm.nih.gov/assembly/GCF_002099405.1, a *C. elegans* genome count more than 5,000 contigs https://www.ncbi.nlm.nih.gov/assembly/GCA_000939815.1/.

But even by applying parsimonious rules, assemblers may make mistakes and produce assembly errors (misassemblies). Some tools, such as QUAST [24], are designed
**Reference genome:**

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

**Misassembled contig:**

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

**Possible assembly graph:**

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

Figure 1.10: A classic error of assembly: the relocation. In the example the assembler produced a contig AC where the genome sequence is ABC. This error may be explained by different scenario, like a missing edge between A and B or a wrong choice of the assembler.

to evaluate the quality of an assembly by mapping the contigs on a reference genome and providing a classification and quantification of the different misassembly types. The classic misassembly is the relocation, when the assembler merges two sequences that are not consecutive in the genome into a contig. An example of relocation is shown on Figure 1.10: the assembler outputs the contig AC where the actual genome sequence is ABC. If the size of B is below a threshold (1 kilobase by default for QUAST) it is considered as a minor mistake as the contigs will still be mapped in their entirety on the reference genome and B will be considered as a deletion error. If B is larger this may be problematic as it induces a chimeric genome structure.

QUAST detects two other types of error (Figure 1.11).

* The translocation, where two (or more) parts of a contig come from different chromosomes

* The inversion, when a part of the contig is the reverse complement of the actual genome sequence

Such errors may appear for various reasons and may be due to the data or the assembler strategy. The connection AB may have been missed, leading the assembler to produce AC or a heuristic choice may have chosen the path AC over the path AB.

Users want assemblers both to produce long contigs and as few misassemblies as possible. The overlap graph is able to produce less errors than the greedy approaches by considering a more global information. But it may lead to very heavy structures with huge numbers of nodes and edges. In order to cope with this problem, a new improved model, the String graph [25], has been proposed. The essential conceptual difference between overlap graphs and string graphs is the transitive reduction of edges. When A overlaps B and C and B overlaps C, the edge A to C is removed because it can be "deduced" by transition from A to B and B to C (Figure 1.12). Later, [26] and [27] proposed further computational improvements to compute and represent a string graph.
<figure><img src="image_15.png" /><figcaption>Figure 1.11: Assembly error types defined by QUAST. From http://quast.bioinf.spbau.ru/manual.html.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 1.12: Transitive reduction of edges, the main difference between overlap graph and string graph. In the left part the connection between a and c can be deduced by the connections ab and bc and is removed in the right part.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 1.13: The previous overlap graph (Figure 1.7) toy example made into a string graph by removing transitively-inferrible edges.</figcaption></figure>

This paradigm was used a lot with long Sanger sequences and for relatively small genomes. Because of the cost of the pairwise overlaps computation, the OLC is too time consuming on high number of short reads from NGS. Thus, other solutions had to be found to be able to deal with the amount of reads to assemble large genomes.

### 1.2.3 De Bruijn graphs

**De Bruijn graph usage** The de Bruijn graph is a directed graph representing overlaps between sequences of symbols, named after Nicolass Govert de Bruijn [28]. Given an alphabet $\sigma$ of $m$ symbols, a $k$ dimensional de Bruijn graph has the following properties.

* $m^k$ vertices produced by all words of length $k$ from the alphabet $\sigma$

* Two vertices A and B are connected by an edge from A to B if and only if the $k-1$ suffix of A is equal to the $k-1$ prefix of B.

This graph has interesting properties and several applications in networking [29], hashing [30] and bioinformatics for genome assembly. Even if the graph used for assembly is called a de Bruijn graph, it is not exactly a de Bruijn graph as defined above.

The first application of the de Bruijn graph in genome assembly was introduced into the EULER assembler [31] in order to tackle assembly complexity. The idea was to consider a partial de Bruijn graph on the alphabet (A,C,T,G) constructed only with the vertices whose words of length $k$, called kmers, appeared in the sequencing data. The intuition of this approach is the following (Figure 1.14):

* A read is represented as a path in the graph

* Reads that overlap with more than $k$ nucleotides will share some kmers

* Extracting paths of such graph will produce assembled reads
<figure><img src="image_18.png" /><figcaption>Figure 1.14: Intuition of how the de Bruijn graph is able to merge overlapping reads. In this figure, the two read share an overlap of 5 bases. In the de Bruijn graph they share 2 kmer and are part of a path of the graph whose sequence is the concatenation of the two reads.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 1.15: Classic representation of canonical kmers and their different kinds of overlaps. For each kmer, the (lexicographically) smaller kmer between the sequence and its reverse complement is chosen as representative and is called the canonical kmer. Since a kmer and its reverse complement are indistinguishable, several types of k-1 overlap have to be considered in a de Bruijn graph. The first one is the forward forward (FF), that comes from canonical to canonical. The second one is the forward reverse (FR), from the canonical to the reverse. The third one is the reverse forward (RF) from the reverse to the canonical.</figcaption></figure>

Another technical point comes from the fact that the sequencing data is not stranded. As previously mentioned, since a read can be extracted from a strand or another, a sequence and its reverse-complement should be considered as the same object. In order to handle this, the principle of canonical kmer is used. A kmer *w* is called "canonical" if and only if it is smaller lexicographically than its reverse complement *v*. The kmers that are canonical are just inserted in the graph when the reverse complement of non canonical kmers are inserted instead of them. This modification leads to a change in the definition of how two kmers overlap and therefore are connected in the graph. In this graph four connections are possible (Figure 1.15):

* Suffix → prefix

* Suffix → reverse prefix

* Prefix → suffix

* Prefix → reverse suffix

This is rather a technical problem than an algorithmic one. Practical and theoretical results show that the difficulty of assembly is not due to this stranding problem [32]. For the sake of clarity we will ignore this property in most cases.
**De Bruijn graph and overlap graph** The de Bruijn graph theoretically achieves the same tasks than the overlap graph, while being conceptually simpler and much more efficient for the three reasons detailed in the following:

* No alignment

* Abstracted coverage

* No consensus

The Figure 1.14 shows how the de Bruijn graph finds (exact) overlaps of length superior to $k - 1$ between two reads. The de Bruijn graph does not explicitly compact reads together. However, selecting long paths from the de Bruijn graph is very similar to compacting overlapping reads in the OLC.

The de Bruijn graph became widely used when the short reads from NGS appeared, as it was better suited than the OLC to handle this kind of sequencing data. The OLC approach did not scale well on the high number of sequences generated by NGS. The use of the de Bruijn graph is very interesting for short read assembly for its ability to deal with the high redundancy of such sequencing in a very efficient way. Indeed a kmer presents dozens of time in the sequencing dataset appears only once in the graph. This makes the de Bruijn graph structure not very sensible to the high coverage, unlike the OLC. The de Bruijn graph was first proposed as an alternative structure [31] because it was less sensible to repeats. Repeats that were problematic in the OLC, creating very complex and edges heavy zones, are collapsed in the de Bruijn graph.

**Redundancy** The high redundancy in the sequencing data can also be used to efficiently filter the sequencing errors. The first highly used de Bruijn graph short reads assembler was Velvet [33]. It introduced core notions of de Bruijn graph assembly, as the idea of error filtering based on kmer abundances. With a high coverage, we expect a high abundance for most kmers. A kmer with very low amount of occurrences is therefore very likely not a genomic kmer (a kmer present in the genome) but rather an erroneous one (not present in the genome and coming from a sequencing error). By admitting in the de Bruijn graph only kmers with a coverage above a threshold, we can get rid of most erroneous kmers almost without losing genomic kmers, but the ones that has an unexpectedly low abundance. Kmers whose abundance are over the abundance threshold (or solidity threshold) are called "solid". In Figure 1.16 we can see that sequencing errors generate a huge number of low abundance kmers.

**De Bruijn graph patterns** De Bruijn graph assembly basically consists in graph simplification by applying heuristics on knowns patterns (Figure 1.17). They rely on path exploration, applying different strategies to handle recurring motifs. After graph simplification, they output the long simple paths as their contigs (Figure 1.18).

The simplest pattern is the tip (or dead end), a short path in the graph that is not extensible because its last kmer has no successor. If the tip is short, shorter than the read length for example, then it is very likely to be due to a sequencing error. But if it is large it could just be the beginning or an end of a chromosome. A basic assembly step is to recognize those short tips and remove them from the graph to simplify it and to allow longer contigs.

Another frequent pattern is the bubble. A bubble arises when several paths start from a kmer and end in another kmer. A sequencing error can create a bubble if the
<figure><img src="image_20.png" /></figure>

Figure 1.16: Comparison of kmer spectra with and without sequencing errors. The two kmer spectra are computed from simulated reads with 1% and 0% errors respectively from *E. coli* with a coverage of 100X. We observe two similar curves but the green one presents a huge quantity of low abundance kmers due to the presence of sequencing errors.

sequencing error is positioned *k* nucleotides apart from the end and the start of a read.

Another source of bubbles can be heterozygosity. For organisms with (at least) a pair of homologous chromosomes, heterozygosity defines the scenario of owning two different versions of a sequence (called alleles) on each chromosome. A well known example is blood type:, for instance an AB individual is heterozygous. A diploid (polyploid) genome is a genome containing two (multiple) complete sets of chromosomes. Within a diploid (polyploid) genome, heterozygosity can happen at various positions. When a diploid (polyploid) genome is sequenced, two (multiple) close sequences are simultaneously sequenced and the original allele of each read is unknown. When the two sequences are almost identical with some minor differences, this creates bubbles in the de Bruijn graph.

A last important source of bubbles are the quasi-repeats. When two almost identical sequences appear in the genome, it can create the X pattern specific of the repeats (Figure 1.19), with bubbles inside the central repeated sequence (Figure 1.20).

Those patterns highly depend on the size of *k*. The *k* parameter, called the "order" of the de Bruijn graph, is a key factor in de Bruijn graph assembly. For a given dataset, distinct orders can lead to the construction of extremely dissimilar graphs.

**Repeats** The first point about the size of *k* is that the smallest *k* is, the more complex the graph will become. With a small *k*, the probability of a kmer (or a *k* − 1mer) to be present at multiple positions in the genome is high. If this happens, then the two occurrences of the repeated kmer are collapsed in the graph and create a "X" structure (Figure 1.19). When such a motif is encountered, the assembler has no mean *a priori* to know how to continue and stop. A *k* too small may results into a profusion of such scheme and produce a fragmented assembly. In a similar way, a repeat of size *k* − 1 can also create an edge between unrelated kmers called "spurious edge". In other words, the smaller the *k* value is, the less significant the connections between kmers are. To see how the size of *k* affects the graph complexity, we present in Table 1.2 the number of repeated kmers in the de Bruijn graph created from reference genomes according to the
<figure><img src="image_21.png" /></figure>

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

<figure><img src="image_23.png" /><figcaption>Figure 1.17: Patterns in a de Bruijn graph with $k = 4$. The first pattern (A) is a tip generated by a sequencing error at the end of a read. The second one (B) is a bubble generated by a sequencing error in the middle of a read. The third one (C) is a bubble generated by an actual variation in the sequenced genome.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 1.18: Example of contigs generation. In this toy example, the tip in red is removed, the bubble is crushed and the green path is chosen over the orange one. But the assembler is not able to choose between the two gray paths to extend its contig so it stopped. The two gray paths are output as contigs with the large blue one.</figcaption></figure>

k=4

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

k=6

<figure><img src="image_26.png" /><figcaption>Figure 1.19: Example of spurious edge and X pattern due to repeats of size k and *k* − 1. With *k*=4 we got a repeat of size 4 that creates the X pattern framed in green. A repeat of size 3 also create a spurious edge between GCAT and CATA. With *k*=6 the graph is linear since no repeat of size 5 exist.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 1.20: Example of quasi-repeats creating a bubble. The two regions R and R' are almost identical, but the small difference between the two occurrences create a bubble inside the X pattern.</figcaption></figure>

<table><caption>Table 1.2: Number of kmers that appear multiple times in the reference genome according to the size of *k*.</caption><thead><tr><th>Reference genome</th><th>kmer size</th><th>Repeated kmer</th></tr></thead><tbody><tr><td><i>E. coli</i></td><td>11</td><td>976,822</td></tr><tr><td><i>E. coli</i></td><td>15</td><td>104,501</td></tr><tr><td><i>E. coli</i></td><td>21</td><td>33,745</td></tr><tr><td><i>E. coli</i></td><td>31</td><td>30,273</td></tr><tr><td><i>E. coli</i></td><td>61</td><td>25,575</td></tr><tr><td><i>E. coli</i></td><td>101</td><td>22,345</td></tr><tr><td><i>E. coli</i></td><td>201</td><td>17,985</td></tr><tr><td><i>E. coli</i></td><td>301</td><td>15,421</td></tr><tr><td><i>C. elegans</i></td><td>15</td><td>14,299,107</td></tr><tr><td><i>C. elegans</i></td><td>21</td><td>3,290,873</td></tr><tr><td><i>C. elegans</i></td><td>31</td><td>2,475,913</td></tr><tr><td><i>C. elegans</i></td><td>61</td><td>1,892,518</td></tr><tr><td><i>C. elegans</i></td><td>101</td><td>1,524,266</td></tr><tr><td><i>C. elegans</i></td><td>201</td><td>1,052,332</td></tr><tr><td><i>C. elegans</i></td><td>301</td><td>824,254</td></tr><tr><td>Human</td><td>15</td><td>250,505,977</td></tr><tr><td>Human</td><td>21</td><td>148,690,202</td></tr><tr><td>Human</td><td>31</td><td>122,846,758</td></tr><tr><td>Human</td><td>61</td><td>102,800,268</td></tr><tr><td>Human</td><td>101</td><td>92,294,300</td></tr><tr><td>Human</td><td>201</td><td>85,392,911</td></tr><tr><td>Human</td><td>301</td><td>81,190,289</td></tr></tbody></table>

size of *k*.

We observe that a low *k* value creates a de Bruijn graph with many repeated kmers that will create X patterns and therefore a complex graph. This can be explained by
<table><caption>Table 1.3: Proportion of kmer that appear a single time in the reference genome according to the size of *k*.</caption><thead><tr><th>Reference genome</th><th>kmer size</th><th>Percent unique kmer</th></tr></thead><tbody><tr><td><i>E. coli</i></td><td>9</td><td>1.1</td></tr><tr><td><i>E. coli</i></td><td>11</td><td>33.2</td></tr><tr><td><i>E. coli</i></td><td>15</td><td>97.7</td></tr><tr><td><i>E. coli</i></td><td>21</td><td>99.3</td></tr><tr><td><i>C. elegans</i></td><td>11</td><td>0.6</td></tr><tr><td><i>C. elegans</i></td><td>15</td><td>76.4</td></tr><tr><td><i>C. elegans</i></td><td>21</td><td>96.4</td></tr><tr><td><i>C. elegans</i></td><td>31</td><td>97.4</td></tr><tr><td>Human</td><td>15</td><td>25</td></tr><tr><td>Human</td><td>21</td><td>93.4</td></tr><tr><td>Human</td><td>31</td><td>95.1</td></tr></tbody></table>

the fact that small kmers are expected to appear multiple times even on a random word. There are 262,144 different 9mers (without accounting for reverse-complements) while the *E. coli* genome is a word of almost 5 millions nucleotides. We expect most 9mers to be present multiple times in the genome. The Table 1.3 confirms this fact. With *k* = 11 (4,194,304 11mers) most kmers are still repeated but with *k* = 15 and higher (1,073,741,824 15mers), most kmers are unique (present only one time in the genome). In practice we need to use a *k* such that $4^k$ is order of magnitude larger than the genome size in order to get a de Bruijn graph with mostly unique kmers. A higher value reduces the number of repeated kmer, with a decreasing efficiency because of the existence of large repeats in genomes. *C. elegans* reference genome presents almost a million repeats larger than 300 nucleotides.

**Overlap detection** The second point about kmer size is that a de Bruijn graph will only (implicitly) detect overlaps larger than *k* − 1 between reads. A high *k* value imposes that reads have to share large overlaps, and a high redundancy may be necessary to do so. A too large *k* value will create holes in the de Bruijn graph (Figure 1.21).

**Error removal** The third point is about sequencing errors. A read of size *L* has *L*−*k*+1 kmers. Thus a high *k* value means less kmers per reads. For example a read of length 100 has 70 31mers or 40 61mers. With a higher *k* value, the abundance of all kmers is lower and it is more difficult to differentiate genomic from erroneous kmers, based on their abundance. Another problem is that sequencing errors "destroy" kmers (Figure 1.22). When *L* ≤ 2 * k* − 1 a sequencing error can make erroneous all kmers from a read. The Table 1.4 shows the effect of the errors and the kmer size on the presence/absence of the genomic kmers without any kmer filtering. 6,716 missed kmers can seem like a low rate, around 1/1000 of the genomic kmers. But each of these "holes" fragments the assembly as seen in Figure 1.21. We observe that reads with a reduced error rate (substitutions) allows for the usage of a higher kmer size with the same coverage.

Given those effects summarized in Table 1.5, the size of *k* should be chosen according
<figure><img src="image_28.png" /></figure>

Figure 1.21: How the use of high values of *k* and solidity thresholds can fragment the de Bruijn graph. The red lines represent the paths of the de Bruijn graph according to the *k* value and the filter usage. The filter consists in removing the unique kmers. We can see that using a high *k* value can fragment the assembly because no large enough overlap exists. Unique kmers removal can also create holes when removing errorless kmers. Using a too large *k* size or a too high solidity threshold creates holes in the de Bruijn graph and results into a fragmented assembly.

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

Figure 1.22: How sequencing errors can produce many erroneous kmers. Red positions are sequencing errors and kmer in red are therefore erroneous kmers while green kmers are genomic. In the left figure, *k* = 8 thus the read has 15 kmers. Because of the sequencing error, 6 of them are erroneous. In the left figure, *k* = 12 thus the read has only 9 kmers. Because of the sequencing error 8 of them are erroneous. If the sequencing error was one base to the left, all kmer would have been erroneous.
<table><caption>Table 1.4: Number of genomic kmers not present in a set of 100 base pairs simulated reads with a coverage of 30X from *E. coli*, according to the error rate.</caption><thead><tr><td>Error rate</td><td>kmer size</td><td>Missed kmers</td></tr></thead><tbody><tr><td>1%</td><td>31</td><td>3</td></tr><tr><td>1%</td><td>61</td><td>6,716</td></tr><tr><td>1%</td><td>91</td><td>1,374,153</td></tr><tr><td>0.1%</td><td>31</td><td>2</td></tr><tr><td>0.1%</td><td>61</td><td>53</td></tr><tr><td>0.1%</td><td>91</td><td>293,675</td></tr></tbody></table>

<table><caption>Table 1.5: Effect of *k* parameter on a de Bruijn graph.</caption><thead><tr><th></th><th>Low k</th><th>High k</th></tr></thead><tbody><tr><th>Benefits</th><td>Detect even small overlap between reads kmers are more covered (more kmer by reads)</td><td>Less repeats in the graph Less spurious edges</td></tr><tr><th>Drawbacks</th><td>More repeats appear in the graph More spurious edges</td><td>Only large overlap are detected kmers are less covered</td></tr></tbody></table>

mainly to the coverage of the sequencing dataset and to its error rate. Some methods concentrate on automated selection of this crucial factor, as kmerGenie [34] in order to produce the best assembly possible given a read set.

**Assembly quality** Ideally we would want to use the highest possible kmer size in order to reduce the number of repeated kmers and spurious edges. To do so without having a graph abounding with holes, multiple solutions exist:

* Get a huge coverage

* Correct sequencing errors

* Use several kmers size

The first solution has no real drawback, but is not satisfying because of the supplementary expenditure. Even if sequencing costs are lowering, sequencings still represent important investments. In some cases the bottleneck may be the quantity of DNA material accessible, that would make impossible to get more sequencing.

The second solution, used by many assemblers [33, 35, 36, 37, 38], is to use a preprocessing to correct sequencing errors. Those techniques use the redundancy of the read set to correct errors based on statistical methods. It presents the interest to use reads with less errors and therefore containing longer errorless sequences. Read correction step can be criticized as it may modify correct bases or may correct errors incorrectly.

The third solution, used in state of the art assembly approaches (called multi-k approaches), is to make use of multiple kmer sizes. The idea of those approaches is to get the graph connectivity from the low kmer sizes, and the low amount of repeats from the high kmer sizes. The precise method may vary among tools but it generally follows the
<figure><img src="image_30.png" /><figcaption>Figure 1.23: Solving a repeat consists in the determination of the possible contexts around the repeat and therefore limits the fragmentation. To solve the presented toy repeat, the assembler must link the green parts and the orange parts together.</figcaption></figure>

pattern described by IDBA [39] in order to construct the best possible de Bruijn graph:

1. Start with a low $k$

2. Create a de Bruijn graph from the reads with $k$

3. Produce contigs using this graph

4. Raise the value of $k$

5. Create a de Bruijn graph from the reads and the contigs with $k$

6. Repeat from step 3 until $k \le k_{max}$

This method allows the use of large $k$ values even if the coverage is not sufficient everywhere. The assembly step with a low $k$ value detects the possible compactions between reads that share low overlaps. The large kmer size allows to "solve" some repeats as presented in Figures 1.23 and 1.19. Those approaches are still not perfect since complex or poorly covered regions can be absent from the contigs. Another problem of those techniques is that they may be computationally challenging since they perform several assembly steps, in a sequential way (time consuming) or simultaneously (memory consuming).

A plethora of de Bruijn graph based assemblers have been developed, ABYSS [40], SOAPdenovo [41], IDBA [39], ALLPATHSLG [42], gossamer [43], minia [44], SPADES [36]. Each brings to the table extremely different properties and trade-off between resources needed, efficiency and assembly quality.

But even the best combinations of those strategies are not able to use the full read length. Despite being conceptually simpler than OLC, the de Bruijn graph is inferior in terms of information usage because $k \le \text{read size}$. Even if we use a very large $k$ value, we may encounter repeats with $k < \text{repeats size} < \text{read size}$ that could be solved by a string graph but not by the de Bruijn graph. The multi-k de Bruijn graphs tend to
come close to the string graph (various overlap sizes detection, usage of a large part of the reads) but also tend to lose their performances advantages.

### "Genome assembly" core messages:

* Assemblers goal is to order the reads and concatenate them into larger sequences while removing the sequencing errors

* In most cases, the assembler is not able to completely recover the genome and will output "contigs" that are sequences supposed to be substring of the genome larger than the reads

* Greedy approaches are efficient but are not suited to large or complex genomes as they handle repeats very poorly, resulting into erroneous assembly.

* Overlap graph is a framework where the overlaps between the reads are computed and put in a graph where reads are nodes, connected if they overlap.

* The de Bruijn graph is a more restricted framework where the words of length *k* from the reads are nodes, connected if they share an overlap of *k* − 1.

* OLC and de Bruijn graph approaches select paths from their graphs using heuristics in order to produce contigs based on the graph structure

* The overlap graph approaches relies on heavy data structures while the de Bruijn graph, conceptually simpler, scales better on large datasets

## 1.3 Assembly hardness

Assemblers are supposed to produce the largest possible contigs while making as few assembly errors as possible from the available read information. Here we describe the challenges that genome assembly may present and the limitations of existing approaches.

### 1.3.1 Repeats

As we have seen, one of the problems of assembly is the repeated sequences. In fact repeats may be the fundamental problem in genome assembly. An analysis of the complexity of the assembly problem according to the type of data available has been made in [45]. The principal conclusion of the paper is that repeats larger than reads are impossible to solve. In fact, for most genomes, the perfect assembly is not achievable using only NGS reads. It may seem odd that such large repeats exist. The fact is that genomes are absolutely not random sequences. When we compare Table 1.6 with Table 1.2, we see that there is no large repeated kmers in random sequences whereas *E. coli* genome has thousands of repeated 301mers. The conclusion is that NGS reads sequences are not enough to produce a complete assembly. [46] has shown the effort necessary to finish short reads assembly. Therefore longer range information is required in order to solve the large repeats and produce more continuous assembly.
<table><caption>Table 1.6: Number of repeats of size *k* in a random genome of the same size than *E. coli* or *C. elegans*.</caption><thead><tr><th>Reference genome</th><th>kmer size</th><th>Repeated kmer</th></tr></thead><tbody><tr><td>Random Genome (<i>E. coli</i>)</td><td>15</td><td>84,148</td></tr><tr><td>Random Genome (<i>E. coli</i>)</td><td>21</td><td>28</td></tr><tr><td>Random Genome (<i>E. coli</i>)</td><td>31</td><td>0</td></tr><tr><td>Random Genome (<i>C. elegans</i>)</td><td>15</td><td>34,050,959</td></tr><tr><td>Random Genome (<i>C. elegans</i>)</td><td>21</td><td>17,953</td></tr><tr><td>Random Genome (<i>C. elegans</i>)</td><td>31</td><td>0</td></tr></tbody></table>

<figure><img src="image_31.png" /><figcaption>Figure 1.24: Differences between Paired End sequencing and Mate Pairs sequencing. In Paired End, small DNA fragment are selected and sequenced from both end. For Mate Pairs, longer fragment are selected, both ends are marked and connected into a circular sequence. The sequence is then cut and the marked part is sequenced from both and generates an opposite sequencing orientation.</figcaption></figure>
1.3.2 Scaffolding

As denoted before, most assemblers use the reads sequences to produce contigs, that are a set of fragments of the genome. In order to improve such genome drafts, other informations and techniques can be used. It is possible via sequencing techniques to obtain short reads that are associated by pairs coming from related positions of the genome. There are two main sequencing techniques called "Paired End" and "Mate Pairs" sequencing (Figure 1.24). In both cases, we get a pair a short reads with the same characteristics than short reads described previously. The additional information is that we have an estimation of the distance covered by the read pair, called the fragment size, since we know that both come from the same DNA fragment. For example, a paired end sequencing of 2*250 base pairs with an fragment size of 800 will produce pairs of reads of 250 bases pair reads spaced by 300 nucleotides in the genome (on average). The paired end sequencing can handle fragment size up to 2,000 bases while mate pair can produce fragment size from 2 to 20 kilo-bases. Another difference is that in Paired End sequencing, the first fragment is read forward and the second is reversed while in Mate Pairs this is the opposite, as shown in Figure 1.24. This can be explained by the difference between the two protocols. In paired end sequencing, fragments of the desired size are selected and sequenced from both ends. In mate pairs sequencing long fragments end are marked, circularized, fragmented, and sequences with marker are sequenced from both ends.

The standard way to use those linking information is to align paired reads on the contigs and use the pair of reads that map on different contigs to order and orient them with respect to each other. This task is frequently called scaffolding, because contigs are arranged together into "scaffolds", based on estimations of the distance between them (Figure 1.25). A complete assembly workflow is summarized in Figure 1.26.

[47] formally defined the problem and showed its NP-completeness and therefore its potential intractability. Most scaffolders try to optimize the number of satisfied links in order to produce scaffolds. The main problem of scaffolds is that they may contain "holes". Scaffolds are basically ordered contigs. When the different contigs used are not overlapping, the scaffolder usually estimates the distance between them and fills the "holes" with 'N' characters to inform that the sequences here are unknown. The scaffolders may be able to order large contigs, but can fail to find the contigs separating them. This may be the case if those contigs are too small or simply if they were not output by the contig generation step. Since the scaffold approach is based on mapping the reads on the contigs, very short contigs (shorter than reads) cannot be considered this way because of multiple mapping problems (when a read maps on several contigs).

Recent scaffolders like SSPACE [48] try to extend contigs using the pairs information earlier to connect them in order to reduce the potential gaps in the later phase. Some standalone tools called gap-giller, as MindTheGap [49], are specifically designed to fill such holes. Other kinds of information start to be used for scaffolding such as long reads from third generation sequencing [50] and 3C [51] a sequencing technology liking reads coming from the same chromosome.
<figure><img src="image_32.png" /><figcaption>Figure 1.25: How the paired reads information can be used to solve some complex scenario. Paired reads show that blue and yellow contigs are separated by one contig, as yellow and red and green and pink. Using this information, we can eliminate solutions where those distances are not correct. The first solution does not respect any links (Blue and yellow are too far apart as green and pink and yellow and read are too close) and is therefore rejected. In this example the solution respecting the link is equal to the genome.</figcaption></figure>

### 1.3.3 Multiple genomes

The assembly problem is defined as a sequence reconstruction. In many cases we are interested in assembling multiple genomes at once. We briefly introduced heterozygosity and highlighted that for many species, as human, two sets of chromosomes, or more, are sequenced. We will refer to a complete set of chromosomes coming from an individual as an haplotype. A human has therefore two haplotypes because it owns two versions of each chromosome (but the sexual ones). Intuitively, we are sequencing two similar genomes at once. We can define the heterozygosity rate by the ratio of the edit distance between the two alleles over the size of the region. Regions identical across haplotypes will be called homozygous while regions with variations will be called heterozygous. Since for many model species such as the human the heterozygosity rate is very low (around 0.1% for the human), most assemblers do not try to assemble both sequences and rather ignore the heterozygosity information. They usually produce a sequence presenting a mix of the different haplotypes as it is shown in Figure 1.27. In an assembly graph such variants may appear as bubbles and most assemblers will select a path through the bubble without trying to conserve haplotype information. Some tools try to assemble separately the different haplotypes when possible as Platanus [52]. In such cases, all homozygous regions can be considered as repeats as they appear identically in two different chromosomes. But other scenarios can force users to assemble multiple genomes at once.

Some species cannot be grown and sequenced in wet labs, therefore their whole environment (sea water for instance) is picked and sequenced, then we must assemble the
<figure><img src="image_33.png" /><figcaption>Figure 1.26: Summary of assembly steps. The global picture of an assembly of *E. coli* from simulated reads with SPAdes. Some refer to assembly for contigs generation and scaffolding steps while other call the contigs generation step as assembly and separate it from the scaffolding.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 1.27: Example on how an assembler may crush haplotypes. The assembler detects two bubbles and crushes them, choosing one path over another. The resulting contigs may contain sequences from both haplotype mixed.</figcaption></figure>
whole data (meta-genomic). Some assemblers are designed to assemble meta-genome as
meta-IDBA [53] or meta-velvet [54]. As their names suggest, they are usually based on
a regular genome assembler adapted to fit to the proposed scenario. In some cases, a
single species is sequenced but with several individuals having different, however closely
related, genomes (pool-seq). Since repeats assess the toughness of assembly, we can sort
scenarii by their apparent hardness:

* Single haploid genomes (some repeats)

* Heterozygous genomes (each homozygous regions is a repeat)

* Pool-seq, multiple related genomes (genome sized repeat)

* Meta-genomic, multiple non-related genomes (genome size repeats, some repeat shared among species)

It seems obvious that perfect assembly and distinguishing closely related genomes are
currently out of reach. But this gives objectives to be met in the assembly field.

"Assembly challenges" core messages:

* Finished assembly is often impossible because of repeats longer than the reads

* Other kinds of data may be used to order contigs in larger sequences called "scaffolds"

* Polyploid, meta-genome and pool-seq assembly are even harder than regular assembly because of systematic repeats

1.4 Outline of the thesis

In this introduction, we highlighted three "challenges" of the assembly process of NGS
data:

* The high amount of resources required

* The fragmentation of produced assembly

* The hardness to assemble complex (large and/or repeat-rich) genomes

To address those issues, we will present new approaches based on the de Bruijn graph.
In the first chapter we present new ways to represent and construct the de Bruijn graph
in order to make it more scalable and more efficient. In the second chapter we introduce
new methods to use the de Bruijn graph as a reference and show the interest of such
structures over a set of fragmented contigs. In the third chapter we propose techniques
to construct high order de Bruijn graphs, allowing the use of the reads information, as
the string graph does, without losing the de Bruijn graph efficiency.
# Chapter 2

## Handling assembly
In this chapter we will assess the computational aspect of assembly. We show why the scalability of methods used for genome assembly can be critical (Section 1). We provide an overview of the state of the art of efficient methods and structures to address such problems for both overlap graph (Section 2) and de Bruijn graph (Section 3) frameworks. Then we present our theoretical and practical contributions. We show that de Bruijn graph assembly may be done using a Navigational Data Structure, a novel model that we introduce, which shows lower theoretical memory bounds. We also propose a new proof of concept assembler that shows practical resources improvement over state of the art tools (Section 4). Those implementations are based on the idea to enumerate and index simple paths of the graph. We argue that such enumeration is a bottleneck in many assembly methods and we provide resources efficient methods to answer this need (Section 5). To do so, we make use of minimal perfect hashing functions as very efficient indexes and provide a new method to compute such functions on very large sets of keys (Section 6).

## 2.1 The assembly burden

In the last section we described the sequencing data but we gave no clue about the typical size of a genome. The fact is that this size can vary a lot among living organisms (Figure 2.2). Though we can indicate some orders of magnitude on genomes size:

* Virus: Thousands base pairs

* Bacteria: Millions base pairs

* Mammals: Billions base pairs

Some species, yet to be sequenced, are expected to present order of magnitude larger genome (Paris japonica, Polychaos dubium) [55] with presumed hundred of billions base pairs. The struggle for assembling such genomes may not be seen directly. Since high coverages are needed for assembly, sequencing datasets represent huge amount of information. Most assemblies typically rely on a mean coverage ranging from 30X to 100X and higher. Consequently, assemblers have to handle millions of reads, even for bacterial genomes. Larger genomes can count hundreds of millions base pairs for most prokaryotes, and up to billions for some mammals or plants genomes. Such sets can reach billions of reads and represent hundreds of gigabytes or even terabytes of data.

Dealing with this tremendous amount of information require either to use huge computational resources or to conceive specific algorithms and data structures designed for resource efficiency. As sequencing cost continues to decrease, sequencing very large genomes becomes affordable, but assembly of such genome is barely possible. If the running time is an important concern, it is usually not the source of intractability. Most assemblers rely on very large graph structure and index, that can require terabytes of RAM on large datasets. Most of the time, such memory requirement is more concerning than running time, as large scale clusters may not be easily accessible. Besides facing the challenge to produce correct assemblies, the future assemblers will have to handle larger and larger datasets, in order to deal with large genomes or meta-genomes while providing a high throughput to follow the sequencing rate. As sequencing costs are dropping (Figure 2.1), the computational resources necessary to treat them is becoming the financial bottleneck.
<figure><img src="image_35.png" /><figcaption>Cost to sequence a human genome (USD)</figcaption></figure>

Figure 2.1: Evolution of the cost of human sequencing. From https://upload.wikimedia.org/wikipedia/commons/e/e7/Historic_cost_of_sequencing_a_human_genome.svg.
<table><thead><tr><th>organism</th><th>genome size (base pairs)</th></tr></thead><tbody><tr><td><strong>model organisms</strong></td><td></td></tr><tr><td>model bacteria <i>E. coli</i></td><td>4.6 Mbp</td></tr><tr><td>budding yeast <i>S. cerevisiae</i></td><td>12 Mbp</td></tr><tr><td>fission yeast <i>S. pombe</i></td><td>13 Mbp</td></tr><tr><td>amoeba <i>D. discoideum</i></td><td>34 Mbp</td></tr><tr><td>nematode <i>C. elegans</i></td><td>100 Mbp</td></tr><tr><td>fruit fly <i>D. melanogaster</i></td><td>140 Mbp</td></tr><tr><td>model plant <i>A. thaliana</i></td><td>140 Mbp</td></tr><tr><td>moss <i>P. patens</i></td><td>510 Mbp</td></tr><tr><td>mouse <i>M. musculus</i></td><td>2.8 Gbp</td></tr><tr><td>human <i>H. sapiens</i></td><td>3.2 Gbp</td></tr><tr><td><strong>viruses</strong></td><td></td></tr><tr><td>hepatitis D virus (smallest known animal RNA virus)</td><td>1.7 Kb</td></tr><tr><td>HIV-1</td><td>9.7 kbp</td></tr><tr><td>influenza A</td><td>14 kbp</td></tr><tr><td>bacteriophage λ</td><td>49 kbp</td></tr><tr><td>Pandoravirus <i>salinus</i> (largest known viral genome)</td><td>2.8 Mbp</td></tr><tr><td><strong>organelles</strong></td><td></td></tr><tr><td>mitochondria - <i>H. sapiens</i></td><td>16.8 kbp</td></tr><tr><td>mitochondria - <i>S. cerevisiae</i></td><td>86 kbp</td></tr><tr><td>chloroplast - <i>A. thaliana</i></td><td>150 kbp</td></tr><tr><td><strong>bacteria</strong></td><td></td></tr><tr><td><i>C. ruddii</i> (smallest genome of an endosymbiont bacteria)</td><td>160 kbp</td></tr><tr><td><i>M. genitalium</i> (smallest genome of a free living bacteria)</td><td>580 kbp</td></tr><tr><td><i>H. pylori</i></td><td>1.7 Mbp</td></tr><tr><td>Cyanobacteria <i>S. elongatus</i></td><td>2.7 Mbp</td></tr><tr><td>methicillin-resistant <i>S. aureus</i> (MRSA)</td><td>2.9 Mbp</td></tr><tr><td><i>B. subtilis</i></td><td>4.3 Mbp</td></tr><tr><td><i>S. cellulosum</i> (largest known bacterial genome)</td><td>13 Mbp</td></tr><tr><td><strong>archaea</strong></td><td></td></tr><tr><td>Nanoarchaeum equitans (smallest parasitic archaeal genome)</td><td>490 kbp</td></tr><tr><td>Thermoplasma acidophilum (flourishes in pH<1)</td><td>1.6 Mbp</td></tr><tr><td>Methanocaldococcus (Methanococcus) jannaschii (from ocean bottom hydrothermal vents; pressure >200 atm)</td><td>1.7 Mbp</td></tr><tr><td>Pyrococcus furiosus (optimal temp 100°C)</td><td>1.9 Mbp</td></tr><tr><td><strong>eukaryotes - multicellular</strong></td><td></td></tr><tr><td>pufferfish <i>Fugu rubripes</i> (smallest known vertebrate genome)</td><td>400 Mbp</td></tr><tr><td>poplar <i>P. trichocarpa</i> (first tree genome sequenced)</td><td>500 Mbp</td></tr><tr><td>corn <i>Z. mays</i></td><td>2.3 Gbp</td></tr><tr><td>dog <i>C. familiaris</i></td><td>2.4 Gbp</td></tr><tr><td>chimpanzee <i>P. troglodytes</i></td><td>3.3 Gbp</td></tr><tr><td>wheat <i>T. aestivum</i> (hexaploid)</td><td>16.8 Gbp</td></tr><tr><td>marbled lungfish <i>P. aethiopicus</i> (largest known animal genome)</td><td>130 Gbp</td></tr><tr><td>herb plant <i>Paris japonica</i> (largest known genome)</td><td>150 Gbp</td></tr></tbody></table>

Figure 2.2: Genome size for a variety of selected organisms. From http://book.bionumbers.org/how-big-are-genomes/.
This fast evolving sequencing environment pushes assembly designers to conceive increas-
ingly efficient dedicated methods and data structures. In this chapter, we will present a
state of the art of such techniques with a focus on the scalability of the de Bruijn graph
assembly.

"The assembly burden" core messages:

* Genomes may count billions of bases

* Sequencing datasets can be composed of billions of reads and represent terabytes of data

* Genome assembly is a High Performance Computing issue

* Efficient and dedicated structures are needed especially to reduce the memory usage

* Decreasing sequencing costs raise the need to efficient assembly

2.2 Overlap graph scalability

The problem of assembly scalability appeared with the overlap graphs and large genomes.
Tools such as TIGR [17] or CAP3 [18] can be considered as brute force approaches and
were suitable only for small genomes. When the concern about repeat misassemblies
arose, the overlap graph approaches have become preferred. The core operation of the
overlap graph is to find all overlaps between reads from pairwise alignment. Therefore
standard overlap graph assembly methods have a worst time complexity quadratic with
the number of reads. This fact may be moderated as only significant overlaps are consid-
ered. Heuristics are used in order to quickly find significant alignment. Tools mainly use
BLAST [22], or other aligners following the "seed and extend" paradigm. The principle
of the seed and extend paradigm is to look for words in common between the query and
the reference (called seeds or anchors) and to try to extend the alignment from the shared
words. Unlike Smith-Waterman algorithm [56], such approaches are not guaranteed to
find all optimal alignments. But they are dramatically faster while providing almost iden-
tical results. Using such anchoring heuristics implies that finding the overlaps of a read
is roughly linear in the number of similar reads (according to the method and parameter
used) instead of being linear in the total number of reads.

However, with high coverage, one read presents a high amount of high quality overlaps
(Figure 2.3). This will impact the running time of each query trying to find the overlaps
shared by a read. Such alignment based approaches were quite efficient with Sanger
sequencing, since the number of reads remained relatively low. With a very low error
rate, low coverage is sufficient to complete correct assemblies, around 10X or 20X. Since
the sequences are also longer (around thousands of base pairs), the number of Sanger reads
is lower than the number of NGS reads for a given coverage. In practice the number of
reads, even for large genome assembly, was not above the million. The number of overlaps
needed to be stored in order to work on the graph could still be challenging. Furthermore,
repeated regions resulted in dense zones (called hub) with very high numbers of overlaps.
Celera [21] and Arachne [57] improved the running time by masking repeated regions to
<figure><img src="image_36.png" /><figcaption>Figure 2.3: How the coverage affects the number of overlaps in an overlap graph. On this toy example, with a 3X uniform coverage, a read shares an overlap superior to 4 bases with 2 other reads in green. With a coverage of 10X each read shares a significant overlap with 6 other reads. The reads in red are overlapping but not detected because of a small overlap with the query read.</figcaption></figure>
simplify the overlap graph. Later, the string graph [25] was proposed to abstract the redundancy present in the sequencing data, as does the de Bruijn graph, using transitive reduction. The assembler BOA [25] was designed to assemble Sanger reads with such string graphs. However with NGS technologies, deeper sequencing was affordable and helped to cope with the error rate. This raised up the number of reads by orders of magnitude, themselves leading to orders of magnitude higher numbers of overlaps to handle. The overlap graph and string graph need enormous resources to process such datasets. Thus de Bruijn graphs were preferred to handle this kind of data.

However, because of the fact that the string graph should theoretically lead to better assembly by using the whole reads instead of kmers, several contributions were proposed in order to allow the string graph to scale up with the tremendous number of short reads. The introduction of the very popular FM-index structure [58] in bioinformatics allowed memory efficient indexing of large DNA sequences. It was first used to index references and enabled very efficient read mapping (fast and stringent local alignment of short reads on a reference). Notable examples of such tools based on a FM-index are Bowtie [59] and BWA [60]. [26] used the same idea to improve the string graph performances. The rationale is to index the reads with a FM-index and to use the research operation to detect overlaps between sequences. The memory usage of the assembly of a human genome has been extrapolated by the author to 700GB, which is not intractable but still very expensive. The time required for assembly is also tractable, days for 100 mega-bases genome, but this is yet very slow compared to most de Bruijn graph assemblers. A more efficient implementation called SGA [27] proposed a distributed FM-index construction. Using multiple small indexes allows for the construction of a global FM-index using less than 60 GB of RAM for a human sequencing dataset of more than 1.2 billions reads. The efficient use of the FM-index allowed the string graph to scale up to large genomes with moderate memory usage, but the computational resources remained high. The human genome assembly of SGA used more than 1500 CPU hours over 6 days.

"Overlap graph scalability" core messages:

* Overlap graphs and string graphs usually rely on memory heavy data structure

* Compressed data structures may be used but still lead to very long assembly processes

* String graphs grow roughly with the size of the dataset when de Bruijn graph grow with the size of the sequenced genome

## 2.3 The scaling story of the de Bruijn graph representation

The main difference between the overlap graph and the de Bruijn graph is the memory usage. The overlap graph grows linearly with the size of the sequencing data. The de Bruijn graph memory usage can (with kmer filtering) be linear in the size of the genome, since we only have roughly *GenomeSize* different kmers in the genome, estimated through their abundances. A very low abundance threshold would keep too much kmers and
therefore use too much memory. However, even with a very low threshold the number of non genomic kmers is hardly very high (Table 2.1). Thus the memory usage used by a de Bruijn graph is rather function of the size of the genome than of the size of the sequencing dataset.

Largest genomes still lead to graphs with billions of nodes and edges, that are not easy to represent with a decent amount of memory. In this section we concentrate on the problem of storing a huge de Bruijn graph in memory. We present here an overview of the main milestones and a state of the art of the data structures used to represent a de Bruijn graph. We describe the first representations that were based on classical graph implementation and exhibit very high memory footprints. Next we introduce the work of Conway and Bromage that proposed a lower bound on the necessary memory to represent a de Bruijn graph. Then we see how the use of external memory highly impacted the memory usage of most assemblers. And finally we present very efficient techniques based on probabilistic data structures. We also present several works on the scaling of de Bruijn graph assembly on a massive number of cores.

### 2.3.1 Generic graphs

The earliest representations of a de Bruijn graph were based on classical and generic graph representations and were not memory efficient. The first graph representations

<table><caption>Table 2.1: Number of erroneous and missing kmers in a de Bruijn graph after an abundance filtering. On simulated reads with a 100X coverage, 1% error rate and k=51 from the *E. coli*. The *E. coli* genome presents 4,554,207 genomic 51mers. The *C. elegans* genome presents 100,286,051 genomic 51mers.</caption><thead><tr><th>Reference genome</th><th>Abundance Threshold</th><th>Erroneous kmers</th><th>Missing kmer</th></tr></thead><tbody><tr><td><i>E. coli</i></td><td>1</td><td>89,044,977</td><td>5</td></tr><tr><td><i>E. coli</i></td><td>2</td><td>3,889,707</td><td>6</td></tr><tr><td><i>E. coli</i></td><td>3</td><td>200,173</td><td>11</td></tr><tr><td><i>E. coli</i></td><td>4</td><td>15,198</td><td>11</td></tr><tr><td><i>E. coli</i></td><td>5</td><td>1,913</td><td>13</td></tr><tr><td><i>E. coli</i></td><td>6</td><td>312</td><td>16</td></tr><tr><td><i>E. coli</i></td><td>7</td><td>33</td><td>20</td></tr><tr><td><i>E. coli</i></td><td>8</td><td>0</td><td>21</td></tr><tr><td><i>C. elegans</i></td><td>1</td><td>1,909,728,320</td><td>6</td></tr><tr><td><i>C. elegans</i></td><td>2</td><td>89,087,143</td><td>7</td></tr><tr><td><i>C. elegans</i></td><td>3</td><td>7,294,417</td><td>11</td></tr><tr><td><i>C. elegans</i></td><td>4</td><td>2,003,946</td><td>15</td></tr><tr><td><i>C. elegans</i></td><td>5</td><td>1,049,812</td><td>21</td></tr><tr><td><i>C. elegans</i></td><td>6</td><td>643,918</td><td>26</td></tr><tr><td><i>C. elegans</i></td><td>7</td><td>425,435</td><td>26</td></tr><tr><td><i>C. elegans</i></td><td>8</td><td>293,457</td><td>75</td></tr><tr><td><i>C. elegans</i></td><td>9</td><td>208,539</td><td>213</td></tr><tr><td><i>C. elegans</i></td><td>10</td><td>150,939</td><td>589</td></tr></tbody></table>
were based on hash tables indexing kmers associated with informations allowing graph traversal, as in EULER [31] or Velvet [33] assemblers. The use of a dynamic hash table for indexing millions of elements rapidly lead to very high memory usage. Velvet's paper reports a memory usage of 2GB for a Streptococcus genome of 2.2 millions bases. The proposed solution considered for scaling was the use of disk in order to store the graph structure. However this kind of solutions was never really put in practice since disk accesses are very slow compared to RAM accesses. Thus such methods would present orders of magnitude longer running time. Such data structure is adapted to bacterial genomes assembly but impracticable on larger one.

In order to address this scalability issue, Abyss [40] introduced several solutions. First the indexing method changed, using an interesting property of the de Bruijn graph: the numbers of fathers and of sons of a node is bounded by four (the size of the nucleic alphabet). In Abyss, all kmers are indexed and the value associated is just eight bits coding the presence or absence of its eight possible neighbors. Secondly, the main cost of a hash table is usually the pointers used for its internal structure. Instead of this, they used open addressing hash tables that represent a low overhead of a few bits by elements. Thirdly they proposed the distribution of the graph on several machines. The idea of their distributed de Bruijn graph is to use a first hash to decide on which hashtable/machine the kmer will be inserted. On large genome, memory usage was still high and some informations had to be transferred via the network between machines. Still, Abyss was the first assembler scalable enough to assemble a human genome. Today some assemblers still require Terabyte level of memory to assemble large genomes as ABYSS that assembled the white spruce using 4.3TB of memory distributed on large memory servers [61] or DiscovarDenovo that used almost 2 terabytes for a human assembly [62]. The main reason is their too memory expensive indexing methods. In the following we present several works that assess the problem of designing smaller data structures to make large genomes assembly more tractable.

## 2.3.2 Theoretical limits

Conway and Bromage [63] analyzed the theoretical minimal number of bits necessary to encode a de Bruijn graph. They first noted that only the nodes need to be encoded, since edges could simply be deduced from the presence of nodes. Thus a de Bruijn graph can be defined as its set of kmers. Therefore the minimal number of bits needed to encode a de Bruijn graph of $n$ kmers is

$$ bits = \log\left(\binom{4^k}{n}\right) $$

That is $\Omega(n.log(n))$ if $4^k > n$. The paper gives the example of a human genome with 5 billions 25-mer that needs at least 12 GB of memory. A regular bit array that may contain all possible kmers would be impossible to store for $k$ larger than 20. Such bit array would be extremely sparse. They therefore propose the use of succinct data structures to represent such arrays that can be extremely compressed based on Elias-Fano encoding (Sarray) [64]. They provide an implementation of this technique in a tool called Gossamer [43]. It achieved a memory usage close to their proposed bound and a proof of concept assembly of a human genome.
### 2.3.3 Kmer Counting

One may be surprised by the apparent contradiction between the assertion that a de Bruijn graph assembler use memory according to the number of genomic kmer and the fact that previously presented tools showed so high memory usage. The problem is that such tools have to deal with all kmers from the reads, not just the genomic ones. Because of sequencing errors, the number of different kmers before filtering is huge. For example a simulated 100X coverage dataset of *E. coli* with 1% error rate contains 79,917,279 erroneous kmers for only 4,554,207 genomic ones. In order to avoid the indexing of a huge amount of erroneous kmer, a new category of tools called kmer counters has been proposed. The goal of a kmer counter is to associate to each kmer of a dataset its abundance, in order to keep only the solid ones for the assembly, using the lowest amount of resources. Jellyfish [65] proposed the use of an open addressing hash table [66]. The use of this pointer-less table allows the counting operation to use order of magnitude less memory.

DSK [67] proposed a disk based algorithm with an extremely low memory footprint and achieved to count the kmers of a human genome using only 4 GB of memory where JellyFish used 70GB. Although the use of disk operations slowed the process, using 18 hours instead of 3.5 for Jellyfish. However the use of Solid State Drive (SSD) allows DSK to be as fast as jellyfish. KMC2 [68] improved the DSK paradigm by compacting some kmers together reducing the disk footprint. This approach enables a great speedup and is even faster than JellyFish.

These algorithms allowed the filtering of a huge part of the erroneous kmers with a low memory usage. Following assemblers could rely on such techniques in order to work on a set of mostly genomic kmers and therefore present order of magnitude lower memory usage.

### 2.3.4 Probabilistic de Bruijn graphs

Several approaches proposed the use of probabilistic structures in order to leverage the memory footprint. SparseAssembler [69] or LightAssembler [70] propose a sub-sampling of kmers based on the fact that most genomic kmers will appear a high number of times. They index only a subsample, storing 1 out of *g* read kmers. It ensues that most genomic kmer are still contained in the sub-sample, while the number of erroneous kmer is roughly divided by *g*. A more aggressive sub-sampling can be done on high coverage. LightAssembler proposes by default *g* = 3 with 25X and a sub-sampling up to *g* = 25 for 280X. Even if those parameters are conservative, such techniques are still more likely to lose kmers than exact approaches because they highly rely on uniform coverage and error distribution. Such techniques may be criticized for not using the whole information present in their data and possibly missing genomic kmers.

Another proposition was to index the kmers thanks to a Bloom filter [71]. The Bloom filter is a probabilistic set presenting false positives but no false negatives, and having the property of using very low memory. A query to a Bloom filter is of the type "Is this kmer in the set ?". If the answer is "no", then we are sure that the kmer is not in the set, meaning there is no false negative. If the answer is "yes", then the queried kmer is likely to be in the set. A query on a kmer not inserted in the set may return "yes". Such errors
<figure><img src="image_37.png" /><figcaption>Figure 2.4: How Minia represent an exact de Bruijn graph storing a set of critical false positive. Source http://minia.genouest.org/. The Bloom filter creates false positives in black and red. False positives not connected to the graph do not matter since they will not be queried. Only false positives that are connected to the graph are indexed. This way an assembly on the real kmer can be performed.</figcaption></figure>

are false positives. The probability of such an error to appear can be controlled by using
more memory, a dozen of bits per kmer can lead to a false positive rate around 1%. In the
context of assembly, kmers are usually inserted into the Bloom filter during an indexing
phase and the filter is queried during the assembly step. Most of the false positives are
harmless to the assembly process, most of them are kmers disconnected from the graph
and will never be visited or queried. But in some rare cases, some false positive kmers
can be connected to the graph and induce errors or fragmentations.

An improvement was proposed in Minia [72], that stores those "critical false positives" in order to obtain a lossless de Bruijn graph representation (Figure 2.4). This structure was itself improved using several Bloom filters in order to store the critical false positives [73]. Minia successfully assembled a human genome with less than 6 GB and less than one day with one processor. On the same dataset, Gossamer used more than 30 GB and 50 hours and AByss more than 300 GB and 15 hours.

2.3.5 Navigational data structures

One may be surprised to see that Minia achieved the assembly of a human genome with
only 6 GB (less than 16 bits per kmer) where the theoretical limit was supposed to be the
double. The Conway and Bromage theoretical limit was the number of bits necessary to
index a set of kmer. The lower memory usage of Minia is due to the fact that it relies on a
data structure specific to assembly. Minia structure may answer wrongly to membership
query, because the storing of the critical false positives are only designed to block the
assembler to "exit" from the de Bruijn graph. The Minia structure is therefore not an
exact representation of the set of kmers which explains why it is able to use less memory
than the theoretical limit. Minia takes advantage of the fact that during the assembly step, kmer membership queries are not random as we are interested in the existence of a neighbor in order to go through the paths of the graph. The membership operation itself is thus not primordial. Minia is able to answer correctly to neighbor-ship queries (membership of neighbors of a node), used for assembly, with a lower memory usage than needed to index exactly all kmers. We call such representation a "navigational data structure" as introduced by [74] as opposed to a "membership data structure". This paper introduces the notion of Navigational Data Structure (NDS) and shows that it can be used to perform de Bruijn graph assembly. The interest of such structures is the lower bound on their sizes, a NDS needs at least 3.24 bits per kmer to represent a de Bruijn graph. Future assemblers could rely on implementation of such data structure to provide highly reduced memory usage.

### 2.3.6 Massively parallel assembly

In previous sections we talked mostly about the memory bottleneck of assembly. But the time required to perform the assembly step can also be a limitation in many analyses, especially for large genomes or large collections of genomes.

Several contributions were proposed in order to offer fast genome assembly thought the use of massive multi-core servers. Since most assembly algorithms consist mainly in graph traversal operations, it is difficult to really increase the throughput by optimizing these operations since they essentially rely on memory accesses. The main way to improve the wall clock time is to perform the graph traversal operations in parallel. Some assemblers are able to use multiple threads to accelerate their processing but the gain is usually limited. In its benchmark, HipMer [75] performed a human assembly in less than one hour on a thousand cores server while other parallel assemblers like Abyss were not able to make an efficient use of such architectures and were more than ten time slower. The challenge in order to enable massive parallelization is to limit the threads communications. The question of memory access is also critical, when the massive parallelization requires the use of multiple machines communicating via very slow network accesses. Massively parallel assemblers rely on efficient graph partitioning strategies following the divide and conquer paradigm in order to improve data locality and to limit synchronization messages between threads dealing with a local task. If tools like HipMer or SWAP [76, 77] are able to use thousands of cores, those optimizations lead to a very high memory usage and more generally very high resource consumption.

"The scaling story of the de Bruijn graph representation" core messages:

* Regular data structures like hash table do not scale on large genome

* Efficient use of external memory to filter erroneous kmer help assemblers to index mostly genomic kmer

* Memory efficient and specific data structures can lead to moderate memory usage to assemble large genomes
<table><caption>Table 2.2: Statistics on unitigs created from references genomes with k=63. The N50 value is a metric to evaluate an assembly. Given a set of sequences of varying lengths, the N50 length is defined as the length *N* for which 50% of all bases in the sequences are in a sequence of length *L* < *N*. Intuitively, we can cover 50% of the assembly with sequences larger than the N50 value.</caption><thead><tr><th>Reference genome</th><th># Unitig</th><th>N50</th><th># kmer</th></tr></thead><tbody><tr><td><i>E. coli</i></td><td>1,520</td><td>67,344</td><td>4,567,544</td></tr><tr><td><i>C. elegans</i></td><td>127,106</td><td>13,605</td><td>96,501,920</td></tr><tr><td>Human</td><td>2,755,964</td><td>4,967</td><td>2,768,098,045</td></tr></tbody></table>

## 2.4 Efficient de Bruijn graph representation

In this section we describe our contribution in the paper "On the representation of de Bruijn graphs" [74]

### 2.4.1 Compacted de Bruijn graph

In practice, an assembly de Bruijn graph can often be decomposed into a set of long simple paths. A simple path is a path whose nodes exhibit in and out degrees of one, except the first and last nodes. A maximal simple path is a simple path that is not included in a larger simple path. A nodes-disjoint set of such maximal simple paths that cover the de Bruijn graph can represent its set of kmers (Figure 2.5). We call the simple paths from such a set "unitigs". But each unitig composed of *p* kmers can be represented as a string of length *p* + *k* − 1. We can therefore represent a de Bruijn graph of *n* kmers with less than

$$2(n + (k - 1)\#\text{unitigs}) \text{ bits}$$

$2n$ bits for the genome sequences and an overhead of $2(k - 1)$ bits for each unitig. This graph of unitigs is called the compacted de Bruijn graph. Unitigs construction is often the first step of the assembly process. The unitig set can be considered as a safe assembly, since graph modification heuristics are applied on the unitigs in order to get longer contigs sequences. Two points can be observed in Table 2.2:

* A genome presents orders of magnitude less unitigs than kmers

* Unitigs can represent quite large sequences

Those points show the interest to consider unitigs instead of kmers as nodes of the graph.

### 2.4.2 De Bruijn graph construction

The question asked here is "Can we take advantage of the unitig representation in order to propose a data structure that achieves a memory consumption close to 2 bits per kmers?" The paper [74] propose a data structure called "DBGFM" to represent a De Bruin graph in low memory. The idea is to index the de Bruijn graph unitigs in a FM-index [58]. The
<figure><img src="image_38.png" /><figcaption>Figure 2.5: Toy example of a de Bruijn graph (left) and the corresponding compacted de Bruijn graph (right). The nodes of the compacted de Bruijn graph are no longer kmers but sequences of length ≥ *k* and called unitigs.</figcaption></figure>
FM-index is a lightweight full-text index based on the Burrows-Wheeler transform [78]. Such a structure can answer membership and neighborhood queries. It allows to count the number of occurrences of a pattern *q* in *O(|q|)*. We can therefore perform membership queries of a kmer in *O(|k|)*. But the main interest is that we can access the in-neighbors of a node in constant time by querying the symbols preceding a kmer. Hence, **DBGFM** is a membership data structure that also provides fast in-neighbors query.

The **DBGFM** structure was integrated into the ABySS assembler as a proof of concept of low memory assembly. But straightforward compacted de Bruijn graph construction requires to store the de Bruijn graph itself in memory. This de Bruijn graph compaction step would therefore become the bottleneck of the assembly workflow. To fulfill this need, we proposed a new external memory based algorithm in order to compute the compacted de Bruijn graph with low memory called **BCALM**.

BCALM algorithm is based on the idea of minimizer [79]. The *l* minimizer of a string *u* is the smallest *l*mer that is a substring of *u* according to a total ordering of the strings e.g. lexicographical. We define the left minimizer Lmin(*u*) and the right minimizer Rmin(*u*) as the minimizer of the *k* − 1 prefix of *u* and the minimizer of the *k* − 1 suffix of *u* respectively. We use the binary relation *u* → *v* between two strings that denotes an exact suffix-prefix overlap of length *k* − 1 between *u* and *v*. The use of minimizers is motivated by the following property: For two strings *u* and *v*, if *u* → *v* then Rmin(*u*)=Lmin(*v*)

For describing the BCALM algorithm we will rely on the notions of compatibility introduced in the next paragraph.

Given a set of strings $S$, we say that $(u, v) \in S^2$ are *compactable* in a set $V \subseteq S$ if $u \to v$ and, $\forall w \in V$, if $w \to v$ then $w = u$ and if $u \to w$ then $w = v$. The compaction operation is defined on a pair of compactable strings $u, v$ in $S$. It replaces $u$ and $v$ by a single string $w = u \cdot v[k+1...|v|]$ where '.' is the string concatenation operator. We say that two strings $(u, v)$ are $m$-compactable in $V$ if they are compactable in $V$ and if $m = Rmin(u) = Lmin(v)$. The $m$-compaction of a set $V$ is obtained by applying the compaction operation as much as possible in any order to all pairs of strings that are $m$-compactable in $V$. It is easy to show that the order in which strings are compacted does not lead to different $m$-compactions. Compaction is a useful notion because a simple way to obtain the unitigs is to greedily perform compaction as long as possible.

2.4.3 Unitig enumeration in low memory: BCALM

The BCALM algorithm is described in Figure 2.6, and a step by step execution is pre-
sented in Figure 2.7. The first step consists into placing the input kmers into different
files $F_m$ according to their minimizers. Then, each file is processed starting from the file
with the smaller minimizer in increasing order. Each file $F_m$ is loaded in memory and
all possible $m$-compaction among the sequences of the file are applied. The idea is that
each file will be small enough to be compacted efficiently using low amounts of memory.
Each sequence of the file are thereafter placed in the output file or in another file to be
eventually further compacted. The rules of choosing which file to write a sequence is
based on its minimizers. If both the left and right minimizers are below $m$, the sequence
is output as a unitig in the output file. Otherwise we write the sequence in the file $F_{m'}$
where $m'$ is the smallest minimizer bigger than $m$ of the sequence. Finally the file $F_m$ is
discarded and the next file is processed.
1: **Input:** Set of kmers *S*, minimizer size *l* < *k*

2: **Output:** Sequences of all simple paths in the de Bruijn graph of *S*

3: Perform a linear scan of *S* to get the frequency of all *l*-mers (in memory)

4: Define the ordering of the minimizers, given by their frequency in *S*

5: Partition *S* into files *F*<sub>*m*</sub> based on the minimizer *m* of each *k*-mer

6: **for** each file *F*<sub>*m*</sub> in increasing order of *m* **do**

7: &nbsp;&nbsp;&nbsp;&nbsp;*C*<sub>*m*</sub> ← *m*-compaction of *F*<sub>*m*</sub> (performed in memory)

8: **for** each string *u* of *C*<sub>*m*</sub> **do**

9: &nbsp;&nbsp;&nbsp;&nbsp;*B*<sub>min</sub> ← min(Lmin(*u*), Rmin(*u*))

10: &nbsp;&nbsp;&nbsp;&nbsp;*B*<sub>max</sub> ← max(Lmin(*u*), Rmin(*u*))

11: **if** *B*<sub>min</sub> ≤ *m* and *B*<sub>max</sub> ≤ *m* **then**

12: &nbsp;&nbsp;&nbsp;&nbsp;Output *u*

13: **else if** *B*<sub>min</sub> ≤ *m* and *B*<sub>max</sub> > *m* **then**

14: &nbsp;&nbsp;&nbsp;&nbsp;Write *u* to *F*<sub>*B*<sub>*max*</sub></sub>

15: **else if** *B*<sub>min</sub> > *m* and *B*<sub>max</sub> > *m* **then**

16: &nbsp;&nbsp;&nbsp;&nbsp;Write *u* to *F*<sub>*B*<sub>min</sub></sub>

17: **end if**

18: **end for**

19: Delete *F*<sub>*m*</sub>

20: **end for**

<caption>Figure 2.6: BCALM: Enumeration of all maximal simple paths in the De Bruijn graph</caption>

The idea of the algorithm is to perform the *m*-compactions for each *m* in order to perform all compactions. The intuition of the rule to place sequences in the correct file is that, if the suffix (or the prefix) has a overlap with a minimizer *n* superior to *m*, then it will eventually be compacted in the file *F*<sub>*n*</sub>. If both minimizers are below *m*, no further compactions can be made on the sequence and it is therefore an unitig. If both minimizers are above *m* then the smallest is chosen, as the algorithm treats minimizers in increasing order.

Several implementation details allow this algorithm to be practical. Reverse comple-
ments are handled as presented before by identifying each kmer with its reverse com-
plement and letting the minimizer be the smallest lmer in both of them. To avoid the
creation of too many files, we encode several virtual files in one physical one, this allowed
to use l = 10 in our experiments. We also avoid to load in memory the whole sequences of
a file as only the k - 1 prefix and suffix are used by the compaction detection algorithm.
For a fixed input S, the number of strings in a file F m depends on the minimizer length
l and the ordering of minimizers. When l increases, the number of k − 1mers in S that
share a given minimizer decreases. Thus, increasing l yields less strings per file, which
decreases the memory usage. We realized that, when highly-repeated lmers are less likely
to be chosen as minimizers, the sequences are more evenly distributed among files. We
therefore perform in-memory lmer counting (line 3) to obtain a sorted frequency table
of all lmers. Each lmer is then mapped to its rank in the frequency array, to create a
total ordering of minimizers (line 4). Our experiments showed a drastic improvement
over lexicographic ordering (Table 2.3).

A proof of correctness of this algorithm is showed in [74].
<figure><img src="image_39.png" /><figcaption>Figure 2.7: An execution of the BCALM algorithm on a toy example. (Panel A) The set S = ACA, AAC, CAA, CAG of k-mers is represented as a de Bruijn graph (k = 3). (Panel B) The frequencies of each l-mer present in the kmers is given in increasing order (l = 2). (Panel C) Steps of the algorithm; initially, each set Fm contains all kmers having minimizer m. For each minimizer m (in order of frequency), perform the m-compaction of Fm and store the result in Cm; the grey arrows indicate where each element of Cm is written to, either the file of a higher minimizer or the output. The row G(Fm) shows a graph where solid arrows indicate m-compactable pairs of strings in Fm. The dash-dot arrow in Step 3 indicates that the two strings are compactable in Fm but not m-compactable; in fact, they are not compactable in S.</figcaption></figure>
<table><caption>Table 2.3: Memory usage of BCALM with three different minimizer orderings: lexicographical, uniformly random, and according to *l*mer frequencies. The dataset used is the human chromosome 14 with *k* = 55 and *l* = 8</caption><thead><tr><td>Ordering used</td><td>Lexicographical</td><td>Uniformly random</td><td>lmer frequency</td></tr></thead><tbody><tr><td>Memory usage</td><td>840 MB</td><td>222 MB</td><td>19 MB</td></tr></tbody></table>

### 2.4.4 Assembly in low memory using BCALM

We tested the efficiency of our algorithm in two datasets. The first dataset is 36 million 155bp Illumina human chromosome 14 reads (2.9 GB compressed fastq) from the GAGE benchmark [80]. The second dataset is 1.4 billion Illumina 100bp reads (54 GB compressed) from the NA18507 human genome (SRX016231). We first processed the reads with kmer counting software. We used DSK for this operation for its low memory usage. We set *k* to 55 and applied a solidity threshold of 5 for chr14 and 3 for the whole genome.

First we evaluate (Table 2.3) the memory usage of BCALM on the chr14 dataset, according to the minimizer ordering used. We observe that the lexicographical order do not perform well. The lexicographical smallest *l*mer is *A*<sup>*l*</sup>, which is abundant in the human chromosomes, resulting in a large *F*<sub>*m0*</sub> and a high memory usage. A second attempt has been made with a uniformly random ordering. If it performs better than the lexicographical one, it is still sensible to repeated *l*mer in the first files resulting in high memory usage. Finally, we ordered the *l*mers according to their frequency in the dataset. This resulted in a memory usage of 19 MB, a 40-fold improvement over the initial lexicographical ordering. The running times of all three orderings were comparable.

We also evaluated the effect of the size of the minimizer on the BCALM performances in Table 2.4. Large *l* generally lead to small memory usage, however we did not see much improvement past *l* = 8 on this dataset.

Using BCALM, we evaluated the performances of the low memory assembly workflow with DSK, BCALM and DKGFM in Table 2.5. For the whole genome dataset, BCALM used only 43 MB of memory to compact a set of 2.5 * 10<sup>9</sup> 55mers and output 40 million sequences of total length 4.6 billions base pairs. DKGFM represented these paths in an index of size 1.5 GB. The overall construction time, including DSK, was roughly 24 hours. In comparison, a subset of this dataset was used to construct the data structure of Salikhov et al [73] in 30.7 hours.

We also compared the memory usage of the proposed workflow with the other state of the art approaches in Table 2.6, showing the interest the proposed low memory assembly workflow.

Those results show that indexing unitigs instead of kmers can lead to a memory efficient assembly steps. The BCALM algorithm, allows an memory efficient construction of the unitig set by relying on external memory.
<table><caption>Table 2.4: Memory usage and wall clock time of BCALM with increasing values of minimizer sizes *l* on the chr14 data. With the system of virtual files, these values of *l* require respectively 4,16,64,256 and 1024 physical files. The used ordering is the one based on *l*-mer counts.</caption><tr><td>Minimizer size <i>l</i></td><td>2</td><td>4</td><td>6</td><td>8</td><td>10</td></tr><tr><td>Memory usage</td><td>9,879 MB</td><td>3,413 MB</td><td>248 MB</td><td>19 MB</td><td>19 MB</td></tr><tr><td>Running time</td><td>27m19s</td><td>22m2s</td><td>20m5s</td><td>18m39s</td><td>21m4s</td></tr></table>

<table><caption>Table 2.5: Running times (wall-clock) and memory usage of DSK, BCALM andDBGFM construction on the human chromosome 14 and whole human genome datasets (*k* = 55 and *l* = 10 for both).</caption><thead><tr><th>Dataset</th><th>DSK</th><th>BCALM</th><th>DBGFM</th></tr></thead><tbody><tr><td rowspan="2">Chromosome 14</td><td>43 MB</td><td>19 MB</td><td>38 MB</td></tr><tr><td>25 mins</td><td>15 mins</td><td>7 mins</td></tr><tr><td rowspan="2">Whole human genome</td><td>1.1 GB</td><td>43 MB</td><td>1.5 GB</td></tr><tr><td>5 h</td><td>12 h</td><td>7 h</td></tr></tbody></table>

<table><caption>Table 2.6: Memory usage of de Bruijn graph data structures, on the human chromosome 14 and whole human genome datasets (*k* = 55 for both). We did not run the algorithm of Conway and Bromage because our machine does not have sufficient memory for the whole genome. Instead, we report the theoretical size of their data structure, assuming that it would be constructed from the output of DSK. As described in [81], this gives a lower bound on how well their implementation could perform.</caption><thead><tr><th></th><th>DBGFM</th><th>Salikhov et al.</th><th>Conway & Bromage</th></tr></thead><tbody><tr><td>chr14</td><td>38.0 MB</td><td>94.9 MB</td><td>&gt; 875 MB</td></tr><tr><td>Full human dataset</td><td>1,462 MB</td><td>2,702 MB</td><td>&gt; 22,951 MB</td></tr></tbody></table>
"Efficient de Bruijn graph representations" core messages:

* Unitigs (maximal simple paths) are great abstractions to represent the de Bruijn graph

* Unitigs can be constructed with a low memory usage

* Indexing unitigs with a compressed index achieves very low memory assembly (less than 1.5GB for a human assembly)

## 2.5 Efficient de Bruijn graph construction

> In this section we describe our contribution in the paper "Compacting de Bruijn graphs from sequencing data quickly and in low memory" [82]

We introduced an algorithm to construct the de Bruijn graph with very low memory usage relying on external memory. But the approach does not exhibit straightforward parallel patterns to exploit and is quite slow. In the proof of concept human assembly, BCALM took 12 hours to complete for a total running time of 24 hours for the whole workflow. As we said before, the kmer counting step tool reached a very high throughput and an extremely reduced memory footprint. The graph compaction step has therefore become the bottleneck for assembly. We showed before that unitigs could be used as atoms for assembly instead of kmers and it appears that most assemblers rely on them or even, in some cases, explicitly compute such sets of sequences [83, 40]. Some papers proposed solutions for fast graph compaction but are mostly based on parallel graph traversal and no memory efficient method have been proposed. This is why we propose a fast and parallel algorithm for graph compaction presenting a low memory footprint called BCALM2.

### 2.5.1 BCALM2 algorithm

We give here a high level overview of the BCALM2 algorithm (Figure 2.8). As in BCALM the first step of the algorithm is to distribute the input set of kmer into files. In BCALM2 a kmer can be inserted in two different files, if it has different left and right minimizers. The second step will consist into performing compactions on all the files in parallel according to the compaction procedure described in Figure 2.9. The last stage will need to glue back together the kmers that were duplicated in two files. The Figure 2.12 shows an execution of BCALM2 on a toy example.

Since the size of the files are small, the compactions can be made with a in-memory algorithm in a parallel way. The resulting sequences are written on disk and are processed in the third stage when all compactions are finished. At the end of the compaction step, we notice that some kmers exist in two copies in the sequences produced. Such kmers are always the prefix or the suffix of the compacted sequences. We record the sequences ends that have a doubled kmers has "lonely" as they should be reunited in the third stage. The other sequences are output as unitigs.
In the last stage, we process the sequences that need to be reunited with the Reunite procedure (Figure 2.10). The goal of Reunite is to associate the sequences sharing a suffix or prefix kmer and to glue them together (Figure 2.11) by reuniting the two occurrences of a doubled kmer. The lonely mark are transfered to the produced sequences and this process is repeated until a sequence has a lonely end. After the reunite operation, there is no duplicated kmer left and the output correspond to the unitig set. In order to be efficient, the reunite operation partitions the sequences in order to process the partitions in a parallel and memory efficient way.

A proof of correctness of this algorithm is available in [82]

## 2.5.2 Implementation details

In this section we describe optimizations and important implementation details. As in BCALM, we do not use the lexical ordering for minimizers but the frequency based minimizer order described previously. Files used in the algorithm are virtual files organized into groups, in order to introduce natural checkpoints in BCALM 2 in between parallel sections. BCALM2 iterates sequentially through the groups, but parallelizes the processing within a group. The For loop at line 1 of Figure 2.8 is executed in parallel within a

**Input:** the set of k-mers K.

1. **for all parallel** $x \in K$ **do**

2. Write $x$ to $F(\text{Lmin}(x))$.

3. **if** $\text{Lmin}(x) \neq \text{Rmin}(x)$ **then**

4. Write $x$ to $F(\text{Rmin}(x))$.

5. **end if**

6. **end for**

7. **for all parallel** $i \in \{1, \dots, 4^l\}$ **do**

8. Run *CompactBucket*($i$)

9. **end for**

10. *Reunite()*

<caption>Figure 2.8: BCALM2($K$)</caption>

1. Load $F(i)$ into memory.

2. $U \leftarrow i$-compaction of $F(i)$.

3. **for all** strings $u \in U$ **do**

4. Mark $u$'s prefix as “lonely” if $i \neq \text{Lmin}(u)$.

5. Mark $u$'s suffix as “lonely” if $i \neq \text{Rmin}(u)$.

6. **if** $u$'s prefix and suffix are not lonely **then**

7. Output $u$.

8. **else**

9. Place $u$ in the Reunite file

10. **end if**

11. **end for**

<caption>Figure 2.9: *CompactBucket*($i$)</caption>
**Input:** the set of strings $R$ from the Reunite file.

1. UF $\leftarrow$ Union find data structure whose elements are the distinct $k$-mer extremities in $R$.

2. **for all parallel** $u \in R$ **do**

3. **if** both ends of $u$ are lonely **then**

4. $UF.union(suf_k(u), pre_k(u))$

5. **end if**

6. **end for**

7. **for all parallel** classes $C$ of UF **do**

8. $P \leftarrow$ all $u \in R$ that have a lonely extremity in $C$

9. **while** $\exists u \in P$ that does not have a lonely prefix **do**

10. Remove $u$ from $P$

11. Let $s = u$

12. **while** $\exists v \in P$ such that $suf_k(s) = pre_k(v)$ **do**

13. $s \leftarrow Glue(s, v)$

14. Remove $v$ from $P$

15. **end while**

16. Output $s$

17. **end while**

18. **end for**

<caption>Figure 2.10: Reunite()</caption>

**Input:** strings $u$ and $v$, such that $suffix_k(u) = prefix_k(v)$.

1. Let $w = u \odot^k v$.

2. Set lonely prefix bit of $w$ to be the lonely prefix bit of $u$.

3. Set lonely suffix bit of $w$ to be the lonely suffix bit of $v$.

4. **return** $w$

<caption>Figure 2.11: Glue($u, v$)</caption>
<figure><img src="image_40.png" /><figcaption>Figure 2.12: Execution of BCALM2 on a small example, with $k = 4$ and $l = 2$. On the top left, we show the input de Bruijn graph. The maximal unitig corresponds to the path from CCCT to TCTA (spelling CCCTCTA), and to the kmers CCCC, CCCA, CTAC, CTAA. In this example, minimizers are defined using a lexicographic ordering of lmers. In the top right, we show the contents of the files. Only five of the files are non-empty, corresponding to the minimizers CC, CT, AA, AC and CA. The doubled kmers are italicized. Below that, we show the set of strings that each i-compaction generates. For example in the file CC, the kmers CCCT and CCTC are compacted into CCCTC, however CCCC and CCCT are not compactable because CCCA is another out-neighbor of CCCC. The lonely ends are denoted by $\bullet$. In the bottom half we show the execution steps of the Reunite algorithm. Nodes in bold are output</figcaption></figure>

group, with each thread given a subset of $K$. Kmers are distributed only to those files that are in the group, with other buckets being ignored. After the kmers are distributed, files from a group are compacted in parallel. The CompactBucket routines are independent of each other, and hence we run CompactBucket(i) in parallel using all available processors. After BCALM2 finishes processing a group, it moves on to the next group. For the Reunite operation, we used a minimal perfect hash function of all distinct kmer extremities in order to perform the necessary operations. The possibly high amount of key to index lead us to the conception of a new method to construct such functions that will be described in the next section. The presented algorithm takes as input a set of kmers, but in our implementation, BCALM2 is based on the GATB [84] library and make use of its kmer counting operation to handle directly read files. In this kmer counter, kmers are divided into partitions according to their minimizer, then each partition is counted independently. We modified the GATB kmer counting algorithm so that partition files correspond exactly to file groups.

### 2.5.3 Large genome de Bruijn graphs

We evaluated BCALM2 performances and its comparison with other tool for compacting a De Bruijn graph. We used two human datasets from GAGE [80] and two larger datasets from the spruce and pine sequencing projects [61, 85].

In Figure 2.13 we show how BCALM2 is affected by the changes in the parameters $k$ and $l$. It shows that BCALM2 has almost identical running times for $l$ from 6 to 10.
<figure><img src="image_41.png" /><figcaption>Figure 2.13: BCALM2 wall-clock running times with respect to parameters *l* (left) and *k* (using 4 cores) and number of cores (using *k* = 55 and *l* = 8), on the chromosome 14 dataset (right).</figcaption></figure>

Short minimizer sizes create fewer files and therefore limit parallel operations. We also show how BCALM2 scales with multiples of cores. We observe that the compactions scale in a nearly linear manner with the number of thread used. There remains overhead related to disk operations.

We also compared the BCALM2 performances with BCALM, ABYSS and Merac-
ulous2 graph compaction steps in Table 2.7. We observe that BCALM2 is the fastest
available tool while using order of magnitude less memory but compared to BCALM.
We further evaluated BCALM 2 on two very large sequencing datasets: Illumina reads
from the 20 Gbp Picea glauca genome (8.5 billion reads, 152–300 bp each, 1.1 TB com-
pressed FASTQ, SRA056234), and Illumina paired-end reads from the 22 Gbp Pinus taeda
genome (9.4 billion reads, 128–154 bp each, 1.2 TB compressed FASTQ, SRX016231).
The k-mer counting step took around a day and less than 40 GB of memory for each
dataset.

Table 2.8 shows the performance of BCALM 2 on these two datasets, as well as
unitigs statistics. Graph construction of the spruce dataset previously required 4.3 TB
of memory and 2 days on a 1380-core cluster [61], while the assembly of the pine dataset
previously required 800 GB of memory and 3 months on a single machine [85]. Although
we used the same sequencing datasets, several parameters differ between these previous
reports and our results (e.g. k value, abundance cutoff, and whether reads were error-
corrected). Hence, run time, memory usage, and unitigs statistics cannot be directly
compared. However, it seems reasonable to infer that BCALM 2 would remains 1–2
orders of magnitude more efficient in time and memory.
<table><caption>Table 2.7: Running times (wall-clock) and memory usage of compaction algorithms for the human datasets. For BCALM2 and BCALM we used $k = 55$, and $l = 8$ and $l = 10$, respectively; abundance cutoffs were set to 5 for Chr 14 and 3 for whole human. We used 16 cores for the parallel algorithms ABySS, Meraculous 2, and BCALM2. Meraculous 2 aborted with a validation failure due to insufficient peak k-mer depth when we ran it with abundance cutoffs of 5. We were able to execute it on chromosome 14 with a cutoff of 8, but not for the whole genome. The exact memory usage was unreported there but is less than < 1 TB. Meraculous 2 was executed with 32 prefix blocks.</caption><thead><tr><th>Dataset</th><th>BCALM2</th><th>BCALM</th><th>ABySS-P</th><th>Meraculous 2</th></tr></thead><tbody><tr><td>Chr 14</td><td>5 mins<br>400 MB</td><td>15 mins<br>19 MB</td><td>11 mins<br>11 GB</td><td>62 mins<br>2.35 GB</td></tr><tr><td>Whole<br>human</td><td>1.2 h<br>2.8 GB</td><td>12 h<br>43 MB</td><td>6.5 h<br>89 GB</td><td>16 h *<br>unreported *</td></tr></tbody></table>

<table><caption>Table 2.8: Performance of BCALM2 on the loblolly pine and white spruce datasets. The kmer size was 31 and the abundance cutoff for kmer counting was 7.</caption><thead><tr><th>Dataset</th><th></th><th>Loblolly<br>pine</th><th>White<br>spruce</th></tr></thead><tbody><tr><td>Distinct<br>(×10<sup>9</sup>)</td><td>k-mers</td><td>10.7</td><td>13.0</td></tr><tr><td>Num threads</td><td></td><td>8</td><td>16</td></tr><tr><td>CompactBucket()</td><td>time</td><td>4 h 40 m</td><td>3 h 47 m</td></tr><tr><td>CompactBucket()</td><td>mem</td><td>6.5 GB</td><td>6 GB</td></tr><tr><td>Reunite file size</td><td></td><td>85 GB</td><td>140 GB</td></tr><tr><td>Reunite()</td><td>time</td><td>4 h 32 m</td><td>3 h 08 m</td></tr><tr><td>Reunite()</td><td>memory</td><td>31 GB</td><td>39 GB</td></tr><tr><td>Total time</td><td></td><td>9 h 12 m</td><td>6 h 55 m</td></tr><tr><td>Total max memory</td><td></td><td>31 GB</td><td>39 GB</td></tr><tr><td>Unitigs (×10<sup>6</sup>)</td><td></td><td>721</td><td>1200</td></tr><tr><td>Total length</td><td></td><td>32.3 Gbp</td><td>49.0 Gbp</td></tr><tr><td>Longest unitig</td><td></td><td>11.2 Kbp</td><td>9.0 Kbp</td></tr></tbody></table>
"Efficient de Bruijn graph construction" core messages:

* Graph construction is the bottleneck of most assemblers

* Unitigs can be constructed in parallel and in low memory

* Very large genomes can be addressed without huge need of RAM

## 2.6 Indexing large sets

We presented the interest of working on unitigs for assembly and how to efficiently compute them. We proposed DBGFM, a data structure able to index the set of unitigs using less memory than previously best known data structures and fast enough to allow efficient assembly. Since DBGFM is more powerful than a NDS because it is a membership data structure, can we propose a better structure that is just a NDS?

We are interested in a structure using less memory and providing faster access to neighbors. We remark that the unitig sequence representation is close to the 2 bits per kmer limit. Each unitig represents an overhead of $2 * (k - 1)$ nucleotides to add to the 2 bits per kmer. Raw unitig sequences is not an interesting structure in practice because a query would be linear. Indexing the last and first $k-1$mers of each unitig and associating to each $k-1$mer the unitigs that start and end with it could lead to an NDS. The cost of such a structure would be at least $8 * \log(\#\text{unitigs})$ bits per unitig, as a $k-1$mer can be associated to up to 8 unitigs ($\log(\#\text{unitigs})$ bits are needed to encode the index of an unitig), plus the indexing cost.

A better scheme would be to index the first and last kmer of each unitig and to associate to each indexed kmer its unique unitig. Indeed a kmer appears in at most one unitig due to the definition of the compacted de Bruijn graph. This way, an unitig can know its sons by querying the index with the four possible kmers. In order to avoid performing useless queries, we can add 8 bits to each unitig to precise which neighbors exist. This leads to a structure of at least $2 * \log(\#\text{unitigs})$ per unitig (to encode the unitig index twice) bits plus 8 bits from the proposed optimization plus the cost of indexing. This kind of structure would allow constant and fast queries, with a very low memory overhead. A regular hashing would be costly in memory and in query time. The best solution to keep a low memory usage would be an open addressing hash table that would offer a $2 * \log(\#\text{unitigs}) * (1/\text{loadFactor})$ bits. In practice the load factor is set below 0.8. Practically we can consider that the unitig identifier would be represented by a 32 bits integer, thus the open addressing would cost roughly one byte by unitig and would produce long queries when collisions occur. Can we do better ?

One may argue that we do not use the fact that our set of keys is static and that we could use a Minimal Perfect Hash Function (MPHF). Such structures use very few memory (as low as 3 bits per key) and guarantee constant query time and in practice very fast access (hundred of nanoseconds). Beside, many applications in bioinformatics could benefit from an efficient and lightweight structure to associate information to large sets of keys as kmers or unitigs: coverage, origin of the sequence or any properties. But MPHF are costly to construct for very large sets of keys, both in running time and in memory usage. We therefore present a technique able to efficiently construct efficient
MPHF from huge sets of keys with low memory and time.

> This method is presented in the paper "Fast and scalable minimal perfect hashing construction for massive key sets" [86] included at the end of the chapter.

This paper observes that most methods to construct MPHF are designed to provide the smallest structures. However, most methods use way more memory during construction than the size of the MPHF itself. Generally, available implementations showed very high running time and memory requirement on large sets of keys, making them impracticable to index billions of objects. The paper therefore proposes a new implementation of a simple algorithm providing efficient MPHF that scales on very large sets.

The algorithm used is pretty straightforward. First keys are hashed into a bit array. Second the positions in the bit array that received an unique key are marked with a one and the corresponding keys are removed from the key set. The remaining keys are hashed in a second bit array and the process is repeated until there is not key left. At the end of the algorithm each key is associated to a unique position in one of the bit array. To find this position the key is hashed to the different bit array until a one is found.

The main interests of the method are the amount of memory needed for construction that is similar to the size of the MPHF and the parallel construction allowing a very short running time. BBHASH constructed a MPHF from a billion key dataset in a minute with 2GB of RAM were the best know MPHF library constructs it in more than 20 minutes using more than 18GB. BBHASH also achieved to construct a trillion key MPHF on a 750GB machine. This experiment is, to our knowledge, the largest MPHF constructed.

Applications BBhash MPHF has been used in several bioinformatics tools and was integrated in the GATB [84] library. It is used to associate information to kmers such as coverage or origin. The unification operation of BCALM2 also relies on BBhash MPHF structure. Another tool, called SRC [87, 88], integrated the BBhash MPHF in order to propose an efficient data structure to index reads for read set comparisons and read coverage estimation.

"Indexing large sets" core messages:

* We can build MPHF from huge sets of keys with low resources

* Such MPHF are used for assembly and also for various applications in bioinformatics
# Fast and scalable minimal perfect hashing for massive key sets

Antoine Limasset¹, Guillaume Rizk¹, Rayan Chikhi², and Pierre Peterlongo¹

1 IRISA Inria Rennes Bretagne Atlantique, GenScale team, Campus de Beaulieu 35042 Rennes, France

2 CNRS, CRISTAL, Université de Lille, Inria Lille - Nord Europe, France

## Abstract

Minimal perfect hash functions provide space-efficient and collision-free hashing on static sets. Existing algorithms and implementations that build such functions have practical limitations on the number of input elements they can process, due to high construction time, RAM or external memory usage. We revisit a simple algorithm and show that it is highly competitive with the state of the art, especially in terms of construction time and memory usage. We provide a parallel C++ implementation called *BBhash*. It is capable of creating a minimal perfect hash function of $10^{10}$ elements in less than 7 minutes using 8 threads and 5 GB of memory, and the resulting function uses 3.7 bits/element. To the best of our knowledge, this is also the first implementation that has been successfully tested on an input of cardinality $10^{12}$. Source code: https://github.com/rizkg/BBHash

1998 ACM Subject Classification H.3.1 E.2

Keywords and phrases Minimal Perfect Hash Functions, Algorithms, Data Structures, Big Data

Digital Object Identifier 10.4230/LIPIcs...

## 1 Introduction

Given a set $S$ of $N$ elements (*keys*), a minimal perfect hash function (MPHF) is an injective function that maps each key of $S$ to an integer in the interval $[1, N]$. In other words, an MPHF labels each key of $S$ with integers in a collision-free manner, using the smallest possible integer range. A remarkable property is the small space in which these functions can be stored: only a couple of bits per key, independently of the size of the keys. Furthermore, an MPHF query is done in constant time. While an MPHF could be easily obtained using a key-value store (e.g. a hash table), such a representation would occupy an unreasonable amount of space, with both the keys and the integer labels stored explicitly.

The theoretical minimum amount of space needed to represent an MPHF is known to be $\log_2(e)N \approx 1.44N$ bits [10, 14]. In practice, for large key sets (billions of keys), many implementations achieve less than $3N$ bits per key, regardless of the number of keys [2, 9]. However no implementation comes asymptotically close to the lower bound for large key sets. Given that MPHFs are typically used to index huge sets of strings, e.g. in bioinformatics [6, 7, 8], in network applications [12], or in databases [5], lowering the representation space is of interest. We observe that in many of these applications, MPHFs are actually used to construct static dictionaries, i.e. key-value stores where the set of keys is fixed and never updated [6, 8]. Assuming that the user only queries the MPHF to get values corresponding to keys that are guaranteed to be in the static set, the keys themselves do not necessarily need to be stored in memory. However the associated values in the dictionary typically do need to be stored, and

<figure><img src="image_42.png" /></figure>
XX:2 Fast and scalable minimal perfect hashing for massive key sets

they often dwarf the size of the MPHF. The representation of such dictionaries then consists of
two components: a space-efficient MPHF, and a relatively more space-expensive set of values.
In such applications, whether the MPHF occupies 1.44 bits or 3 bits per key is thus arguably
not a critical aspect.

In practice, a significant bottleneck for large-scale applications is the construction step of
MPHFs, both in terms of memory usage and computation time. Constructing MPHFs efficiently
is an active area of research. Many recent MPHF construction algorithms are based on efficient
peeling of hypergraphs [1, 3, 4, 11]. However, they require an order of magnitude more
space during construction than for the resulting data structure. For billions of keys, while
the MPHF itself can easily fit in main memory of a commodity computer, its construction
algorithm requires large-memory servers. To address this, Botelho and colleagues [4] propose
to divide the problem by building many smaller MPHFs, while Belazzougui et al. [1] propose
an external-memory algorithm for hypergraph peeling. Very recently, Genuzio et al. [11]
demonstrated practical improvements to the Gaussian elimination technique, that make
it competitive with [1] in terms of construction time, lookup time and space of the final
structure. These techniques are, to the best of our knowledge, the most scalable solutions
available. However, when evaluating existing implementations, the construction of MPHFs for
sets that significantly exceed a billion keys remains prohibitive in terms of time and space
usage.

A simple idea has been explored by previous works [6, 12, 16] for constructing PHFs
(Perfect Hash Functions, non minimal) or MPHFs using arrays of bits, or fingerprints. However,
it has received relatively less attention compared to other hypergraph-based methods, and
no implementation is publicly available in a stand-alone MPHF library. In this article we
revisit this idea, and introduce novel contributions: a careful analysis of space usage during
construction, and an efficient, parallel implementation along with an extensive evaluation
with respect to the state of the art. We show that it is possible to construct an MPHF using
almost as little memory as the space required by the final structure, without partitioning the
input. We propose a novel implementation called *BBhash* ("Basic Binary representAtion of
Successive Hashing") with the following features:

- construction space overhead is small compared to the space occupied by the MPHF,

- multi-threaded,

- scales up to to very large key sets (tested with up to 1 trillion keys).

To the best of our knowledge, there does not exist another usable implementation that
satisfies any two of the features above. Furthermore, the algorithm enables a time/memory
trade-off: faster construction and faster query times can be obtained at the expense of a
few more bits per element in the final structure and during construction. We created an
MPHF for ten billion keys in 6 minutes 47 seconds and less than 5 GB of working memory,
and an MPHF for a trillion keys in less than 36 hours and 637 GB memory. Overall, with
respect to other available MPHF construction approaches, our implementation is at least two
orders of magnitudes more space-efficient when considering internal and external memory
usage during construction, and at least one order of magnitude faster. The resulting MPHF
has slightly higher space usage and faster or comparable query times than other methods.

2 Efficient construction of minimal perfect hash function

2.1 Method overview

Our MPHF construction procedure revisits previously published techniques [6, 12]. Given a set
$F_0$ of keys, a classical hash function $h_0$ maps keys to an integer in $[1, |F_0|]$. A bit array $A_0$ of
<figure><img src="image_43.png" /><figcaption>■ **Figure 1** MPHF construction and query example. The input is a set $F_0$ composed of $N = 6$ keys ($k_1$ to $k_6$). All keys are hashed using a hash function $h_0$ and are attempted to be placed in an array $A_0$ at positions given by the hash function. The keys $k_3$ and $k_6$ do not have collisions in the array, thus the corresponding bits in $A_0$ are set to '1'. The other keys from $F_0$ that are involved in collisions are placed in a new set $F_1$. In the second level, keys from $F_1$ are hashed using a hash function $h_1$. Keys $k_1$ and $k_5$ are uniquely placed while $k_2$ and $k_4$ collide, thus they are then stored in the set $F_2$. With the hash function $h_2$, the keys from $F_2$ have no collision, and the process finishes. The MPHF query operation is very similar to the construction algorithm. Let $A$ be the concatenation of $A_0, A_1, A_2$ (see bottom part of the figure). To query $k_2$, the key is first hashed with $h_0$. The associated value in $A_0$ is '0', so $k_2$ is then hashed with $h_1$. The value associated in $A_1$ is again '0'. When finally hashed with $h_2$, the value associated in $A_2$ is '1' and thus the query stops here. The index returned by the MPHF is the rank of this '1' (here, 5) in $A$. In this example, the MPHF values returned when querying $k_1, k_2, k_3, k_4, k_5$ and $k_6$ are respectively 4,5,2,6,3, and 1.</figcaption></figure>

size $|F_0|$ is created such that there is a 1 at position $i$ if and only if exactly one element of $F_0$ has a hash value of $i$. We say that there is a *collision* whenever two keys in $F_0$ have the same hash value. Keys from $F_0$ that were involved in a collision are inserted into a new set $F_1$. The process repeats with $F_1$ and a new hash function $h_1$. A new bit array $A_1$ of size $|F_1|$ is created using the same procedure as for $A_0$ (except that $F_1$ is used instead of $F_0$, and $h_1$ instead of $h_0$). The process is repeated with $F_2, F_3, \dots$ until one of these sets, $F_{\text{last}+1}$, is empty.

We obtain an MPHF by concatenating the bit arrays $A_0, A_1, \dots, A_{\text{last}}$ into an array $A$. To perform a query, a key is hashed successively with hash functions $h_0, h_1, \dots$ as long as the value in $A_i$ ($i \ge 0$) at the position given by the hash function $h_i$ is 0. Eventually, by construction, we reach a 1 at some position of $A$ for some $i = d$. We say that the *level* of the key is $d$. The index returned by the MPHF is the rank of this one in $A$. See Figure 1 for an example.

## 2.2 Algorithm details

### 2.2.1 Collision detection

During construction at each level $d$, collisions are detected using a temporary bit array $C_d$ of size $|A_d|$. Initially all $C_d$ bits are set to '0'. A bit of $C_d[i]$ is set to '1' if two or more keys from
XX:4 Fast and scalable minimal perfect hashing for massive key sets

$F_d$ have the same value $i$ given by hash function $h_d$. Finally, if $C_d[i] = 1$, then $A_d[i] = 0$. Formally:

$$C_d[i] = 1 \Rightarrow A_d[i] = 0;$$

$$(h_d[x] = i \text{ and } A_d[i] = 0 \text{ and } C_d[i] = 0) \Rightarrow A_d[i] = 1 \text{ (and } C_d[i] = 0);$$

$$(h_d[x] = i \text{ and } A_d[i] = 1 \text{ and } C_d[i] = 0) \Rightarrow A_d[i] = 0 \text{ and } C_d[i] = 1.$$

### 2.2.2 Queries

A query of a key $x$ is performed by finding the smallest $d$ such that $A_d[h_d(x)] = 1$. The (non minimal) hash value of $x$ is then $(\sum_{i<d} |F_i|) + h_d(x)$.

### 2.2.3 Minimality

To ensure that the image range of the function is $[1, |F_0|]$, we compute the cumulative rank of each '1' in the bit arrays $A_i$. Suppose, that $d$ is the smallest value such that $A_d[h_d(x)] = 1$. The minimal perfect hash value is given by $\sum_{i<d}(\text{weight}(A_i) + \text{rank}(A_d[h_d(x)]), where \text{weight}(A_i)$ is the number of bits set to '1' in the $A_i$ array, and $\text{rank}(A_d[y])$ is the number of bits set to 1 in $A_d$ within the interval $[0, y]$, thus $\text{rank}(A_d[y]) = \sum_{j<y} A_d[j]$. This is a classic method also used in other MPHFs [3].

### 2.2.4 Faster query and construction times (parameter $\gamma$)

The running time of the construction depends on the number of collisions on the $A_d$ arrays, at each level $d$. One way to reduce the number of collisions, hence to place more keys at each level, is to use bit arrays ($A_d$ and $C_d$) larger than $|F_d|$. We introduce a parameter $\gamma \in \mathbb{R}$, $\gamma \ge 1$, such that $|C_d| = |A_d| = \gamma|F_d|$. With $\gamma = 1$, the size of $A$ is minimal. With $\gamma \ge 2$, the number of collisions is significantly decreased and thus construction and query times are reduced, at the cost of a larger MPHF structure size. The influence of $\gamma$ is discussed in more detail in the following analyses and results.

## 2.3 Analysis

Proofs of the following observations and lemma are given in the Appendix.

### 2.3.1 Size of the MPHF

The expected size of the structure can be determined using a simple argument, previously made in [6]. When $\gamma = 1$, the expected number of keys which do not collide at level $d$ is $|A_d|e^{-1}$, thus $|A_d| = |A_{d-1}|(1 - e^{-1}) = |A_0|(1 - e^{-1})^d$. In total, the expected number of bits required by the hashing scheme is $\sum_{d \ge 0} |A_d| = N \sum_{d \ge 0} (1 - e^{-1})^d = eN$, with $N$ being the total number of input keys ($N = |F_0|$). Note that consequently the image of the hash function is also in $[1, eN]$, before minimization using the rank technique. When $\gamma \ge 1$, the expected proportion of keys without collisions at each level $d$ is $|A_d|e^{-\frac{1}{\gamma}}$. Since each $A_d$ no longer uses one bit per key but $\gamma$ bits per key, the expected total number of bits required by the MPHF is $\gamma e^{\frac{1}{\gamma}}N$.

### 2.3.2 Space usage during construction

We analyze the disk space used during construction. Recall that during construction of level $d$, a bit array $C_d$ of size $|A_d|$ is used to record collisions. Note that the $C_d$ array is only
needed during the *d*-th level. It is deleted before level *d* + 1. The total memory required during level *d* is $\sum_{i \le d} (|A_i|) + |C_d| = \sum_{i < d} (|A_i|) + 2|A_d|$.

**► Lemma 1.** For $\gamma > 0$, the space of our MPHF is $S = \gamma e^{\frac{1}{\gamma}} N$ bits. The maximal space during construction is $S$ when $\gamma \le \log(2)^{-1}$, and $2S$ bits otherwise.

A full proof of the Lemma is provided in the Appendix.

# 3 Implementation

We present *BBhash*, a C++ implementation available at http://github.com/rizkg/BBHash.
We describe in this section some design key choices and optimizations.

## 3.1 Rank structure

We use a classical technique to implement the rank operation: the ranks of a fraction of the '1's present in $A$ are recorded, and the ranks in-between are computed dynamically using the recorded ranks as checkpoints.

In practice 64 bit integers are used for counters, which is enough for realistic use of an MPHF, and placed every 512 positions by default. These values were chosen as they offer a good speed/memory trade-off, increasing the size of the MPHF by a factor 1.125 while achieving good query performance. The total size of the MPHF is thus $(1 + \frac{64}{512})\gamma e^{\frac{1}{\gamma}} N$.

## 3.2 Parallelization

Parallelization is achieved by partitioning keys over several threads. The algorithm presented in Section 2 is executed on multiple threads concurrently, over the same memory space. Built-in compiler functions (e.g. `sync_fetch_and_or`) are used for concurrent access in the $A_i$ arrays. The efficiency of this parallelization scheme is shown in the Results section, but note that it is fundamentally limited by random memory accesses to the $A_i$ arrays which incur cache misses.

## 3.3 Hash functions

The MPHF construction requires classical hash functions. Other authors have observed that common hash functions behave practically as well as fully random hash functions [2]. We therefore choose to use xor-shift based hash functions [13] for their efficiency both in terms of computation speed and distribution uniformity [15].

## 3.4 Disk usage

In the applications we consider, key sets are typically too big to fit in RAM. Thus we propose to read them on the fly from disk. There are mainly two distinct strategies regarding the disk usage during construction: 1/ during each level *d*, keys that are to be inserted in the set $F_{d+1}$ are written directly to disk. The set $F_{d+1}$ is then read during level *d* + 1 and erased before level *d* + 2; or 2/ at each level all keys from the original input key file are read and queried in order to determine which keys were already assigned to a level *i* < *d*, and which would belong to $F_d$. When the key set becomes small enough (below user-defined threshold) it is loaded in ram to avoid costly re-computation from scratch at each level.

The first strategy obviously provides faster construction at the cost of temporary disk usage. At each level *d* > 0, two temporary key files are stored on disk: $F_d$ and $F_{d+1}$. The highest disk
XX:6 Fast and scalable minimal perfect hashing for massive key sets

usage is thus achieved during level 1, i.e. by storing $|F_1| + |F_2| = |F_0|((1-e^{-1/\gamma})+(1-e^{-1/\gamma})^2)$ elements. With $\gamma = 1$, this represents $\approx 1.03N$ elements, thus the construction overhead on disk is approximately the size of the input key file. Note that with $\gamma = 2$ (resp. $\gamma = 5$), this overhead diminishes and becomes a ratio of $\approx 0.55$ (resp. $\approx 0.21$) the size of the input key file.

The first strategy is the default strategy proposed in our implementation. The second one has also been implemented and can be optionally switched on.

## 3.5 Termination

The expected number of unplaced keys decreases exponentially with the number of levels but is not theoretically guaranteed to reach zero in a finite number of steps. To ensure termination of the construction algorithm, in our implementation a maximal number $D$ of levels is fixed. Then, the remaining keys are inserted into a regular hash table. Value $D$ is a parameter, its default value is $D = 25$ for which the expected number of keys stored in this hash table is $\approx 10^{-5}N$ for $\gamma = 1$ and becomes in practice negligible for $\gamma \ge 2$, allowing the size overhead of the final hash table to be negligible regarding the final MPHF size.

# 4 Results

We evaluated the performance *BBhash* for the construction of large *MPHFs*. We generated files containing various numbers of keys (from 1 million to 1 trillion keys). In our tests, a key is a binary representation of a pseudo-random positive integer in $[0; 2^{64}]$. Within each file, each key is unique. We also performed a test where input keys are strings (n-grams) to ensure that using integers as keys does not bias results. Tests were performed on a cluster node with a Intel® Xeon® CPU E5-2660 v3 2.60GHz 20-core CPU, 256 GB of memory, and a mechanical hard drive. Except for the experiment with $10^{12}$ keys, running times include the time needed to read input keys from disk. Note that files containing key sets may be cached in memory by the operating system, and all evaluated methods benefit from this effect during MPHF construction. We refer to the Appendix for the specific commands and parameters used in these experiments.

We first analyzed the influence of the $\gamma$ value (the main parameter of *BBhash*), then the effect of using multiple threads depending on the parallelization strategy. Second, we compared *BBhash* with other state-of-the-art methods. Finally, we performed an MPHF construction on $10^{12}$ elements.

## 4.1 Influence of the $\gamma$ parameter

We report in Figure 2 (left) the construction times and the mean query times, as well as the size of the produced MPHF, with respect to several $\gamma$ values. The main observation is that $\gamma \ge 2$ drastically accelerates construction and query times. This is expected since large $\gamma$ values allow more elements to be placed in the first levels of the MPHF; thus limiting the number of times each key is hashed to compute its level. In particular, for keys placed in the very first level, the query time is limited to a single hashing and a memory access. The average level of all keys is $e^{(1/\gamma)}$, we therefore expect construction and query times to decrease when $\gamma$ increases. However, larger $\gamma$ values also incur larger MPHF sizes. One observes that $\gamma > 5$ values seem to bring very little advantage at the price of higher space requirements. A related work used $\gamma = 1$ in order to minimize the MPHF size [6]. Here, we
<figure><img src="image_44.png" /><figcaption>Figure 2 Left: Effects of the gamma parameter on the performance of *BBhash* when run on a set composed of one billion keys, when executed on a single CPU thread. Times and MPHF size behave accordingly to the theoretical analysis, respectively $O(e^{(1/\gamma)})$, and $O(\gamma e^{(1/\gamma)})$. Right: Performance of the *BBhash* construction time according to the number of cores, using $\gamma = 2$.</figcaption></figure>

argue that using $\gamma$ values larger than 1 has significant practical merits. In our tests, we often used $\gamma = 2$ as it yields an attractive time/space trade-off during construction and queries.

## 4.2 Parallelization performance

We evaluated the capability of our implementation to make use of multiple CPU cores. In Figure 2 (right), we report the construction times with respect to the number of threads. We observe a near-ideal speed-up with respect to the number of threads with diminishing returns when using more than 10 threads, which is likely due to cache misses that induce a memory access bottleneck.

In addition to these results, we applied *BBhash* on a key set of 10 billion keys and on a key set of 100 billion keys, again using default parameters and 8 threads. The memory usage was respectively 4.96GB and 49.49GB, and the construction time was respectively 462 seconds and 8913 seconds, showing the scalability of *BBhash*.

## 4.3 Comparisons with state of the art methods

We compared *BBhash* with state-of-the-art MPHF methods. CHD (http://cmph.sourceforge.net/) is an implementation of the compressed hash-and-displace algorithm [2]. EMPHF [1] is based on random hypergraph peeling, and the HEM [4] implementation in EMPHF is based on partitioning the input data. Sux4J is a Java implementation of [11]. We did not include other methods cited earlier because they do not provide an implementation [12, 16] or the software integrates a non-minimal perfect hash function that is not stand-alone [6]. However single-threaded results presented in [16] show that construction times and MPHF sizes are comparable to ours, query times are significantly longer, and no indication is provided about the memory usage during construction. Our benchmark code is available at https://github.com/rchikhi/benchmphf.

Figure 3 shows that all evaluated methods are able to construct MPHFs that contain a billion elements, but only *BBhash* scales up to datasets that contain $10^{11}$ elements and more. Overall, *BBhash* shows consistently better time and memory usage during construction.

We additionally compared the resulting MPHF size, i.e. the space of the data structure
XX:8 Fast and scalable minimal perfect hashing for massive key sets

<figure><img src="image_45.png" /><figcaption>■ **Figure 3** Memory footprint and construction time with respect to the number of keys. All libraries were run using default parameters, including $\gamma = 2$ for *BBhash*. For a fair comparison, *BBhash* was executed on a single CPU thread. Except for Sux4J, missing data points correspond to runs that exceeded the amount of available RAM. Sux4J limit comes from the disk usage, estimated at approximately 4TB for $10^{11}$ keys.</figcaption></figure>

returned by the construction algorithm, and the mean query time across all libraries on a dataset consisting of a billion keys (Table 1). MPHFs produced by *BBhash* range from 2.89 bits/key (when $\gamma = 1$ and ranks are sampled every 1024 positions) to 6.9 bits/key (when $\gamma = 5$ and a rank sampling of 512). The 0-0.8 bits/key size difference between our implementation and the theoretical space usage of the *BBhash* structure size is due to additional space used by the rank structure. We believe that a reasonable compromise in terms of query time and structure size is 3.7 bits/key with $\gamma = 2$ and a rank sampling of 512, which is marginally larger than the MPHF sizes of other libraries (ranging from 2.6 to 3.5 bits/key). As we argued in the Introduction, using one more bit per key seems to be a reasonable trade-off for performance.

Construction times vary by one or two orders of magnitude across methods, *BBhash* being the fastest. With default parameters ($\gamma = 2$, rank sampling of 512), *BBhash* has a construction memory footprint 40× to 60× smaller than other libraries except for Sux4j, for which *BBhash* remains 4× smaller. Query times are roughly within an order of magnitude (179 – 1037 ns) of each other across methods, with a slight advantage for *BBhash* when $\gamma \ge 2$. Sux4j achieves an attractive balance with low construction memory and query times, but high disk usage. In our tests, the high disk usage of Sux4j was a limiting factor for the construction of very large MPHFs.
<table><thead><tr><th>Method</th><th>Query time (ns)</th><th>MPHF size (bits/key)</th><th>Const. time* (s)</th><th>Const. memory**</th><th>Disk. usage (GB)</th></tr></thead><tbody><tr><td>BBhash γ = 1</td><td>271</td><td>3.1</td><td>60 (393)</td><td>3.2 (376)</td><td>8.23</td></tr><tr><td>BBhash γ = 1 minirank</td><td>279</td><td>2.9</td><td>61(401)</td><td>3.2 (376)</td><td>8.23</td></tr><tr><td>BBhash γ = 2</td><td>216</td><td>3.7</td><td>35 (229)</td><td>4.3 (516)</td><td>4.45</td></tr><tr><td>BBhash γ = 2 nodisk</td><td>216</td><td>3.7</td><td>80 (549)</td><td>6.2 (743)</td><td>0</td></tr><tr><td>BBhash γ = 5</td><td>179</td><td>6.9</td><td>25 (162)</td><td>10.7 (1,276)</td><td>1.52</td></tr><tr><td>EMPHF</td><td>246</td><td>2.9</td><td>2,642</td><td>247.1 (29,461)†</td><td>20.8</td></tr><tr><td>EMPHF HEM</td><td>581</td><td>3.5</td><td>489</td><td>258.4 (30,798)†</td><td>22.5</td></tr><tr><td>CHD</td><td>1037</td><td>2.6</td><td>1,146</td><td>176.0 (20,982)</td><td>0</td></tr><tr><td>Sux4J</td><td>252</td><td>3.3</td><td>1,418</td><td>18.10 (2,158)</td><td>40.1</td></tr></tbody></table>

■ Table 1 Performance of different MPHF algorithms applied on a key set composed of 10⁹ 64-bits random integers, of size 8GB. Each time result is the average value over three tests. The 'nodisk' row implements the second strategy described in Section 3.4, and the 'minirank' row samples ranks every 1024 positions instead of 512 by default. *The column "Const. time" indicates the construction time in seconds. In the case of BBhash, the first value is the construction time using eight CPU threads and the second value in parenthesis is the one using one CPU thread. **The column "Const. memory" indicates the RAM used during the MPHF construction, in bits/key and the total in MB in parenthesis. † The memory usages of EMPHF and EMPHF HEM reflect the use of memory-mapped files (mmap scheme).

Note that EMPHF, EMPHF HEM and Sux4j implement a disk partitioning strategy, that could in principle also be applied to others methods, including ours. Instead of creating a single large MPHF, they partition the set of input keys on disk and construct many small MPHFs independently. In theory this technique allows to engineer the MPHF construction algorithm to use parallelism and lower memory, at the expense of higher disk usage. In practice we observe that the existing implementations that use this technique are not parallelized. While EMPHF and EMPHF HEM used relatively high memory in our tests (around 30 GB for 1 billion elements) due to memory-mapped files, they also completed the construction successfully on another machine that had 16 GB of available memory. However, we observed what appears to be limitations in the scalability of the scheme: we were unable to run EMPHF and EMPHF HEM on an input of 10 billion elements using 256 GB of memory. Regardless, we view this partitioning technique as promising but orthogonal to the design of efficient "monolithic" MPHFs constructions such as *BBhash*.

## 4.4 Performance on an actual dataset

In order to ensure that using pseudo-random integers as keys does not bias results, we ran *BBhash* using strings as keys. We used n-grams extracted from the Google Books Ngram dataset<fnref n="1" />, version 20120701. On average the n-gram size is 18. We also generated random words of size 18. As reported in Table 2, we obtained highly similar results to those obtained with random integer keys.

<footnote n="1">¹ http://storage.googleapis.com/books/ngrams/books/datasetsv2.html</footnote>
XX:10 Fast and scalable minimal perfect hashing for massive key sets

<table><caption>Table 2 Performance of *BBhash* ($\gamma$ = 2, 8 threads) when using ASCII strings as keys.</caption><thead><tr><th>Dataset</th><th>Query time (ns)</th><th>MPHF size (bits/key)</th><th>Const. time (s)</th></tr></thead><tbody><tr><td>10<sup>8</sup> Random strings</td><td>325</td><td>3.7</td><td>35</td></tr><tr><td>10<sup>8</sup> Ngrams</td><td>296</td><td>3.7</td><td>37</td></tr></tbody></table>

## 4.5 Indexing a trillion keys

We performed a very large-scale test by creating an MPHF for $10^{12}$ keys. For this experiment, we used a machine with 750 GB of RAM. Since storing that many keys would require 8 TB of disk space, we instead used a procedure that deterministically generates a stream of $10^{12}$ pseudo-random integers in $[0, 2^{64} - 1]$. We considered the streamed values as input keys without writing them to disk. In addition, key sets of cardinality below 20 billion (2% of the input) were stored in memory to avoid re-computation from scratch at each subsequent level. Thus, the reported computation time should not be compared to previously presented results as this experiment has no disk accesses. The test was performed using $\gamma$ = 2, 24 threads.

Creating the MPHF took 35.4 hours and required 637 GB RAM. This memory footprint is roughly separated between the bit arrays ($\approx$ 459 GB) and the memory required for loading 20 billion keys in memory ($\approx$ 178 GB). The final MPHF occupied 3.71 bits per key.

## 5 Conclusion

We have proposed a resource-efficient and highly scalable algorithm for constructing and querying MPHFs. Our algorithmic choices were motivated by simplicity: the method only relies on bit arrays and classical hash functions. While the idea of recording collisions in bit arrays to create MPHFs is not novel [6, 12], to the best of our knowledge *BBhash* is the first implementation that is competitive with the state of the art. The construction is particularly time-efficient as it is parallelized and mainly consists in hashing keys and performing memory accesses. Moreover, the additional data structures used during construction are provably small enough to ensure a low memory overhead during construction. In other words, creating the MPHF does not require much more space than the resulting MPHF itself. This aspect is important when constructing MPHFs on large key sets in practice.

Experimental results show that *BBhash* generates MPHFs that are slightly larger to those produced by other methods. However *BBhash* is by far the most efficient in terms of construction time, query time, memory and disk footprint for indexing large key sets (of cardinality above $10^9$ keys). The scalability of our approach was confirmed by constructing MPHFs for sets as large as $10^{12}$ keys. To the best of our knowledge, no other MPHF implementation has been tested on that many keys.

A time/space trade-off is achieved through the $\gamma$ parameter. The value $\gamma$ = 1 yields MPHFs that occupy roughly $3N$ bits of space and have little memory overhead during construction. Higher $\gamma$ values use more space for the construction and the final structure size, but they achieve faster construction and query times. Our results suggest that $\gamma$ = 2 is a good time-versus-space compromise, using 3.7 bits per key. With respect to hypergraph-based methods [1, 3, 4, 11], *BBhash* offers significantly better construction performance, but the resulting MPHF size is up to 1 bit/key larger. We however argue that the MPHF size, as long as it is limited to a few bits per key, is generally not a bottleneck as many applications use MPHFs to associate much larger values to keys. Thus, we believe that this work will unlock many high performance computing applications where the possibility to index billions keys
and more is a huge step forward.

An interesting direction for future work is to obtain more space-efficient MPHFs using our method. We believe that a way to achieve this goal is to slightly change the hashing scheme. We would like to explore an idea inspired by the CHD algorithm for testing several hash functions at each level and selecting (then storing) one that minimizes the number of collisions. At the price of longer construction times, we anticipate that this approach could significantly decrease the final structure size.

## Acknowledgments

This work was funded by French ANR-12-BS02-0008 Colib'read project. We thank the GenOuest BioInformatics Platform that provided the computing resources necessary for benchmarking. We thank Djamal Belazzougui for helpful discussions and pointers.

## References

1. Djamal Belazzougui, Paolo Boldi, Giuseppe Ottaviano, Rossano Venturini, and Sebastiano Vigna. Cache-oblivious peeling of random hypergraphs. In *Data Compression Conference (DCC), 2014*, pages 352–361. IEEE, 2014.

2. Djamal Belazzougui, Fabiano C Botelho, and Martin Dietzfelbinger. Hash, displace, and compress. In *European Symposium on Algorithms*, pages 682–693. Springer, 2009.

3. Fabiano C Botelho, Rasmus Pagh, and Nivio Ziviani. Simple and space-efficient minimal perfect hash functions. In *Algorithms and Data Structures*, pages 139–150. Springer, 2007.

4. Fabiano C Botelho, Rasmus Pagh, and Nivio Ziviani. Practical perfect hashing in nearly optimal space. *Information Systems*, 38(1):108–131, 2013.

5. Chin-Chen Chang and Chih-Yang Lin. Perfect hashing schemes for mining association rules. *48(2):168–179, 2005*. doi:10.1093/comjnl/bxh074.

6. Jarrod A Chapman, Isaac Ho, Sirisha Sunkara, Shujun Luo, Gary P Schroth, and Daniel S Rokhsar. Meraculous: de novo genome assembly with short paired-end reads. *PloS one*, 6(8):e23501, 2011.

7. Yupeng Chen, Bertil Schmidt, and Douglas L Maskell. A hybrid short read mapping accelerator. *BMC Bioinformatics*, (1):67. doi:10.1186/1471-2105-14-67.

8. Rayan Chikhi, Antoine Limasset, and Paul Medvedev. Compacting de bruijn graphs from sequencing data quickly and in low memory. *Bioinformatics*, 32(12):i201–i208, 2016.

9. Zbigniew J Czech, George Havas, and Bohdan S Majewski. Perfect hashing. *Theoretical Computer Science*, 182(1):1–143, 1997.

10. Michael L Fredman and János Komlós. On the size of separating systems and families of perfect hash functions. *SIAM Journal on Algebraic Discrete Methods*, 5(1):61–68, 1984.

11. Marco Genuzio, Giuseppe Ottaviano, and Sebastiano Vigna. Fast scalable construction of (minimal perfect hash) functions. In V. Andrew Goldberg and S. Alexander Kulikov, editors, *Experimental Algorithms: 15th International Symposium, SEA 2016, St. Petersburg, Russia, June 5-8, 2016, Proceedings*, pages 339–352. Springer International Publishing, Cham, 2016. doi:10.1007/978-3-319-38851-9_23.

12. Yi Lu, Balaji Prabhakar, and Flavio Bonomi. Perfect hashing for network applications. In *2006 IEEE International Symposium on Information Theory*, pages 2774–2778. IEEE, 2006.

13. George Marsaglia et al. Xorshift rings. *Journal of Statistical Software*, 8(14):1–6, 2003.

14. Kurt Mehlhorn. On the program size of perfect and universal hash functions. In *Foundations of Computer Science, 1982. SFCS'08. 23rd Annual Symposium on*, pages 170–175. IEEE, 1982.
15 Michael Mitzenmacher and Salil Vadhan. Why simple hash functions work: exploiting the entropy in a data stream. In *Proceedings of the nineteenth annual ACM-SIAM symposium on Discrete algorithms*, pages 746–755. Society for Industrial and Applied Mathematics, 2008.

16 Ingo Müller, Peter Sanders, Robert Schulze, and Wei Zhou. *Retrieval and Perfect Hashing Using Fingerprinting*, pages 138–149. Springer International Publishing, Cham, 2014. doi: 10.1007/978-3-319-07959-2_12.
# Appendix

## Proofs of MPHF size and memory required for construction

MPHF size with $\gamma = 1$.

$$ 
\begin{aligned}
\sum_{d \ge 0} |A_d| &= N \sum_{d \ge 0} (1 - e^{-1})^d \\
&= N \frac{1}{1 - (1 - e^{-1})} && \text{as } \lim_{d \to +\infty} (1 - e^{-1})^d = 0 \\
&= eN
\end{aligned}
$$

MPHF size using any $\gamma \ge 1$. With $\gamma \ge 1$: $|A_d| = \gamma|A_{d-1}|(1 - e^{-1/\gamma}) = \gamma|A_0|(1 - e^{-1/\gamma})^d = \gamma N(1 - e^{-1/\gamma})^d$

$$ \text{Thus, } \sum_{d \ge 0} |A_d| = \gamma N \sum_{d \ge 0} (1 - e^{-1/\gamma})^d $$

Moreover, as $\lim_{d \to +\infty} (1 - e^{-1/\gamma})^d = 0$ since for $\gamma > 0, 0 < 1 - e^{-1/\gamma} < 1$, on has:

$$ \sum_{d \ge 0} |A_d| = \gamma N \frac{1}{1 - (1 - e^{-1/\gamma})} = \gamma e^{\frac{1}{\gamma}} N $$

Note that this proof stands for any $\gamma$ value $> 0$, but that with $\gamma < 1$ the theoretical and practical MPHF sizes increase exponentially as $\gamma$ get close to zero.

**Lemma 1.** Let $m(d)$ be memory required during level $d$ and let $R$ be the ratio between the maximal memory needed during the MPHF construction and the MPHF total size denoted by $S$. Formally,

$$ R = \frac{\max_{d \ge 0}(m(d))}{S} = \frac{\max_{d \ge 0}(m(d))}{\gamma e^{\frac{1}{\gamma}} N} $$

First we prove that $\lim_{d \to \infty} \frac{m(d)}{S} = 1$.

$$ m(d) = \sum_{i<d} |A_i| + 2|A_d| = \gamma N \left( \frac{1 - (1 - e^{-1/\gamma})^d}{e^{-1/\gamma}} + 2(1 - e^{-1/\gamma})^d \right) $$

Since for $\gamma > 0, 0 < 1 - e^{-1/\gamma} < 1$, then $\lim_{d \to \infty} m(d) = \gamma e^{\frac{1}{\gamma}} N$. Thus $\lim_{d \to \infty} \frac{m(d)}{S} = 1$. Before going further, we need to compute $m(d+1) - m(d)$:

$$ 
\begin{aligned}
m(d+1) - m(d) &= \sum_{i<d+1} |A_i| + 2|A_{d+1}| - \sum_{i<d} |A_i| + 2|A_d| \\
&= |A_d| + 2|A_{d+1}| - 2|A_d| = 2|A_{d+1}| - |A_d| \\
&= 2\gamma N(1 - e^{-1/\gamma})^{d+1} - \gamma N(1 - e^{-1/\gamma})^d \\
&= \gamma N(1 - e^{-1/\gamma})^d(2(1 - e^{-1/\gamma}) - 1) \\
&= \gamma N(1 - e^{-1/\gamma})^d(1 - 2e^{-1/\gamma})
\end{aligned}
$$

We now prove $R \le 1$ when $\gamma \le \frac{1}{\log(2)}$ and also, $R < 2$ when $\gamma > \frac{1}{\log(2)}$.
XX:14 Fast and scalable minimal perfect hashing for massive key sets

■ Case 1: $\gamma \le \frac{1}{\log(2)}$

We have $\frac{m(0)}{S} = 2e^{-\frac{1}{\gamma}} \le 2e^{-\log(2)} = 1$.

Moreover, as $m(d+1) - m(d) = \gamma N(1 - e^{-\frac{1}{\gamma}})^d (1 - 2e^{-\frac{1}{\gamma}})$ and as, with $\gamma \le \frac{1}{\log(2)}$:
$1 - e^{-\frac{1}{\gamma}} \ge 0.5$, and $1 - 2e^{-\frac{1}{\gamma}} \ge 0$ then $m(d+1) - m(d) \ge 0$, thus, $m$ is an increasing function.

To sum up, with $\gamma \le \frac{1}{\log(2)}$, we have **1/** that $\frac{m(0)}{S} \le 1$, **2/** that $\lim_{d \to \infty} \frac{m(d)}{S} = 1$, and **3/** that $m$ is increasing, then $R \le 1$.

■ Case 2: $\gamma > \frac{1}{\log(2)}$ We have $\frac{m(0)}{S} = 2e^{-\frac{1}{\gamma}}$. With $\gamma > \frac{1}{\log(2)}$, $1 < \frac{m(0)}{S} < 2$. Moreover, $m(d+1) - m(d) = \gamma N(1 - e^{-\frac{1}{\gamma}})^d (1 - 2e^{-\frac{1}{\gamma}})$ is negative as: $1 - e^{-\frac{1}{\gamma}} > 0$ and $1 - 2e^{-\frac{1}{\gamma}} < 0$ for $\gamma > \frac{1}{\log(2)}$. Thus $m$ is a decreasing function with $d$.

With $\gamma > \frac{1}{\log(2)}$, we have **1/** that $\frac{m(0)}{S} < 2$, **2** that $\lim_{d \to \infty} \frac{m(d)}{S} = 1$ and **3** that $m$ is decreasing. Thus $R < 2$.

<figure><img src="image_46.png" /></figure>
Algorithms pseudo-codes

<caption>**Algorithm 1:** MPHF construction.</caption>

**Data:** $F_0$ a set of $N$ keys, integers $\gamma$ and $last$
**Result:** array of bit arrays $\{A_0, A_1, \dots, A_{last}\}$, hash table $H$
$i=0$;
**while** $F_i$ not empty and $i \le last$ **do**
| $A_i$ = ArrayFill($F_i$, $\gamma$);
| **foreach** key $x$ of $F_i$ **do**
| | $h$ = hash($x$) mod $(\gamma * N)$;
| | **if** $A_i[h] == 0$ **then**
| | | $F_{i+1}$.add($x$)
| | $i=i+1$;
| **Construct** $H$ using remaining elements from $F_{last+1}$;
| **Return** $\{A_0, A_1, \dots, A_{last}, H\}$

In practice $F_i$ with $i > 1$ are stored on disk (see Section 3.4). The hash table $H$ ensures that elements in $F_{last+1}$ are mapped without collisions to integers in $[[F_0] - |F_{last+1}| + 1, |F_0|]$

<caption>**Algorithm 2:** ArrayFill</caption>

**Data:** $F$ array of $N$ keys, integer $\gamma$
**Result:** bit array $A$
Zero-initialize $A$ and $C$ two bit arrays with $\gamma * N$ elements;
**foreach** key $x$ of $F$ **do**
| $h$ = hash($x$) mod $(\gamma * N)$;
| **if** $A[h] == 0$ and $C[h] == 0$ **then**
| | $A[h] = 1$;
| **if** $A[h] == 1$ and $C[h] == 0$ **then**
| | $A[h] = 0$;
| | $C[h] = 1$;
| **if** $A[h] == 0$ and $C[h] == 1$ **then**
| | Skip;
| Delete $C$;
| **Return** $A$;

Note that the case $A[h] == 1$ and $C[h] == 1$ never happens.

<caption>**Algorithm 3:** MPHF query</caption>

**Data:** bit arrays $\{A_0, A_1, \dots, A_{last}\}$, hash table $H$, key $x$
**Result:** integer index of $x$
$i=0$;
**while** $i \le last$ **do**
| $h$ = hash$_i(x)$ mod $A_i$.size();
| **if** $A_i[h] == 1$ **then**
| | **return** $\sum_{j<i} |A_j| + \text{rank}(A_i[h])$;
| $i = i + 1$;
| **return** $H[x]$;

Note, when $x$ is not an element from the key set of the MPHF, the algorithm may return a wrong integer index.
XX:16 Fast and scalable minimal perfect hashing for massive key sets

Commands

In this section we describe used commands for each presented result. Time and memory
usages where computed using using "`/usr/bin/time -verbatim`" unix command. The disk usage
was computed thanks to a home made script measuring each 1/10 second the size of the
directory using the "`du -sk`" unix command, and recording the highest value. The `BBhash`
library and its `Bootest` tool are available from https://github.com/rizkg/BBHash.

Commands used for Section 4.1:

```
for ((gamma=1;gamma<11;gamma++)); do
./Bootest 1000000000 1 ${gamma} -bench
done
```

Note that 1000000000 is the number of keys tested and 1 is the number of used cores.
Additional tests, with larger key set and 8 threads:

```
for ((gamma=1;gamma<11;gamma++)); do
./Bootest 1000000000 1 ${gamma} -bench
done
```

Commands used for Section 4.2:

```
for keys in 1000000000 10000000000; do
./Bootest ${keys} 8 2 -bench
done
```

Commands used for Section 4.3:

We remind that our benchmark code, testing EMPHF, EMPHF MEM, CHD, and Sux4J is
available at https://github.com/rchikhi/benchmphf.

BBhash commands:

```
for keys in 1000000 10000000 100000000 10000000000\  
1000000000 100000000000; do
./Bootest ${keys} 1 2 -bench
done
```

BBhash command with nodisk (Table 1) was

```
./Bootest 1000000000 1 2 -bench -nodisk
```

and

```
./Bootest 1000000000 8 2 -bench -nodisk
```

Commands EMPHF & EMPHF HEM:

```
for keys in 1000000 10000000 100000000 10000000000\  
1000000000 100000000000; do
./benchmphf ${keys} -emphf
done
```

EMPHF (resp. EMPHF HEM) is tested by using the #define EMPHF_SCAN macro
ndon to assess the disk size footprint, the line
"`unlink(tmp);`" from file "`emphf/mmap_memory_model.hpp`" was commented.

Commands CHD:
for keys in 1000000 10000000 100000000 10000000000\ 1000000000 100000000000; do ./benchmphf ${keys} -chd done

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

Commands Sux4J:
for each size, the "Sux4J/slow/it/unimi/dsi/sux4j/mph/LargeLongCollection.java" was
modified indicating the used size.

./run-sux4j-mphf.sh

Commands used for Section 4.4:

As explained Section 4.4, the keyString.txt file is composed of n-grams extracted from
the Google Books Ngram dataset<fnref n="2" />, version 20120701.

./BootestFile keyStrings.txt 10 2

Commands used for Section 4.5:

BBhash command for indexing a trillion keys, with keys generated on the fly.

./Bootest 1000000000000 24 2 -onthefly

<footnote n="2">² http://storage.googleapis.com/books/ngrams/books/datasetsv2.html</footnote>
# Chapter 3

The de Bruijn graph as a reference
In this chapter we are interested in considering the de Bruijn graph as a reference. We first describe the different structures used to represent a genome and argue why a de Bruijn graph could be an interesting object to do so (Section 1). Then we present our theoretical and practical contribution about mapping reads on a de Bruijn graph (Section 2). In a last part, we show a direct application of such a method for read correction (Section 3).

# 3.1 Genome representations

## 3.1.1 Reference sequences

The current main representation of a genome is a set of sequences. Fasta files (Figure 3.1) seemed at first adapted to represent a genome. Due to the hardness of producing a genome of reasonable quality via the sequencing data, the loss of allelic information was considered secondary and non-impacting Slightly chimeric sequences were seen as a small cost to pay to access the global structure of the genome. The hardness of the haplotype separation, still present today, tends to show that it was a reasonable choice to postpone such a challenge. But there are several reasons why the set of sequences representation is not ideal and we detail them in the following paragraphs.

**Heterozygosity** Actual reference sequences are generated by assemblers that ignore the heterozygosity information. Sequences produced from several haplotypes contain nucleotides from different alleles (Figure 3.2). When the dissimilarity between alleles is high, or when they present structural variants [89] (such as a large insertion), the output assembly may present chimeric sequences.

We could argue that we could still use different sequences to represent the alleles. The challenge is to be able to link the differences together in order to produce the different allele sequences. This process is called phasing (Figure 3.2). For low heterozygosity rate (of the order of one per thousand for a human), variation may be quite distant. In order to phase the successive variants it would require information with very long range

>seq1
CATGCATCGATGCCATCGATCGATCTAGCTGACTGATCCATGCCTGACCGATCGTAGCTAGCTA
CTACTGTACGTGACTGCTAGTCATCGATCGCATGACTG
>seq2
CATGACCATCGTAGCTAGCGACGTAGCTACTGATCTAGCTGATCGTACGTAGCTGATCTAGCTG
ACTGCT
...

<figure><img src="image_48.png" /><figcaption>Figure 3.1: Example of fasta file. The first line in a FASTA file starts either with a "></figcaption></figure>
<figure><img src="image_49.png" /><figcaption>Figure 3.2: Polyploid assembly opposed to haploid assembly. Most assembler practice "haploid" assembly (top) and crush bubbles and output contigs that may contain sequences from different haplotypes. A polyploid assembly would consist to produce haplotypes consistent contigs (bottom).</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>Figure 3.3: Diploid assembly can result in a suit of bubbles. In this example a diploid genome is assembled, homozygous regions in gray are shared by the two haplotypes. Each difference between them create a bubble in the graph.</figcaption></figure>
(Figure 3.3). This fact explains why it appears almost impossible to phase genomes with reads of hundreds of bases.

Nevertheless we are able to study variants even if we are not able to reconstruct the haplotypes. Some variants can be detected from the assembly graph [90] and reported by the assembler. The question may be "how to represent the variant information across the genome?". To place each variant at a position on a reference may be extremely useful. However it presents a concerning allele bias [91] [92]. If some small differences may not be a problem and can be easily be positioned on the reference, larger structural variations are not straightforward to encode and may be completely absent in the reference [93]. Those points show that the linear representation is not a good support for the heterozygosity information yet.

**Fragmented genomes** Nowadays, most genomes are unfinished and are left at the state of contigs or scaffolds. While this is sufficient for some applications like reads mapping, one can argue that the linear representation loses some information in the case of unfinished assembly. During the assembly process, contigs and scaffolds are extended until an ambiguity on the way to extend is encountered. Even if there are several possibilities that justify the contigs not to be extended, reads mapping could exploit the fact that we know possible extensions for its alignment. Such information could lead to a more efficient mapping, especially around repeats usually responsible for broken contigs. Another problem is that usually, small contigs are not output by assemblers and neither used for scaffolds creation or as references to map on. Again, there are also very good reasons to do so, such as contigs shorter than reads occasioning multiple mapping problems. However, assembly may present "holes" when a complex and hard to assemble region appears. Such a potential source of bias tends to show that a set of sequences representing scaffolds or contigs is not satisfying.

**Redundancy** Even if we were able to perfectly assemble genomes, we would be interested for many biological applications to store several genomes of different individuals. Very wide projects such as the 1000 genomes [94] plan the sequencing of massive amounts of individuals to study genomic variations among them. The raw storing of the individual sequences is highly inefficient in both terms of storage and searching because of the high redundancy due to the similarity among individuals from a same species. Since the differences between individuals may be very low and the amount of data very large, highly efficient structures may be used in order to index or store the amount of data in an entropic efficient way [95, 96].

### 3.1.2 Genome graphs

Graph structure is a natural candidate to represent genome for several reasons. Several current representations of a genome as a graph vastly depend on the planned use, as no data structure prevails for common use yet [97] [98]. To index multiple genomes, first representation used annotations to locate variations with respect to a reference [99]. This approach presents several flaws as scalability issues, reference bias and impossibility to represent large structural variations. For pan-genomes [100] propose an ordering of the set of linear sequences, providing an efficient way to create pan-genome references for
reads mapping. In order to represent the variations among haplotypes, the variation
graph has been introduced, where haplotypes are encoded as a set of walks through a
sequence graph. In order to perform queries on such heavy data structures, the Posi-
tional BWT (PBWT) have been proposed [101] and improved in gPBWT [102] to better
represent the haplotypes added in the graph. Those structures are very similar to linear
sequences since they are supposed to represent finished genomes. Even if these sequences
represent a real "map" of a genome with a coordinate system, this is still a hard task
to obtain such sequences. Several uses can benefit from a non fully positioned reference:
quantification [103] [104], assembly, scaffolding and indexation of repeated sequences or
fragmented references [105]. Assembly graphs are natural candidates for this task as they
are literally conceived to represent fragmented data. In the following, we present methods
and applications relying on assembly graphs as references. Even if such approaches could
be based on any assembly graph, the following will focus on the use of a de Bruijn graph.

3.1.3 De Bruijn graphs as references

This part concentrates on the interest of the de Bruijn graph to represent genomes. We
will present the following arguments:

* Easiness to construct from sequencing data

• Efficient representation of repeated data

* Represents variants and complex regions

**Construct a reference de Bruijn graph** To use a de Bruijn graph as a reference, it has to contain the whole genomic information with almost no sequencing error. Building a de Bruijn graph from a reference genome is easy but has a limited interest since the transformation in a de Bruijn graph is lossy as sequences respective positions may be lost. To represent fragmented references such as contigs, we can argue that a graph may contain more information than the contigs set itself. The real challenge is to produce a good reference from raw sequencing data. We present here different strategies to produce a good reference de Bruijn graph and evaluate them. In order to do this, we simulated reads from a known reference and compare the set of kmers from the reference and from the constructed graph. Two kinds of errors may appear:

* Some kmers are in the graph but not in the reference. They are erroneous kmers, also called "false positives"

* Some kmers are in the reference but not in the graph. They are missed kmers also called "false negatives"

As described in the introduction, the main technique to get rid of erroneous kmers is to
remove low abundance ones. Three zones can be observed in a kmer spectrum:

* A very high peak of low abundance kmers

* A flat zone

* A bell containing high abundance kmers

The essential idea is that genomic kmers will mostly be present in the third zone while
most erroneous kmer will be present in the first zone. By applying a coverage cutoff we
only keep the right part of the spectrum curves. The question is how to decide the cutoff
parameter. Techniques exist to find the end of the first peak by trying to fit models
on the kmer spectrum curve such as kmerGenie [34]. A deeper sequencing makes easier
<figure><img src="image_51.png" /><figcaption>Figure 3.4: Spectrum of 31mer from simulated reads from de Bruijn graph reference (Left). Reads were of length 100 and with coverages of 30 (blue), 50 (green) and 100 (purple) respectively. We observe that higher coverage results in a larger bell zone more distant from the error peak. It is thus easier to separate the erroneous kmers from others with a high coverage. On the right a kmer spectrum of a real sequencing with *k* = 31 and a coverage around 100X from *E. coli*. We can see that unlike spectrum from simulated dataset, the "hole" between the peak and the bell do not reach a value of zero since simulated dataset may not represent all aspect of NGS sequences distribution.</figcaption></figure>

the determination of such a cutoff, as can be seen in Figure 3.4 at left. We note that on simulated datasets, the determination of the cutoff can seem easy because the flat zone has no element. In practice on real dataset (in Figure 3.4 at right), the situation is not this simple since a notable number of element are present with the unexpected medium abundances. The value of *k* also has a high influence on the kmer spectrum. In Figure 3.5 we observe that the larger the *k*, the lesser the genomic kmers are covered and the closer to zone 1 and 3. If *k* is too high compared to the read length it may be difficult to separate the genomic and erroneous kmers since some genomic kmers will be present with a very low abundance. This is mainly due to the fact that when *k* comes close to the read length, very few kmers are produced by each read. To reduce this problem with large kmers, a step of read correction highly reduces the amount of errors present in the reads and allow genomic kmers to be more covered and therefore to "move" the curve to the right as shown in Figure 3.5.

Another way to remove sequencing errors it to remove "tips" or "short dead ends" from the graph. If $k \ge \text{ReadSize}/2$, an isolated sequencing error will not form a bubble but a tip in the graph. As we may be afraid of losing genomic kmers with the solidity cutoff, we advocate the use of low thresholds and therefore apply a tip removal step on the graph. We now evaluate in Tables 3.1 and 3.2 the number of false positives and false negatives in the produced graph based on different strategies using the following methods:

* Kmer filtering with a given threshold

* Using raw reads or reads corrected with the Bloocoo corrector [106]

* Tipping of the graph, removing tips shorter than the read length

We can see that the combination of kmer filtering, read correction and tipping can result in almost perfect reference graphs. Such tasks are the basic steps of assembly. Most
<figure><img src="image_52.png" /><figcaption>Figure 3.5: Raw reads versus corrected reads. Spectrum with various *k* from 21 to 81. Reads are simulated from *E. coli* reference. Reads were of length 100 and with a coverage of 100. The *k* values are 81 (yellow), 61 (blue), 41 (green), 21 (purple). The left plot is the kmer spectrum of reads without correction and the right plot is the kmer spectrum of the reads after correction with Bloocoo. We can see that kmer spectrum from corrected reads are more distant from the error peaks. Reads correction allows an easier erroneous kmers filtering. We also observe that higher kmer size may be difficult to filter as less coverage discrepancy is present to separate erroneous and genomic kmers.</figcaption></figure>

<table><caption>Table 3.1: Evaluation of the de Bruijn graph created from simulated sequencing from *E. coli*, 100x coverage of 100 base pairs reads. We constructed a de Bruijn graph with *k* = 51 and evaluated the number of erroneous kmers (FP) and missing kmers (FN). The abundance threshold (filtering) applied is indicated in the strategy description.</caption><thead><tr><th>Strategy</th><th># Unitig</th><th>FP</th><th>FN</th></tr></thead><tbody><tr><td>Filtering 2</td><td>434,510</td><td>3,876,482</td><td>8</td></tr><tr><td>Filtering 2 + tiping</td><td>2,474</td><td>29,515</td><td>8</td></tr><tr><td>Filtering 3</td><td>28,310</td><td>198,088</td><td>11</td></tr><tr><td>Filtering 3 + tiping</td><td>978</td><td>582</td><td>11</td></tr><tr><td>Filtering 5</td><td>1273</td><td>1716</td><td>19</td></tr><tr><td>Filtering 5 + tiping</td><td>944</td><td>17</td><td>19</td></tr><tr><td>Correction + filtering 2</td><td>6,860</td><td>16,202</td><td>2</td></tr><tr><td>Correction + filtering 2 + tiping</td><td>1,034</td><td>591</td><td>2</td></tr><tr><td>Correction + filtering 3</td><td>1,061</td><td>713</td><td>7</td></tr><tr><td>Correction + filtering 3 + tiping</td><td>956</td><td>73</td><td>7</td></tr><tr><td>Correction + filtering 5</td><td>952</td><td>59</td><td>9</td></tr><tr><td>Correction + filtering 5 + tiping</td><td>944</td><td>17</td><td>9</td></tr><tr><td>Reference</td><td>941</td><td>0</td><td>0</td></tr></tbody></table>