<table><caption>Table 3.2: Evaluation of the de Bruijn graph created from simulated sequencing from *C. elegans*, 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># Unitigs</th><th>FP</th><th>FN</th></tr></thead><tbody><tr><td>Filtering 2</td><td>10,167,089</td><td>89,087,143</td><td>6</td></tr><tr><td>Filtering 2 + tiping</td><td>3,580,486</td><td>2,696,841</td><td>55</td></tr><tr><td>Filtering 3</td><td>1,103,166</td><td>7,294,417</td><td>7</td></tr><tr><td>Filtering 3 + tiping</td><td>432,648</td><td>813,483</td><td>30</td></tr><tr><td>Filtering 5</td><td>303,908</td><td>1,049,812</td><td>15</td></tr><tr><td>Filtering 5 + tiping</td><td>204,481</td><td>231,083</td><td>33</td></tr><tr><td>Reference</td><td>64,044</td><td>0</td><td>0</td></tr></tbody></table>

tools add other simplifications and heuristics to produce larger contigs. While there is a plethora of assemblers, we would find relevant to propose tools in order to create a clean and safe graph. This kind of unitigs generator, "unitigers", could be used as a front end by assemblers that could apply their different post-treatments according to their specificity directly on the cleaned graph. In the following, we will rely on such steps to get reference de Bruijn graphs with almost no sequencing error and no missed kmer.

**Efficient representation of redundancy** We have previously seen that tools exist to construct the de Bruijn graph in an efficient way. The graph constructions of table 3.1 on *E. coli* took less than 10 minutes with BCALM2. The usage of a (non compacted) de Bruijn graph does not seem interesting since a genome will require ≈ *genomeSize* kmer to be stored. On the other hand, the compacted de Bruijn graph can efficiently represent a genome since each nucleotide will require 2 bits plus the global overhead of the unitigs number.

As shown in the table 3.3 that compacted de Bruijn graph is a space efficient genome representation. Based on such observations, some tools were developed in order to efficiently construct the de Bruijn graph as the representation of a genome or multiple genomes [107] [108] [109]. We already presented the efficiency of BCALM2 in order to construct a compacted de Bruijn graph but it was designed to handle very large reads sets and is not necessarily efficient for de Bruijn graph construction from reference genomes. Twopaco [109] proposed an efficient approach to build such graphs from many reference genomes and was able to construct the de Bruijn graph of one hundred human genomes in less than a day. This method to construct a de Bruijn graph from a reference genome without splitting the reference into its set of kmers is way more efficient than BCALM2.

To access the efficiency of the de Bruijn graph to represent multiple genome, we downloaded 100 *E. coli* completed genome from NCBI and created the de Bruijn graph of those merged reference genomes (Table 3.4). We can see that the de Bruijn graph is a compact data structure to store a large number of genomes even if a larger *k* present a less interesting compression ratio. For comparison the reference file zipped represented 153 mega bytes. Indeed the de Bruijn graph representation is lossy and can not be compared to the lossless zip compression but [110] argue that this efficient redundancy factoring
<table><caption>Table 3.3: Nucleotides number in de Bruijn graph representations of a genome.</caption><thead><tr><th>Genome</th><th>Structure</th><th>Number of nucleotides used for representation</th></tr></thead><tbody><tr><td>E. coli</td><td>Reference sequence</td><td>4,639,675</td></tr><tr><td>E. coli</td><td>de Bruijn graph k=51</td><td>232,770,375</td></tr><tr><td>E. coli</td><td>de Bruijn graph k=101</td><td>462,106,108</td></tr><tr><td>E. coli</td><td>Compacted de Bruijn graph k=51</td><td>4,611,175</td></tr><tr><td>E. coli</td><td>Compacted de Bruijn graph k=101</td><td>4,619,908</td></tr><tr><td>C. elegans</td><td>Reference sequence</td><td>100,286,401</td></tr><tr><td>C. elegans</td><td>de Bruijn graph k=51</td><td>4,889,975,931</td></tr><tr><td>C. elegans</td><td>de Bruijn graph k=101</td><td>9,864,383,665</td></tr><tr><td>C. elegans</td><td>Compacted de Bruijn graph k=51</td><td>104,537,731</td></tr><tr><td>C. elegans</td><td>Compacted de Bruijn graph k=101</td><td>103,927,665</td></tr><tr><td>Human</td><td>de Bruijn graph k=51</td><td>138,318,741,180</td></tr><tr><td>Human</td><td>de Bruijn graph k=101</td><td>286,185,334,362</td></tr><tr><td>Human</td><td>Compacted de Bruijn graph k=51</td><td>2,809,115,473</td></tr><tr><td>Human</td><td>Compacted de Bruijn graph k=101</td><td>3,083,875,662</td></tr></tbody></table>

<table><caption>Table 3.4: Size of the de Bruijn graph of 100 E. coli genomes according to k value used for its construction.</caption><thead><tr><th>Dataset</th><th># Nucleotides</th></tr></thead><tbody><tr><td>Reference file</td><td>485,354,804</td></tr><tr><td>k=25</td><td>55,864,796</td></tr><tr><td>k=51</td><td>89,669,228</td></tr><tr><td>k=101</td><td>137,742,996</td></tr><tr><td>k=201</td><td>188,631,137</td></tr><tr><td>k=301</td><td>215,417,352</td></tr></tbody></table>
<figure><img src="image_53.png" /><figcaption>Figure 3.6: Example of nested repeats. In this example, a large repeat in blue contains a small repeat in red. This scenario creates complex cases where the repeats patterns can be nested in one another's.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Figure 3.7: Example of tandem repeats. In this example, multiple occurrences of a repeat are consecutive in the genome. This scenario creates cases where it is hard to estimate the number of times the repeat sequence should be present in the assembly.</figcaption></figure>

propriety could be useful for pan-genome analysis.

**Complete reads information** Assembly, based on the de Bruijn graph or not, may be difficult, especially around repeats or more complex patterns as quasi-repeats (Figure 1.20), nested (Figure 3.6) or tandem repeats (Figure 3.7). Even read mapping is hard around such complex repeats and some regions can show mappability problem [111], since reads can have several possible mappings. The de Bruijn graph representation allow a reduced issue as the repeats of the genome are collapsed. In the cases of an unfinished genome represented by a set of contigs, the problem is different. Contigs are produced from reads, using various heuristics to select path from an assembly graph. This selection may induce biases according to the heuristics used, especially around complex patterns. For example a read not mapped on the contigs set may be due to missing sequences in the contigs outputed. The assemblers biases could be problematic and we argue that mapping directly on a de Bruijn graph can provide a less biased and more complete reference.

In those very complex zones, most assemblers are incapable to output large contigs. Most of the time those regions are absent from the contig set, because they are not outputted by the assembler. Indeed assembler that output all contigs produce small contigs in such cases. One conclusion is that scaffolding may be impracticable on such situations because it rely on read mapping on large contigs. The situation of those region can not be improved by usual scaffolding methods, and are lost in most assembly processes. In such regions, the graph is complex and branching but the genomic sequences are still paths of the graph. The information carried by the graph is then superior to the
one brought by contigs, as the graph contains information either on branching regions or on how the contigs could be ordered.

A second point is that the graph can carry the haplotypes information. In presence of multiple closely related genomes or haplotypes, minor variations will be represented in bubbles in the graph in an efficient way. De Bruijn graph is an efficient structure to capture variant information and several methods use the de Bruijn graph to discover polymorphisms among several individuals [90] [112] [99] [113]. This property has already been used to index several genomes in order to avoid multiple alignment in whole genomes comparisons [114] [107] and metagenomic quantification [103].

### "Genome representations" core messages:

* Linear sequences may not be the ideal representation of a genome

* Several graphs based structures have been proposed to fit various applications to represent finished or almost finished genomes

* The de Bruijn graph may be an efficient structure for applications where a reference genome is not available

## 3.2 Read mapping on the de Bruijn graph

### 3.2.1 An efficient tool for an NP-Complete problem

As we have seen, the de Bruijn graph can be efficiently used as a reference. Read mapping is a core operation we want to perform on a reference and most tools are designed to work on a reference genome represented by flat sequences, fragmented or not [115] [60] [59]. Those tools are able to index and to map reads on a set of sequences. The problematic of read mapping is to know if a read can be aligned on the reference, where and with an indication on the alignment quality. As we mentioned before, some complex repeated regions sometimes occur in genomes, and suffer from low mappability [111]. Reads from repeated region may have multiple matching sites and then be hard to be mapped with high confidence. When a genome is represented as a de Bruijn graph, repeated regions are merged in the graph. This factorization reduces the problem of multiple mapping. This de Bruijn graph property explains its usage in genomic or meta-genomic, even on finished reference genomes [105]. We argue that being able to map reads directly on such structures would be less biased and more complete than mapping on contigs.

Another interest of mapping reads on a de Bruijn graph is to improve the de Bruijn graph itself to produce a better assembly. This read information added in the graph can be used to avoid false connections between nodes or to solve some repeats. Surprisingly, despite the interest in such techniques, no practical solution has been designed for this task. Assemblers using this kind of information do not present generic procedure or rely on alignment on flat sequences. All those different usages motivates the need of such tools.

This paper proposes a formal definition of the problem of mapping reads on a de Bruijn graph and proves its NP-completeness. We also provides a practical solution based on several heuristics, called BGREAT.
Technical results are described in the paper "Read mapping on de Bruijn graphs" [116] included at the end of the chapter.

BGREAT is based on simple heuristics in order to be fast and scalable. It follows the seed and extend paradigm by indexing the *k* − 1 suffix and prefix of the unitigs as anchors. Once an anchor is found on a read, the rest of the read sequence is mapped to the graph unitigs. By indexing the prefix and suffix and associating them the indices of the unitigs sharing them, BGREAT is able to navigate through the graph unitigs as it was presented in the last chapter. Some unitigs may not contain any anchor, for example if they map entirely on large unitigs. But in such cases those reads can be mapped with a regular mapping tool. Another essential heuristic of BGREAT is the greediness of its mapping process. When several paths can be used, the path with the minimal number of mismatches is picked. The use of those greedy heuristics allows very fast mapping of reads even on a large de Bruijn graph. BGREAT is able to map a human dataset of 3 billion reads on a de Bruijn graph created from it, in less than 5 hours using 20 cores and 10GB. Another interesting conclusion is that we are able to map more reads on the de Bruijn graph than on a set of contigs. An experiment on a human dataset showed that only 63% of the reads were mapped on the contig set where 85% were mapped on the de Bruijn graph using the graph alignment method.

### 3.2.2 Mapping refinements

The amount of potential applications of read alignment on de Bruijn graph leads us to the conception of several improvements of the BGREAT tool. We propose an advanced version called BGREAT2 (unpublished) with following improvements that we will detail:

* Improved anchor system

* Indexing based on a MPHF

* Optimal mapping and multi-mapping management

In order to improve the performance of our anchoring scheme, we no longer use the *k*-1mer at the extremities of the unitigs. We index kmers from the unitigs and associate to each kmer its position in the unitigs. This modification allows two significant changes.

First, the size of the anchor can be chosen. Therefore a smaller *k* than the one selected to construct the graph can be used. The order of the de Bruijn graph may be quite large, thus a *k*-1mer would be a bad anchor in terms of sensibility. This way, a high order graph can be used without presenting a change in the anchoring scheme. Indexing all kmers of the graph can be costly, for a human graph indexing all kmers can use up to 200 GB. A parameter can be defined in order to reduce the memory usage of BGREAT2 that regulates which fraction of the kmers are indexed. We can show that even by indexing only a fraction of the kmers (one out of 10 by default), the mapping performances of BGREAT2 are not highly impacted. Indexing 1 out of 10 anchors in each unitig, with at least one anchor indexed by unitigs allowed BGREAT2 to reduce its memory usage from 6,898 MB to 928 MB to index a *C. elegans* reference graph. The indexing time also dropped from 61 to 10 seconds while the ratio of mapped reads went from 99.6% to 99.3%. The mapping time was also impacted, but in a less impressive way and improved from 11k to 15k reads by second on a single thread.
The second interest of the new indexing scheme is that BGREAT2 is self contained
and do not need a mapping tool like Bowtie. Since kmers are indexed along the whole
length of the unitigs, reads mapping inside a unitig can be anchored and mapped. In
order to index kmers and the de Bruijn graph overlaps, BGREAT2 no longer uses a
dynamic hash table but the BBhash MPHF. This choice allows a highly reduced memory
usage and a faster access to the unitigs compared to previously used dynamic hashing.

We also improved the mapping rules. The algorithm tries to find a perfect mapping
(without error), then a mapping with one error and so on until the maximal error number
allowed. This way, the obtained mapping is guaranteed to contain a minimal number of
mismatch. Several multiple mapping strategies are also proposed: output all equivalent
mapping, output one mapping, or none. Albeit some optimizations and features have still
to be added, those improvements allow BGREAT2 to be more than a proof of concept.
We will present in the following sections several uses of this tool.

"Read mapping on the de Bruijn graph" core messages:

* Mapping reads on a de Bruijn graph is a hard problem

* With appropriate heuristics we can propose a satisfying and efficient mapping

* Mapping on the de Bruijn graph may allow more mapping hits than mapping on contigs

## 3.3 De novo, reference guided, read correction

We have shown previously that an almost errorless de Bruijn graph could be constructed.
We therefore propose a new method of correction for short reads by aligning them on
a reference represented by a de Bruijn graph constructed from the reads themselves.
Correction methods based on alignment on a reference genome would have the two same
essential drawbacks than reference guided assembly:

* Need for a reference

* Correction may be biased by the reference used

Here we propose to construct the reference directly from the reads to be corrected in
order to avoid biases. We argue that we do not need to produce a good assembly in order
to obtain a reference good enough to correct the reads. The idea of such a technique can
be summarized in the following steps:

* Graph construction

* Graph cleaning

* Read mapping on the graph

### 3.3.1 Limits of kmer spectrum

Most state of the art correctors are based on kmer spectrum techniques. Other techniques
based on suffix array or on multiple alignment, are presented in [117]. Those techniques
are less used because they rely on memory expensive data structures and do not scale
well on large genomes and datasets. Kmer spectrum correction was proposed when NGS
<figure><img src="image_55.png" /><figcaption>Figure 3.8: Example of read correction by mapping it on a de Bruijn graph. The query read is mapped onto the graph and the bases of the graph are trusted over the read ones. The read is thus corrected according to the reference graph. An alternative mapping for the reads would be ATTCACCCGG but this path includes two mismatches and is therefore not optimal. This example also shows a problem of kmer spectrum technique. In the read to correct, the sequencing error (the red T) is validated by the genomic kmers TTCT, TCTC and CTCC (framed in red). According to heuristics employed by those tools, this base is unlikely to be a sequencing error and will not be corrected.</figcaption></figure>

reads appeared and were assembled with de Bruijn graph [46]. As in de Bruijn graph assembly, a set of "solid" kmers are computed. Those highly abundant kmers are very likely to be errorless and are indexed. When a "non solid" kmer is found in a read, the tool tries to convert it into a solid kmer in the most parsimonious way. Since such correctors only index trusted kmers against low abundance ones, bases of a read covered by trusted kmers will be considered as correct. In the example of Figure 3.8, all bases of the read will be considered correct because they are covered by several solid kmers. When an erroneous base is found, the algorithms try to replace the erroneous kmers by solid kmers. A hard part for those kind of algorithms is when close errors are found. If multiple errors appear on a single kmer, not all the possible alternative kmers can be tested and no real satisfying solution has been found for such cases.

### 3.3.2 Reads correction by de Bruijn graph mapping

Here we present BCOOL, a still unpublished proof of concept short reads corrector based on read mapping on a de Bruijn graph. The idea is to create an almost errorless de Bruijn graph from the reads to be corrected and use it as a reference to correct the reads mapped on it. After the graph cleaning step, we map each read on the graph and when a mismatch is detected, the nucleotide of the graph is trusted over the nucleotide of the read, and the read is modified to match with the de Bruijn graph sequence.

What are the differences between the kmer spectrum correction technique and this proposition? Both kmer spectrum and our approach use a set of trusted kmers as reference to correct the reads. The difference is that we use the structure of the de Bruijn graph, i.e. how the kmers are connected. There are several advantages to use this information. Tip removal can typically remove sequencing errors above the solidity threshold. It allows our set of solid kmers to contain less erroneous kmers. We also have a better comprehension of complex cases. In Figure 3.8 a kmer spectrum method would not correct anything since the sequencing error is covered by multiple solid kmers and would consider it safe. Our
<table><caption>Table 3.5: Correction benchmark on simulated 100bp reads (100X coverage 1% error rate) from *E. coli* by varying the *k* parameter. Solidity threshold to construct the de Bruijn graph was fixed to 5 and tipping length to 100.</caption><thead><tr><th>Kmer size used</th><th>FP</th><th>FN</th><th>TP</th><th>Correction ratio</th><th>% Erroneous read</th></tr></thead><tbody><tr><td>BCOOL k=21</td><td>2,731</td><td>11,235</td><td>4,627,370</td><td>332</td><td>0.16</td></tr><tr><td>BCOOL k=31</td><td>1,092</td><td>2,427</td><td>4,636,805</td><td>1,318</td><td>0.04</td></tr><tr><td>BCOOL k=41</td><td>915</td><td><b>1,593</b></td><td><b>4,637,775</b></td><td><b>1,850</b></td><td><b>0.03</b></td></tr><tr><td>BCOOL k=51</td><td>1,150</td><td>2,183</td><td>4,637,245</td><td>1,392</td><td>0.03</td></tr><tr><td>BCOOL k=61</td><td>2,146</td><td>2,590</td><td>4,636,858</td><td>980</td><td>0.04</td></tr><tr><td>BCOOL k=71</td><td>2,452</td><td>43,202</td><td>4,596,267</td><td>102</td><td>0.58</td></tr><tr><td>BFC</td><td>2,105</td><td>78,679</td><td>4,559,123</td><td>57</td><td>0.67</td></tr><tr><td>Bloocoo</td><td>4,147</td><td>209,071</td><td>4,430,443</td><td>22</td><td>3.49</td></tr><tr><td>Lighter</td><td><b>624</b></td><td>6,538</td><td>4,633,059</td><td>648</td><td>0.09</td></tr><tr><td>Musket</td><td>1,031</td><td>6,456</td><td>4,632,901</td><td>620</td><td>0.13</td></tr></tbody></table>

method maps the whole read, making use of the sequences next to the error to correct it. Close errors can be problematic for other correctors, when multiple errors occur on a kmer the number of kmers able to replace it is exponential. It would be difficult to check the solidity of all possible kmers to replace the erroneous one. Our alignment method is not impacted by such cases since we do not correct kmers but entire reads. In order to assess the potential of such a technique we tested it against state of the art read correctors. In our benchmarks we will compare those four correctors Musket [118], BFC [119], Lighter [120] and Bloocoo [106] for their performances in our tests and in several published benchmarks.

We generated simulated reads with known errors, corrected them and evaluated the correction results. We use several metrics to evaluate the correction quality. The false positives are correct bases that have been modified or incorrect bases that have been wrongly modified. It can be seen as errors introduced by the correction. The false negative are sequencing errors that were not modified. It can be seen as not corrected errors. The true positives are sequencing errors correctly corrected. To provide a global measure of the correction, we evaluate the number of errors before and after the correction, counting FP and FN. A ratio of 10 mean that the amount of error was divided by ten by the correction. We also present the read percentage that contain at least an error.

The first results on the *E. coli* bacteria are reported in Table 3.5. We can observe that most tools propose a very good correction. All tools achieve to divide the error rate by more than 10, showing that they take care of most sequencing errors. Results also show that almost all reads are errorless after the correction. We see that BCOOL is able to provide excellent correction with adapted parameters. Globally, we argue that correctors provide very good corrections on small genomes.

The results on the larger genome of *C. elegans* are presented in Tables 3.6 and 3.7 on reads of length 100 and 250 respectively. Here we see that correctors still provide good correction by dividing the global error rate and proposing a vast majority of perfect reads but in a less impressive way than on *E. coli*. Here again BCOOL provides very interesting results on both read length of 100 and 250.

The results on the human genome are presented in Table 3.8. We see that BCOOL is
<table><caption>Table 3.6: Correction benchmark on simulated 100bp reads (100X coverage 1% error rate) from *C. elegans* by varying the *k* parameter. Solidity was fixed to 5 and tipping length to 100.</caption><thead><tr><th>Kmer size used</th><th>FP</th><th>FN</th><th>TP</th><th>Ratio error rate</th><th>% erroneous read</th></tr></thead><tbody><tr><td>BCOOL k=21</td><td>1,455,866</td><td>6,185,857</td><td>93,792,770</td><td>13</td><td>4.05</td></tr><tr><td>BCOOL k=31</td><td>1,420,047</td><td>1,529,334</td><td>98,558,023</td><td>34</td><td>1.57</td></tr><tr><td>BCOOL k=41</td><td>883,155</td><td>753,403</td><td>99,419,269</td><td>61</td><td>0.87</td></tr><tr><td>BCOOL k=51</td><td>651,040</td><td>498,536</td><td>99,716,829</td><td>87</td><td>0.59</td></tr><tr><td>BCOOL k=61</td><td>625,830</td><td><strong>447,835</strong></td><td><strong>99,795,401</strong></td><td><strong>93</strong></td><td><strong>0.49</strong></td></tr><tr><td>BCOOL k=71</td><td>803,946</td><td>1,214,915</td><td>99,042,681</td><td>50</td><td>0.97</td></tr><tr><td>BFC</td><td><strong>289,081</strong></td><td>4,020,868</td><td>96,187,218</td><td>23</td><td>2.23</td></tr><tr><td>Bloocoo</td><td>800,396</td><td>6,471,767</td><td>93,805,016</td><td>14</td><td>5.23</td></tr><tr><td>Lighter</td><td>613,314</td><td>2,789,652</td><td>97,490,978</td><td>29</td><td>2.16</td></tr><tr><td>Musket</td><td>1,120,401</td><td>4,932,597</td><td>95,250,573</td><td>17</td><td>4.47</td></tr></tbody></table>

<table><caption>Table 3.7: Correction benchmark on simulated 250bp reads (100X coverage 1% error rate) from *C. elegans* by varying the *k* parameter. Solidity was fixed to 5 and tipping length to 100.</caption><thead><tr><th>Kmer size used</th><th>FP</th><th>FN</th><th>TP</th><th>Ratio error rate</th><th>% erroneous read</th></tr></thead><tbody><tr><td>BCOOL k=21</td><td>1,048,454</td><td>7,723,461</td><td>92,278,080</td><td>11</td><td>7.72</td></tr><tr><td>BCOOL k=31</td><td>1,039,760</td><td>3,501,749</td><td>96,582,586</td><td>22</td><td>4.20</td></tr><tr><td>BCOOL k=41</td><td>787,674</td><td>2,099,774</td><td>98,052,234</td><td>35</td><td>2.67</td></tr><tr><td>BCOOL k=51</td><td>608,383</td><td>1,424,096</td><td>98,771,993</td><td>49</td><td>1.87</td></tr><tr><td>BCOOL k=61</td><td>500,090</td><td>1,132,677</td><td>99,092,022</td><td>61</td><td>1.49</td></tr><tr><td>BCOOL k=71</td><td>411,359</td><td>996,580</td><td>99,247,373</td><td>71</td><td>1.28</td></tr><tr><td>BCOOL k=81</td><td>343,784</td><td>950,693</td><td>99,304,838</td><td>77</td><td>1.27</td></tr><tr><td>BCOOL k=91</td><td>301,430</td><td><strong>944,491</strong></td><td><strong>99,318,476</strong></td><td><strong>80</strong></td><td><strong>1.12</strong></td></tr><tr><td>BCOOL k=101</td><td>271,094</td><td>978,287</td><td>99,291,125</td><td>80</td><td>1.12</td></tr><tr><td>Musket</td><td>812,832</td><td>3,691,147</td><td>96,498,970</td><td>22</td><td>7.22</td></tr><tr><td>Lighter</td><td>520,365</td><td>3,022,271</td><td>97,254,509</td><td>28</td><td>4.76</td></tr><tr><td>Bloocoo</td><td>873,717</td><td>5,746,372</td><td>94,519,928</td><td>15</td><td>10.38</td></tr><tr><td>BFC</td><td><strong>176,279</strong></td><td>6,289,736</td><td>93,933,429</td><td>16</td><td>5.26</td></tr></tbody></table>

<table><caption>Table 3.8: Correction benchmark on simulated 100bp reads (100X coverage 1% error rate) from Human by varying the *k* parameter. Solidity was fixed to 5 and tipping length to 100.</caption><thead><tr><th>Kmer size used</th><th>FP</th><th>FN</th><th>TP</th><th>Ratio error rate</th><th>% erroneous read</th></tr></thead><tbody><tr><td>BCOOL k=41</td><td>186,222,378</td><td>165,218,857</td><td>2,346,232,116</td><td>7</td><td>5.16</td></tr><tr><td>BCOOL k=51</td><td>137,556,747</td><td>115,336,678</td><td>2,400,321,919</td><td>10</td><td>3.52</td></tr><tr><td>BCOOL k=61</td><td>99,095,091</td><td>88,713,912</td><td>2,429,964,113</td><td>13</td><td>2.47</td></tr><tr><td>BCOOL k=71</td><td><strong>85,492,564</strong></td><td><strong>79,637,965</strong></td><td><strong>2,440,258,246</strong></td><td><strong>15</strong></td><td><strong>2.10</strong></td></tr><tr><td>Lighter</td><td>55,280,808</td><td>883,251,227</td><td>1,632,909,855</td><td>3</td><td>27.83</td></tr><tr><td>Bloocoo</td><td>80,025,724</td><td>374,402,778</td><td>2,148,076,789</td><td>6</td><td>11.16</td></tr><tr><td>BFC</td><td>88,500,362</td><td>962,587,329</td><td>1,497,440,057</td><td>2</td><td>27.59</td></tr></tbody></table>
able to scale up to large genome and dataset as the human one. Musket needs more than 250 GB and was not able to correct this dataset on our testing server. Once again we see that correctors provide less impressive results on larger genomes. Those results show that BCOOL provides good read correction even compared to state of the art correctors.

Our proposed workflow can be summarized into:

* Best kmer size and abundance threshold selection (Not done yet)

* Graph construction

* Graph cleaning

* Reads mapping on the graph

We see that BCOOL can lead to significant improvements to reads correction and can scale up to very large genomes. But the approach is still parameter dependent. A high k value allows the construction of a high quality graph but can create holes. Furthermore we only evaluated BCOOL on very high coverage dataset and the solidity threshold is also an impactful parameter that can be non trivial to set for users not aware of de Bruijn graph properties. An automated evaluation of the kmer spectrum could infer the good values to set, in order to get rid of most sequencing error without losing genomic kmers. In order to infer those parameters we could, as in Kmergenie, use a fast kmer spectrum estimator such as ntCard [121]. Many improvements and tuning have yet to be performed to propose an efficient and robust corrector but preliminary results are encouraging.

"De novo, reference guided, read correction" core messages:

* Mapping reads on a reference de Bruijn graph enables to correct them

* All operations of such a method are scalable
<figure><img src="image_56.png" /></figure>

# Read mapping on de Bruijn graphs

Antoine Limasset<fnref n="1" />*, Bastien Cazaux²,³, Eric Rivals²,³ and Pierre Peterlongo<fnref n="1" />

**Abstract**

**Background:** Next Generation Sequencing (NGS) has dramatically enhanced our ability to sequence genomes, but not to assemble them. In practice, many published genome sequences remain in the state of a large set of contigs. Each contig describes the sequence found along some path of the assembly graph, however, the set of contigs does not record all the sequence information contained in that graph. Although many subsequent analyses can be performed with the set of contigs, one may ask whether mapping reads on the contigs is as informative as mapping them on the paths of the assembly graph. Currently, one lacks practical tools to perform mapping on such graphs.

**Results:** Here, we propose a formal definition of mapping on a de Bruijn graph, analyse the problem complexity which turns out to be NP-complete, and provide a practical solution. We propose a pipeline called **GGMAP** (Greedy Graph MAPPing). Its novelty is a procedure to map reads on branching paths of the graph, for which we designed a heuristic algorithm called **BGREAT** (de Bruijn Graph REAd mapping Tool). For the sake of efficiency, **BGREAT** rewrites a read sequence as a succession of unitigs sequences. **GGMAP** can map millions of reads per CPU hour on a de Bruijn graph built from a large set of human genomic reads. Surprisingly, results show that up to 22 % more reads can be mapped on the graph but not on the contig set.

**Conclusions:** Although mapping reads on a de Bruijn graph is complex task, our proposal offers a practical solution combining efficiency with an improved mapping capacity compared to assembly-based mapping even for complex eukaryotic data.

**Keywords:** Read mapping, De Bruijn graph, NGS, Sequence graph, path, Hamiltonian path, Genomics, Assembly, NP-complete

**Background**

Next Generation Sequencing technologies (NGS) have drastically accelerated the generation of sequenced genomes. However, these technologies remain unable to provide a single sequence per chromosome. Instead, they produce a large and redundant set of reads, with each read being a piece of the whole genome. Because of this redundancy, it is possible to detect overlaps between reads and to assemble them together in order to reconstruct the target genome sequence.

Even today, assembling reads remains a complex task for which no single piece of software performs consistently well [1]. The assembly problem itself has been shown to be computationally difficult, more precisely NP-hard [2]. Practical limitations arise both from the structure of

genomes (repeats longer than reads cannot be correctly resolved) and from the sequencing biases (non-uniform coverage and sequencing errors). Applied solutions represent the sequence of the reads in an assembly graph: the labels along a path of the graph encode a sequence. Currently, most assemblers rely on two types of graphs: either the de Bruijn graph (DBG) for the short reads produced by the second generation of sequencing technologies [3], or for long reads the overlap graph (which was introduced in the Celera Assembler [4]) and variants thereof, like the string graph [5]. Then, the assembly algorithm explores the graph using heuristics, selects some paths and outputs their sequences. Due to these heuristics, the set of sequences obtained, called contigs, is biased and fragmented because of complex patterns in the graph that are generated by sequencing errors, and genomic variants and repeats. The set of contigs is rarely satisfactory and is

<footnote>*Correspondence: antoine.limasset@irisa.fr</footnote>

<footnote n="1">¹IRISA Inria Rennes Bretagne Atlantique, GenScale team, Campus de Beaulieu, 35042 Rennes, France</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
usually post-processed, for instance, by discarding short contigs.

The most frequent computational task for analyzing a set of reads is mapping them on a reference genome. Numerous tools are available to map reads when the reference genome has the form of a set of sequences (e.g. BWA [6] and Bowtie [7]). The goal of mapping on a finished genome sequence is to say whether a sequence can be aligned to this genome, and in this case, at which location(s). This is mostly done with a heuristic (semi-global) alignment procedure that authorizes a small edit or Hamming distance between the read and genome sequences. Read mapping process suffers from regions of low mappability [8]. Repeated genomic regions may not be mapped precisely since the reads mapping on these regions have multiple matches. When a genome is represented as a graph, the mappability issue is reduced, as occurrences of each repeated region are factorized, limiting the problem of multiple matches of reads.

When the reference is not a finished genome sequence, but a redundant set of contigs, the situation differs. The mapping may correctly determine whether the read is found in the genome, but multiple locations may for instance not be sufficient to conclude whether several true locations exist. Conversely, an unfruitful mapping of a read may be due to an incomplete assembly or to the removal of some contigs during post-processing. In such cases, we argue it may be interesting to consider the assembly graph as a (less biased and/or more complete) reference instead of the set of contigs. Then mapping on the paths of this graph is needed to complement mapping on set of contigs. This motivates the design and implementation of BGREAT.

In this context, we explore the problem of mapping reads on a graph. Aligning or mapping sequences on sequence graphs (a generic term meaning a graph representing sequences along its paths) has already been explored in the literature in different application contexts: assembly, read correction, or metagenomics.

In the context of assembly, once a **DBG** has been built, mapping the reads back to the graph can help in eliminating unsupported paths or in computing the coverage of edges. To our knowledge, no practical solution has been designed for this task. Cerulean assembler [9] mentions this possibility, but only uses regular alignment on assembled sequences. Allpaths-LG [10] also performs a similar task to resolve repeats using long noisy reads from third generation sequencing techniques. Its procedure is not generic enough to suit the mapping of any read set on a **DBG**. From the theoretical view point, the question is related to the NP-hard *read-threading* problem (also termed *Eulerian superpath problem* [2, 11]), which consists in finding a read coherent path in the **DBG** (a path that can be represented as a sequence of reads as

defined in [5]). The assembler called SPADES [12] threads the reads against the **DBG** by keeping track of the paths used during construction, which requires a substantial amount of memory. Here, we propose a more general problem, termed *De Bruijn Graph Read Mapping Problem* (**DBGRMP**), as we aim at mapping to a graph any source of NGS reads, either those reads used for building the graph or other reads.

Recently, the hybrid error correction of long reads using short reads has become a critical step to leverage the third generation of sequencing technologies. The error corrector LoRDEC [13] builds the **DBG** of the short reads, and then aligns each long read against the paths of the **DBG** by computing their edit distance using a dynamic programming algorithm (which is slow for our purposes). For short reads correction, several tools that evaluate the *k*-mer spectrum of reads to correct the sequencing errors use a probabilistic or an exact representation of a **DBG** as a reference [14, 15].

In the context of metagenomics, Wang et al. [16] have estimated the taxonomic composition of a metagenomics sample by mapping reads on a **DBG** representing several genomes of closely-related bacterial species. In fact, the graph collapses similar regions of these genomes and avoids redundant mapping. Their tool maps the read using BWA on the sequence resulting from the random concatenation of unitigs of the **DBG**. Hence, a read cannot align over several successive nodes of the graph (ER: *il y a un pb ce n'est pas vrai*). Similarly, several authors have proposed to store related genomes into a single, less repetitive, **DBG** [17–19]. However, most of these tools are efficient only when applied to very closely related sequences that result in flat graphs. The *BlastGraph* tool [19], is specifically dedicated to the mapping of reads on graphs, but is unusable on real world graphs (see Results section).

Here, we formalize the mapping of reads on a De Bruijn graph and show that it is NP-complete. Then we present the pipeline **GGMAP** and dwell on **BGREAT**, a new tool which enables to map reads on branching paths of the **DBG** (Section **GGMAP**: a method to map reads on de Bruijn Graph). For the sake of efficiency, **BGREAT** adopts a heuristic algorithm that scales up to huge sequencing data sets. In Section Results, we evaluate **GGMAP** in terms of mapping capacity and of efficiency, and compare it to mapping on assembled contigs. Finally, we discuss the limitations and advantages the of **GGMAP** and give some directions of future work (Section Discussion).

## Methods

We formally define the problem of mapping reads on a **DBG** and investigate its complexity (Section Complexity of mapping reads on the paths of a **DBG**). Besides, we propose a pipeline called **GGMAP** to map short
reads on a representation of a `DBG` (Section GGMAP: a method to map reads on de Bruijn Graph). This pipeline includes `BGREAT`, a new algorithm mapping sequences on branching paths of the graph (Section `BGREAT`: mapping reads on branching paths of the `CDBG`).

## Complexity of mapping reads on the paths of a `DBG`

In this section, we present the formal problem we aim to solve and prove its intractability. First, we introduce preliminary definitions, then formalize the problem of mapping reads on paths of a `DBG`, called the De Bruijn Graph Read Mapping Problem (`DBGRMP`), and finally prove it is NP-complete. Our starting point is the well-known Hamiltonian Path Problem (HPP); we apply several reductions to prove the hardness of `DBGRMP`.

**Definition 1** (de Bruijn graph). Given a set of strings $S = \{r_1, r_2, \dots, r_n\}$ on an alphabet $\Sigma$ and an integer $k \ge 2$, the de Bruijn graph of order $k$ of $S$ (`dBG_k(S)`) is a directed graph $(V, A)$ where:

$V = \{d \in \Sigma^k | \exists i \in \{1, \dots, n\} \text{ such that } d \text{ is a substring of } r_i \in S\}, \text{ and}$

$A = \{(d, d') | \text{ if the suffix of length } k-1 \text{ of } d \text{ is a prefix of } d'\}$.

**Definition 2** (Walk and Path of a directed graph). Let $G$ be a directed graph.

* A walk of $G$ is an alternating sequence of nodes and connecting edges of $G$.

* A path of $G$ is a walk of $G$ without repeated node.

* A Hamiltonian path is a path that visits each node of $G$ exactly once.

**Definition 3** (Sequence generated by a walk in a `dBG_k`). Let $G$ be a de Bruijn graph of order $k$. A walk of $G$ composed of $l$ nodes ($v_1, \dots, v_l$) generates a sequence of length $k+l-1$ obtained by the concatenation of $v_1$ with the last character of $v_2$, of $v_3$, \dots, of $v_l$.

We define the *de Bruijn Graph Read Mapping Problem* (`DBGRMP`) as follows:

**Definition 4** (De Bruijn Graph Read Mapping Problem). Given

* $S$, a set of strings over $\Sigma$,

* $k$, an integer such that $k \ge 2$,

* $q := q_1 \dots q_{|q|}$ a word of $\Sigma^*$ such that $|q| \ge k$,

* a cost function $F: \Sigma \times \Sigma \to \mathbb{N}$, and

* a threshold $t \in \mathbb{N}$,

decide whether there exists a path of the dBG_k(S) composed of $|q| - k + 1$ nodes (generating a word $m :=$

$m_1 \dots m_{|q|} \in \Sigma^{|q|})$ such that the cost $C(m, q) := \sum_{i=1}^{|q|} F(m_i, q_i) \le t$.

We recall the definition of the *Hamiltonian Path Problem* (HPP), which is NP-complete [20].

**Definition 5** (*Hamiltonian Path Problem* (HPP)). Given a directed graph $G$, the HPP consists in deciding whether there exists a Hamiltonian path of $G$.

To prove the NP-completeness of `DBGRMP` we introduce two intermediate problems. The first problem is a variant of the Asymmetrical Travelling Salesman Problem.

**Definition 6** (*Fixed Length Asymmetric Travelling Salesman Problem* (FLATSP)). Let

* $l$ be an integer,

* $G := (V, A, c)$ be a directed graph whose edges are labeled with a non-negative integer cost (given by the function $c: A \to \mathbb{N}$),

* $t \in \mathbb{N}$ be a threshold.

`FLATSP` consists in deciding whether there exists a path $p := (v_1, \dots, v_l)$ of $G$ composed of $l$ nodes whose cost $c(p) := \sum_{j=1}^{l-1} c((v_j, v_{j+1}))$ satisfies $c(p) \le t$.

We consider the restriction of `FLATSP` to instances having a unit cost function (i.e., where $c(a) = 1$ for any $a \in A$) and where $l$ equals both the threshold and the number of nodes in $V$. This restriction makes `FLATSP` very similar to HPP, and the hardness result quite natural.

**Proposition 1.** `FLATSP` is NP-complete even when restricted to instances with a unit cost function and satisfying $l = |V| = t$.

*Proof.* We reduce HPP to an instance of `FLATSP` where the cost function $c$ simply counts the edges in the path, and where the path length $l$ equals the threshold $t$ and the number of nodes in $V$.

Let $G = (V, A)$ be a directed graph, which is an instance of HPP. Let $H = (V, A, c: A \to \{1\})$, and $l := |V|$ and $t := l$. Thus $(H, l, t)$ is an instance of `FLATSP`.

Let us now show that there is an equivalence between the existence of a Hamiltonian path in $G$ and the existence of a path $p = (v_1, \dots, v_l)$ of $H$ such that $c(p) \le t$. Assume that $G$ has a Hamiltonian path $p$. In this case, $p$ is also a path in $H$ of length $|V|$, and then the cost of $p$ equals its length, i.e. $c(p) = \sum_{i=1}^{|V|} 1 = |V|$. Hence, there exists a path $p$ of $H$ such that $c(p) \le t = |V|$.

Assume that there exists a path $p = (v_1, \dots, v_{|V|})$ of $H$ such that $c(p) \le t$. As $p$ is a path it has no repeated
nodes, and as by assumption $l = |V|$, one gets that $p$ is a Hamiltonian path of $H$, and thus also a Hamiltonian path of $G$, since $G$ and $H$ share the same set of nodes and edges. $\square$

The second intermediate problem is called the *Read Graph Mapping Problem* (GRMP) and is defined below. It formalizes the mapping on a general sequence graph. Hence, DBGRMP is a specialization of GRMP, since it considers the case of the de Bruijn graph.

**Definition 7 (Graph Read Mapping Problem).** Given

* a directed graph $G = (V, A, x)$, whose edges are labeled by symbols of the alphabet ($x : A \to \Sigma$),

* $q := q_1 \dots q_{|q|}$ a word of $\Sigma^*$,

* a cost function $F: \Sigma \times \Sigma \to \mathbb{N}$,

* a threshold $t \in \mathbb{N}$,

GRMP consists in deciding whether there exists a path $p := (v_1, \dots, v_{|q|+1})$ of $G$ composed of $|q| + 1$ nodes, which generates a word $m := m_1 \dots m_{|q|} \in \Sigma^{|q|}$ such that $m_i := x((v_i, v_{i+1}))$, and which satisfies $\sum_{i=1}^{|q|} F(m_i, q_i) \le t$. Here, $m$ is called the word generated by $p$.

**Proposition 2.** *GRMP is NP-complete.*

*Proof.* We reduce FLATSP to GRMP.

Let $(G = (V, A, c : A \to \mathbb{N}), l \in \mathbb{N}, t \in \mathbb{N})$ be an instance of FLATSP. Let $\Sigma = \{y_1, \dots, y_{|\Sigma|}\}$ an alphabet larger than the largest value of $c(A)$, and let $s$ be the application such that $s: \{0, \dots, |\Sigma|\} \to \Sigma$ and such that for each $i$ in $\{0, \dots, |\Sigma|\}$, $s(i) = y_i$. Let $H = (V, A, x := s \circ c)$ and let $\alpha$ be a letter that does not belong to $\Sigma$, let $q = \alpha^{l-1}$ and $F$ such that for each $i$ in $\{0, \dots, |\Sigma|\}$, $F(\alpha, y_i) = i$. Thus, we obtain $|q| = l - 1$.

Now, let us show that there is an equivalence between the existence of a path $p = (v_1, \dots, v_l)$ of $G$ such that $c(p) \le t$ and the existence of a path $p' = (u_1, \dots, u_{|q|+1})$ of $H$ composed of $|q| + 1$ nodes, which generates a word $m = m_1 \dots m_{|q|}$ of $\Sigma^{|q|}$, where each $m_j = x((u_j, u_{j+1}))$, and such that $\sum_{j=1}^{|q|} F(m_j, q_j) \le t$. Assume that there exists a path $p = (v_1, \dots, v_l)$ of $G$ such that $c(p) \le t$. By definition, $p$ is a path in $H$. Let $m$ be the word generated by $p$. Thus we have $\sum_{j=1}^{|q|} F(m_j, q_j) = \sum_{j=1}^{l-1} F(m_j, \alpha) = \sum_{j=1}^{l-1} c((v_j, v_{j+1})) \le t$.

Now, suppose that there exists a path $p' = (u_1, \dots, u_{|q|+1})$ of $H$ composed of $|q| + 1$ nodes, which generates a word $m = m_1 \dots m_{|q|}$ of $\Sigma^{|q|}$, where each $m_j = x((u_j, u_{j+1}))$, and such that $\sum_{j=1}^{|q|} F(m_j, q_j) \le t$. By the construction of $H$, $p'$ is a path in $G$ of length $|q| + 1 = l$. Hence, we obtain $\sum_{j=1}^{l-1} c((u_j, u_{j+1})) = \sum_{j=1}^{|q|} F(m_j, \alpha) = \sum_{j=1}^{l-1} F(m_j, q_j) \le t$. $\square$

**Theorem 1.** *DBGRMP is NP-complete.*

Figure 1 illustrates the gadget used in the proof of Theorem 1. Basically, the gadget creates a **DBG** node (a word) formed by concatening the labels of the two preceding edges in the original graph.

*Proof.* Let us now reduce GRMP to **DBGRMP**.

Let $(G := (V, A, x : A \to \Sigma), q \in \Sigma^*, F : \Sigma \times \Sigma \to \mathbb{N}, t \in \mathbb{N})$ be an instance of GRMP. Let $\$$ and $\Delta$ be two distinct letters that do not belong to $\Sigma$, and let $\Sigma' := \Sigma \cup \{\$, \Delta\}$. Let $V'$ be a set of words of length 2 defined by

$$ V' := \left\{ \begin{array}{l} \{\alpha_i \beta_j | x(i,j) \\ \qquad = \alpha \text{ and } \exists l \in V \text{ such that } x(j, l) = \beta\} \end{array} \right\} \quad \text{set 1} \\ \cup \quad \left\{ \begin{array}{l} \{\Delta_i \$_i | \exists j \in V, \text{ such that } x(i,j) \\ \qquad = \alpha \text{ and } \nexists l \in V \text{ such that } (l, i) \in A\} \end{array} \right\} \quad \text{set 2} \\ \cup \quad \left\{ \begin{array}{l} \{\$_i \alpha_i | \exists j \in V, \text{ such that } x(i,j) \\ \qquad = \alpha \text{ and } \nexists l \in V \text{ such that } (l, i) \in A\} \end{array} \right\} \quad \text{set 3} \end{aligned} \right. \quad (1) $$

Any letter of a word in $V'$ is a symbol of $\Sigma'$ numbered by a node of $V$. Moreover, if that symbol is taken from $V$ then it labels an edge of $A$ that goes out a node, say $i$, of $V$, and the number associated to that symbol is $i$. In fact, $V'$ is the union of three sets (see Eq. 1):

set 1 considers the cases of an edge of $A$ labeled $\alpha$ followed by an edge labeled $\beta$, sets 2 and 3 contain the cases of an edge of $A$ labeled $\alpha$ that is not preceded by another edge of $A$; for each such edge one creates two words: $\Delta_i \$_i$ in set 2 and $\$_i \alpha_i$ in set 3.

Let $H$ be the 2-dBG of $V'$; note that $\Sigma'$ is the alphabet of the words of $V'$. Now let $z$ be the application from $V'$ to $\Sigma$ that for any $\alpha_i$ of $V'$ satisfies $z(\alpha_i) = \alpha$. (Note that in this equation, the right term is a shortcut meaning the symbol of $\alpha_i$ without its numbering $i$; this shortcut is used only for the sake of legibility, but can be properly written with a heavier notation). Let $F' : \Sigma' \times \Sigma \to \mathbb{N}$ be the application such that $\forall (\alpha_i, \beta) \in \Sigma' \times \Sigma$, $F'(\alpha_i, \beta) = F(z(\alpha_i), \beta) = F(\alpha, \beta)$.

Let us show that this reduction is a bijection that transforms a positive instance of GRMP into a positive instance of **DBGRMP**. Assume there exists a path $p := (v_1, \dots, v_{|q|+1})$ of $G$ which generates a word $m = m_1 \dots m_{|q|} \in \Sigma^{|q|}$ satisfying $m_i = x((v_i, v_{i+1}))$ and such that $\sum_{i=1}^{|q|} F(m_i, q_i) \le t$. We show that there exists a path $p'$ of $G'$ which generates a word $m' = m'_1 \dots m'_{|q|} \in \Sigma'^{|q|}$ such that $\sum_{i=1}^{|q|} F(m'_i, q_i) \le t$.

We build the path $p'$ as the “concatenation” of two paths, denoted $p'_{start}$ and $p'_{end}$, that we define below. Let $\gamma_j := x((v_j, v_{j+1}))_{v_j} = (m_j)_{v_j}$ for all $j$ between 1 and $|q|$. One has that $\gamma_j \in \Sigma'$. Now, let
<figure><img src="image_57.png" /><figcaption>**Fig. 1** Illustration of the gadget used in the proof of Theorem 1. Encoding a directed graph into a CDBG of order 2. The directed graph G (top) admits the same words than the 2-DBG G' (bottom), if we ignore the numbers</figcaption></figure>

$$p'_{start} := \begin{cases} (x((v'_l, v_l))_{v'_l}, x((v_l, v_1))_{v'_l}, x((v_l, v_1))_{v'_l}, x((v_1, v_2))_{v_1}) & \text{if } \exists l, l' \in V \text{ such that } (l, 1) \in A \text{ and } (l', l) \in A \\ (\$v_l, x((v_l, v_1))_{v'_l}, x((v_l, v_1))_{v'_l}, x((v_1, v_2))_{v_1}) & \text{if } \exists l \in V \text{ such that } (l, 1) \in A \text{ and } \nexists l' \in V \text{ such that } (l', l) \in A \\ (\Delta_{v_1}, \$v_1, \$v'_l, x((v_1, v_2))_{v_1}) & \text{otherwise.} \end{cases}$$

and let

$$p'_{end} := (\gamma_1\gamma_2, \dots, \gamma_{|q|-1}\gamma_{|q|}).$$

Let $m'$ denote the word generated by $p'$. Clearly, one sees that $m' = (m_1)_{v_1} \dots (m_{|q|})_{v_{|q|}}$, and since $m_i = z((m'_i)_{v_i})$, one gets that $z(m') = m$ and $\sum_{i=1}^{|q|} F'(m'_i, q_i) = \sum_{i=1}^{|q|} F(m_i, q_i) \le t$.

In the other direction, the proof is similar since our construction is a bijection. □

## GGMAP: a method to map reads on de Bruijn Graph

We propose a practical solution for solving DBGRMP. We consider the case of short (hundred of base pairs) reads with a low error rate (1 % of substitution), which is a good approximation of widely used NGS reads. Since errors are mostly substitutions, mapping is computed using the Hamming distance.

Our solution is designed for mapping on a compacted de Bruijn graph (CDBG) any set of short reads, either those used to build the graph or reads from another individual or species. We recall that a CDBG is representation of a CDBG in which each non branching path is merged into a single node. The sequence of each node is called a *unitig*. Figure 2 shows a CDBG and the associated CDBG.

In a CDBG, the nodes are not necessarily k-mers, words of length k, but *unitigs*, with some unitigs being longer than reads. Thus, while mapping on a CDBG, one distinguishes between two mapping situations: **i/** the reads mapping completely on a unitig of the graph, and **ii/** the reads whose mapping spans two or more unitigs. For the latter, we say that the read *maps on a branching path of the graph*.

Taking advantage of the extensive research carried out for mapping reads on flat strings, GGMAP uses Bowtie2 [7] to map the reads on the unitigs. In addition, GGMAP integrates our proposed new tool, called BGREAT, for mapping reads on branching paths of the CDBG. Figure 3 provides an overview of the pipeline.

GGMAP takes as inputs a query set of reads and a reference CDBG. To avoid including sequencing errors in the CDBG, we construct the reference CDBG after filtering out all k-mers whose coverage lies below a user-defined threshold c. This error removal step is a classical preprocessing step that is performed in k-mer based assemblers. The unitigs of the CDBG are computed using BCALM2 (the parallel version of BCALM [21]), using the k-mers having a coverage ≥ c. GGMAP uses such a set of unitigs as CDBG.

We now propose a detailed description of BGREAT.

## BGREAT: mapping reads on branching paths of the CDBG

As previously mentioned, BGREAT is designed for mapping reads on branching paths of a CDBG, using reasonable resources both in terms of time and memory. Our approach follows the usual “seed and extend” paradigm. More generally, the proposed implementation applies heuristic schemes, both regarding the indexing and the alignment phases.

### Indexing heuristic

We remind that our algorithm maps reads that span at least two distinct unitigs. Such mapped reads inevitably traverse one or more CDBG edge(s). In a CDBG, edges are represented by the prefix and suffix of size $k-1$ of each unitig. We call such sequences the *overlaps*. In order to limit the index size and the computation time, our algorithm indexes only overlaps that are later used as seeds. Those overlaps are good anchors for several reasons: they are long enough ($k-1$) to be selective, they cannot be shared by more than eight unitigs (four starting and four ending with the overlap), and a CDBG usually has a reasonable number of unitigs and then of overlaps. For
<figure><img src="image_58.png" /><figcaption>**Fig. 2** A toy example of a CDBG of order *k* with *k* = 4 (top) and its compacted version (bottom)</figcaption></figure>

instance, the CDBG in our experiment with human data has 70 million unitigs and 87 million overlaps for 3 billion *k*-mers). In our implementation, the index is a minimal perfect hash table indicating for each overlap the unitig(s) starting or ending with this (*k* − 1)-mer. Using a minimal perfect hash function limits the memory footprint, while keeping efficient query times (see Table 3).

**Read alignment**

Given a read, each of its *k* − 1-mers is used to query the index. The index detects which *k* − 1-mers represent an overlap of the CDBG. An example of a read together with the matched unitigs are displayed on Fig. 4. Once the overlaps and their corresponding unitigs have been computed, the alignment of the read is performed from left to right as presented in Algorithm 1. Given an overlap position *i* on the read, the unitigs starting with this overlap are aligned to the sequence of the read starting from position *i*. The best alignment is recorded. In addition, to improve speed, if one of the at most four unitigs ending with the same overlap is the next overlap detected on the

read, then this unitig is tested first, and if the alignment contains less mismatch than the user defined threshold, the other unitigs are not considered. Note that this optimization does not apply for the first and last overlaps of a read.

**Algorithm 1:** Greedy algorithm for mapping a read on multiple unitigs once the potential overlaps present in the read have been detected.

**Data:** Read *r*, Integer *n*
**for** *the n first overlaps of r do**

Find a path *begin* that maps the begin of *r*
**if begin found then**

**for the n last overlaps of r do**
Find a path *end* that maps the end of *r*
**if end found then**

Find (in a greedy way) a path *cover* that
map the read from *begin* to *end*
**if cover found then**

write path;

**return**

<figure><img src="image_59.png" /><figcaption>**Fig. 3** Unitig construction, as used in the proposed experiments (upper part of the figure) and GGMAP pipeline. Reads to be mapped can be distinct from reads used for building the graph. Long unitigs are unitigs longer than the reads. We remind that tools BCALM and BOWTIE2 are respectively published in [7, 21]</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>**Fig. 4** Representation of the mapping of a read (top sequence) on a CDBG, whose nodes are represented on lines 2, 3, and 4. (step 1) the overlaps of the graph that are also present in the read are found (here TACAC, GCTGC, and AGCTA, represented on line 1). (step 2) unitigs that map the beginning and the end of the read are found (those represented on line 2). (step 3) cover the rest of the read, guided by the overlaps (here with unitigs represented on lines 3 and 4)</figcaption></figure>

This mapping procedure is performed only if the two extremities of the read are mapped by two unitigs. The extreme overlaps of the read enables BGREAT to quickly filter out unmappable reads. For doing this, the first (resp. last) overlap of the read is used to align the read to the first (resp. last) unitig. Note that, as polymorphism exists between the read and the graph, some of the overlaps present on the read may be spurious. In this case the alignment fails, and the algorithm continues with the next (resp. previous) overlap. At most *n* alignment failures are authorized in each direction. If a read cannot be anchored neither on the left, nor on the right, it is considered as not aligned to the graph.

Note that the whole approach is greedy: given two or more possible choices, the best one is chosen and back-tracking is excluded. This results in a linear time mapping process, since each position in the read can lead to a maximum of four comparisons, and the algorithm continues as long as the cumulated number of mismatches remains below the user defined threshold. Because of heuristics, a read may be unmapped or wrongly mapped for any of the following reasons.

* All overlaps on which the read should map contain errors, in this case the read is not anchored or only badly anchored and thus not mapped.

* The *n* first or *n* last overlaps of the read are spurious, in this case the *begin* or *end* is not found and the read is not mapped. By default and in all experiments *n* = 2.

* The greedy choices made during the path selection are wrong.

We implemented *BGREAT* as a dependence-free tool in C++ available at github.com/Malfoy/BGREAT.

## Results

Beforehand we give details about the data sets (Subsection Data sets and CDBG construction), then we perform several evaluations of *GGMAP* and of *BGREAT*. First, we compare graph mapping to mapping on the contigs resulting from an assembly (Subsection Graph mapping vs assembly mapping). Second, we assess how many reads are mapped on branching paths vs on unitigs (Subsection Mapping on branching paths usefulness). Third, we

evaluate the efficiency of BGREAT in both terms of throughput and scalability (Subsection *GGMAP* perfor-mances), then assess the quality of the mapping itself (Subsection *GGMAP* accuracy). All *BGREAT* alignments were performed authorizing up to two mismatches.

There are very few published tools to compare *GGMAP* with. Indeed, we found only one published tool, called *BlastGraph* [19], which was designed for mapping reads on a DBG. However, on our simplest data set coming from the *E.coli* genome (see Table 1), *BlastGraph* crashed after ≈ 124 h of computation. Thus, *BlastGraph* was not further investigated here.

## Data sets and CDBG construction

For our experiments we used publicly available Illumina read data sets from species of increasing complexity: from the bacterium *E.coli*, the worm *C.elegans*, and from Human. Detailed information about the data sets are given in Additional file 1: Table S1 (identifiers, read length, read numbers, and coverages – from 70x to 112x–).

For each of these three data sets, we generated a CDBG using BCALM. From the *C.elegans* read set, we additionally generated an artificially complex graph, by using small *k* and *c* values (respectively 21 and 2). This particular graph, called *C.elegans_cpx*, contains lot of small unitigs. We used it to assess situations of highly complex and/or low quality sequencing data. The characteristics of the CDBG obtained on each of these data sets are given in Table 1.

## Graph mapping vs assembly mapping

We compared *GGMAP* to the popular approach consisting in mapping the reads to the reference contigs computed by an assembler. For testing this approach, for each of the three sets used, we first assembled them and then we mapped back the reads on the obtained set of contigs. We used two different assemblers, the widely used Velvet [22], and Minia [23], a memory efficient assembler based on Bloom filters. Finally, we used Bowtie2 for mapping the reads on the obtained contigs.

The results reported in Table 2 show that the number of reads mapped on assembled contigs is smaller than the one obtained with *GGMAP*. We obtained similar results in terms of number of reads mapped on the
<table><caption>**Table 1** CDBG used in this study</caption><thead><tr><th>CDBG Id</th><th>Reads Id</th><th>k</th><th>c</th><th>Number of unitigs</th><th>Mean length of unitigs</th></tr></thead><tbody><tr><td>E.coli</td><td>SRR959239</td><td>31</td><td>3</td><td>42,843</td><td>134</td></tr><tr><td>C.elegans_norm</td><td>SRR065390</td><td>31</td><td>3</td><td>1,627,335</td><td>93</td></tr><tr><td>C.elegans_cpx</td><td>SRR065390</td><td>21</td><td>2</td><td>8,273,338</td><td>34</td></tr><tr><td rowspan="2">Human</td><td>SRR345593</td><td rowspan="2">31</td><td rowspan="2">10</td><td rowspan="2">69,932,343</td><td rowspan="2">70</td></tr><tr><td>SRR345594</td></tr></tbody></table>

C.elegans_cpx and C.elegans_norm are two distinct graphs, constructed using the same read set from C.elegans genome. The suffixes *norm* and *cpx* respectively stand for "normal" (using *c* = 3 and *k* = 31) and for "complex" (using a low threshold *c* = 2 and small value *k* = 21)

assemblies yielded by Velvet and Minia (see Additional file 1: Table S2). Let us emphasize that on the Human dataset, **GGMAP** maps 22 additional percents of reads on the graph than Bowtie2 does on the assembly.

We notice that the more complex the graph, the higher the advantage of mapping on the CDBG. This is due to the inherent difficulty of assembling with huge and highly branching graphs. This is particularly prominent in the results obtained on the artificially complex *C.elegans_cpx* CDBG.

We also highlight that our approach is resource efficient compared to most assembly processes. For instance, Velvet used more than 80 gigabytes of memory to compute the contigs for the *C. elegans* data set with *k* = 31. On this data set, our workflow used at most 4 GB memory (during *k*-mer counting). In terms of throughput, using **BGREAT** and then Bowtie2 on long unitigs is comparable to using Bowtie2 on contigs alone. See Section **GGMAP** performances for more details about **GGMAP** performances.

### Mapping on branching paths usefulness

Mapping the reads on branching paths of the graph is not equivalent to simply mapping the reads on unitigs. Indeed, at least 13% of reads (mapping reads SRR959239 on the *E.coli* CDBG) and up to 66% of reads (mapping reads SRR065390 on *C.elegans_cpx* CDBG) map on the branching paths of the graph (see Fig. 5). These reads cannot be mapped when using only the set of unitigs as a reference. As expected, the more complex the graph, the larger the benefit of **BGREAT**'s approach. On the complex *C.elegans_cpx* graph, only 23% of reads can be

fully mapped on unitigs, while 89% of them are mapped by additionally using **BGREAT**. On a simpler graph as *C.elegans_norm* the gap is smaller, but remains significant (72 vs 93%). Complete mapping results are shown in Additional file 1: Table S3.

### Non reflexive mapping on a CDBG

The **GGMAP** approach is also suitable for mapping a distinct read set from the one used for constructing the CDBG. We mapped another read set from *C.elegans* (SRR1522085) on the *C.elegans_norm* CDBG. Results in this situation are similar to those observed when performing reflexive mapping (i.e., when mapping the reads used to construct this graph): among 89% of mapped reads, 15% were mapped on branching paths of the graph (See Fig. 5).

### **GGMAP** performances

Table 3 presents **GGMAP** time and memory footprints. It shows that **BGREAT** is very efficient in terms of throughput while using moderate resources. Presented heuristics and implementation details allow **BGREAT** to scale up to real-world instances of the problem, being able to map millions of reads per CPU hour on a Human CDBG with a low memory footprint. **BGREAT** mapping is parallelized and can efficiently use dozens of cores.

### **GGMAP** accuracy

To measure the impact of the read alignment heuristics, we forced the tool to explore exhaustively all potential alignment paths once a read is anchored on the graph. Results on the *E.coli* dataset show that the greedy approach is much faster than the exhaustive one (38 × faster), while the mapping capacity is little impacted: the overall number of mapped reads increases by only 0.03% with the exhaustive approach. We thus claim that the choice of the greedy strategy is a satisfying trade-off.

To further evaluate the **GGMAP** accuracy, we assess the recall and mapping quality in the following experiment. We created a CDBG from Human chromosome 1 (hg19 version). Thus, each *k*-mer of the chromosome appears in the graph. Furthermore, from the same sequence, we

<table><caption>**Table 2** Percentage of mapped reads, either mapping on contigs (here obtained thank to the Minia assembler) or mapping on CDBG with **GGMAP**</caption><thead><tr><th>Set</th><th>% mapped on contigs</th><th>% mapped on CDBG</th></tr></thead><tbody><tr><td>E.coli</td><td>95,57</td><td>97,16</td></tr><tr><td>C.elegans_norm</td><td>80,60</td><td>93,24</td></tr><tr><td>C.elegans_cpx</td><td>56,33</td><td>89,15</td></tr><tr><td>Human</td><td>63,16</td><td>85,70</td></tr></tbody></table>
<figure><img src="image_61.png" /></figure>

simulated reads with distinct error rates (0, 0.1, 0.2, 0.5, 1 and 2%). For each error rate value, we generated one million reads. We evaluated the *GGMAP* results by mapping the simulated reads on the graph. As the graph is error free, except in some rare cases due to repetitions, the differences between a correctly mapped read and the path it maps to in the graph occur at erroneous positions of the read. If this is not the case, we say that the read is not mapped at its optimal position. Among the error free positions of a simulated read, the number of mismatches observed between this read and the mapped path is called the “distance to optimal”. Results are reported in Table 4 together with the obtained recall (number of mapped reads over the number of simulated reads). Those results

show the limits of *BGREAT* while mapping reads from divergent individuals. With 2% of substitutions in reads, only 90.85% of the reads are perfectly mapped. Nevertheless, with this divergence rate, 97.28% of reads are mapped at distance at most one from optimum. With over 99% of perfectly mapped reads, these results show that with the current sequencing characteristics, i.e. a 0.1% error rate, the mapping accuracy of *BGREAT* is suitable for most applications.

## Discussion

We proposed a formal definition of the de Bruijn graph Read Mapping Problem (DBGRMP) and proved its NP-completeness. We proposed a heuristic algorithm offering

<table><caption>**Table 3** Time and memory footprints of *BGREAT* and *BOWTIE2*</caption><thead><tr><th rowspan="2">CDBG Id</th><th rowspan="2">Mapped set<br>(nb reads)</th><th colspan="3">BGREAT</th><th colspan="3">BOWTIE2</th></tr><tr><th>Wall clock time</th><th>CPU time</th><th>Memory</th><th>Wall clock time</th><th>CPU time</th><th>Memory</th></tr></thead><tbody><tr><td>E.coli</td><td>SRR959239<br>(5,128,790)</td><td>28 s</td><td>1m40</td><td>19 MB</td><td>1m17</td><td>3m53</td><td>29 MB</td></tr><tr><td>C.elegans_cpx</td><td>SRR065390<br>(67,155,743)</td><td>19m21</td><td>72m31</td><td>975 MB</td><td>8m12</td><td>33m</td><td>1.66 GB</td></tr><tr><td>C.elegans_norm</td><td>"</td><td>13m03</td><td>51m28</td><td>336 MB</td><td>17m49</td><td>72m31</td><td>493 MB</td></tr><tr><td>C.elegans_norm</td><td>SRR1522085<br>(22,509,110)</td><td>1m54</td><td>7m13</td><td>336 MB</td><td>3m29</td><td>14m12</td><td>493 MB</td></tr><tr><td>Human</td><td>SRR345593<br>SRR345594<br>(2,967,536,821)</td><td>4h30</td><td>87 h</td><td>9.7 GB</td><td>4h38</td><td>90h15</td><td>21 GB</td></tr></tbody></table>

Indicated wall clock times use four cores, except for the human samples for which 20 cores were used
<table><caption>**Table 4** *GGMAP* mapping results on simulated reads from the reference of the human chromosome 1 with default parameters</caption><thead><tr><th rowspan="2">% Errors in simulated reads</th><th colspan="5">Distance to optimum of BGREAT mapped reads (percentage)</th></tr><tr><th>0</th><th>1</th><th>2</th><th>3</th><th>≥ 4</th></tr></thead><tbody><tr><td>0</td><td>100</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0.1</td><td>99.31</td><td>0.52</td><td>0.09</td><td>0.04</td><td>0.04</td></tr><tr><td>0.2</td><td>98.79</td><td>0.91</td><td>0.21</td><td>0.07</td><td>0.02</td></tr><tr><td>0.5</td><td>97.2</td><td>2.17</td><td>0.41</td><td>0.17</td><td>0.05</td></tr><tr><td>1</td><td>94.88</td><td>3.72</td><td>0.92</td><td>0.41</td><td>0.07</td></tr><tr><td>2</td><td>90.85</td><td>6.43</td><td>1.79</td><td>0.83</td><td>0.1</td></tr></tbody></table>

Results show the recall of *GGMAP* and the quality of *BGREAT* mapping, as represented by the “distance to optimum” value. For instance 94.88 % of the reads were mapped without error, 3.72 % were mapped with a distance to the optimum of one etc. Due to approximate repeats in human chromosome 1, the reported distance to optimum is an upper bound

a practical solution. We developed a tool called *BGREAT* implementing this algorithm using a compacted de Bruijn graph (CDBG) as a reference.

From the theoretical viewpoint, the problem DBGRMP considers paths rather than walks in the graph. The current proof of its hardness does not seem to be adaptable to the cases of walks. A perspective is to extend the hardness result to that more general case.

We emphasize that our proposal does not enable genome annotation. It has been designed for applications aiming at a precise quantification of sequenced data, or a set of potential variations between the reads and the reference genome. In this context, it is essential to map as much reads as possible. Experiments show that a significant proportion of the reads (between ≈ 13 and ≈ 66% depending on the experiment) can be only mapped on branching paths of the graph. Hence, mapping only on the nodes of the graph or on assembled contigs is thus insufficient. This statement holds true when mapping the reads used for building the graph, but also with reads from a different experiment. Moreover, our results show that a potentially large number of reads (up to ≈ 32%) that are mapped on a CDBG cannot be mapped on a classical assembly.

With *GGMAP*, the mapping quality is very high: using Human chromosome 1 as a reference and reads with a realistic error rate (similar to that of Illumina technology), over 99 % of the reads are correctly mapped. The same experiment also pointed out the limits of mapping reads on a divergent graph reference (≥ 2 % substitutions): approximately 10 % of the reads are mapped at a suboptimal position.

A weak point of *BGREAT* lies in its anchoring technique. Reads mapped with *BGREAT* must contain at least one exact *k* - 1-mer that is an arc of the CDBG, i.e., an overlap between two connected nodes. This may be a serious limitation when the original read set diverges greatly

from the reads to be mapped. Improving the mapping technique may be done by using not only unitig overlaps as anchors at the cost of higher computational resources. Another solution may consist in using a smarter anchoring approach, like spaced seeds, which can accommodate errors in the anchor [24].

A natural extension consists in adapting *BGREAT* for mapping, on the CDBG obtained from short reads, the long (a few kilobases in average) and noisy reads produced by the third generation of sequencers, whose error rate reaches up to 15 % (with mostly insertion and deletion errors for e.g. Pacific Biosciences technology). Such adaptation is not straightforward because of our seeding strategy, which requires long exact matches. The anchoring process must be very sensitive and very specific, while the mapping itself must implement a Blast-like heuristic or an alignment-free method. However, mapping such long reads on a CDBG could be of interest for correcting these reads as in [13], or for solving repeats, if long reads are mapped on the walks (which main include cycles) of the CDBG. Our NP-completeness proof only considers mapping on (acyclic) paths. Proving the hardness of the problem of mapping reads on walks of a CDBG remains open.

Incidentally, using the same read set for constructing the CDBG and for mapping opens the way to major applications. Indeed, the graph and the exact location of each read on it may be used for *i*/ read correction as in [15], by detecting differences between reads and the mapped area of the graph in which low support *k*-mers likely due to sequencing errors are absent, or for *ii*/ read compression by recording additionally the mapping errors, or for *iii*/ both correction and compression by conserving only for each read its mapping location on the graph.

Having for each read (used for constructing the graph or not) its location on the CDBG also provides the opportunity to design algorithms for enriching the graph, for instance enabling a quantification that is sensitive to local variations. This would be valuable for applications such as variant calling, analysis of RNA-seq variants [25], or of metagenomic reads [26].

Additionally, *BGREAT* results provide pieces of information for distant *k*-mers in the CDBG, about their co-occurrences in the mapped read data sets. This offers a way for the resolution, in the de Bruijn graph, of repeats larger than *k*. It could also allow to phase the polymorphisms and to reconstruct haplotypes.

## Conclusion

A take home message is that read mapping can be significantly improved by mapping on the structure of an assembly graph rather than on a set of assembled contigs (respectively ≈ 22 % and ≈ 32 % of additional reads mapped for the Human and a complex *C.elegans* data
sets). This is mainly due to the fact that assembly graphs retains more genomic information than assembled contigs, which also suffer from errors induced by the complexity of assembly. Moreover, mapping on a compacted De Bruijn Graph can be fast. The availability of BGREAT opens the door to its application to fundamental tasks such as read error correction, read compression, variant quantification, or haplotype reconstruction.

## Additional file 1

**Additional file 1:** Read mapping on De Bruijn graphs additional file. Three complementary tables are presented. Main characteristics of data sets used in this study. Assembly and mapping approach comparison. Results of BGREAT on real read sets. (PDF 40 kb)

### Abbreviations

CDBG, Compacted De Bruijn graph; DBG, De Bruijn graph; DBGRMP, De Bruijn graph read mapping problem; FLATSP, fixed length assymetric travelling salesman problem; GRMP, graph read mapping problem; HPP, Hamiltonian path problem

### Acknowledgements

We would like to thank Yannick Zakowski, Claire Lemaitre and Camille Marchet for proofreading the manuscript and discussions.

### Funding

This work was funded by French ANR-12-B502-0008 Colib'read project, by ANR-11-BINF-0002, and by a MASTODONS project.

### Availability of data and materials

Our implementations are available at github.com/Malfoy/BGREAT. In addition to the following pieces of information, Additional file 1: Table S1 presents the main characteristics of these datasets.

SRR959239 http://www.ncbi.nlm.nih.gov/sra/?term=SRR959239
SRR065390 http://www.ncbi.nlm.nih.gov/sra/?term=SRR065390
SRR1522085 http://www.ncbi.nlm.nih.gov/sra/?term=SRR1522085
SRR345593 and SRR345594 http://www.ncbi.nlm.nih.gov/sra/?term=SRR345593.

### Authors' contributions

PP initiated the work and designed the study. AL, BC and ER designed the formalism and the proofs of NP-hardness. AL designed the algorithmic framework, implemented the BGREAT and performed the tests. All authors wrote and accepted the final version of the manuscript.

### Competing interests

The authors declare that they have no competing interests.

### Ethics approval and consent to participate

Not applicable.

### Author details

¹IRISA Inria Rennes Bretagne Atlantique, GenScale team, Campus de Beaulieu, 35042 Rennes, France. ²L.I.R.M.M., UMR 5506, Université de Montpellier et CNRS, 860 rue de St Priest, F-34392 Montpellier Cedex 5, France. ³Institut Biologie Computationnelle, Université de Montpellier, F-34392 Montpellier, France.

Received: 9 December 2015 Accepted: 26 May 2016

Published online: 16 June 2016

### References

1. Bradnam KR, Fass JN, et al. Assemblathon 2: evaluating de novo methods of genome assembly in three vertebrate species. GigaScience. 2013;2:10. doi:10.1186/2047-217X-2-10.

2. Nagarajan N, Pop M. Parametric complexity of sequence assembly: theory and applications to next generation sequencing. J Comput Biol. 2009;16(7):897–908. doi:10.1089/cmb.2009.0005.

3. Chaisson MJ, Pevzner PA. Short read fragment assembly of bacterial genomes. Genome Res. 2008;18(2):324–30. doi:10.1101/gr.7088808.

4. Myers EW, Sutton GG, et al. A whole-genome assembly of Drosophila. Science (New York, N.Y.) 2000;287(5461):2196–204. doi:10.1126/science.287.5461.2196.

5. Myers EW. The fragment assembly string graph. Bioinformatics. 2005;21(Suppl 2):79–85. doi:10.1093/bioinformatics/bti1114.

6. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60. doi:10.1093/bioinformatics/btp324.

7. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9. doi:10.1038/nmeth.1923.

8. Lee H, Schatz MC. Genomic dark matter: the reliability of short read mapping illustrated by the genome mappability score. Bioinformatics. 2012;28(16):2097–105. doi:10.1093/bioinformatics/bts330.

9. Deshpande V, Fung EDK, Pham S, Bafna V. Cerulean: A Hybrid Assembly Using High Throughput Short and Long Reads. In: Lecture Notes in Computer Science vol. 8126 LNBI. Springer; 2013. p. 349–63. doi:10.1007/978-3-642-40453-5_27.

10. Ribeiro FJ, Przybylski D, Yin S, Sharpe T, Gnerre S, Abouelleil A, Berlin AM, Montmayeur A, Shea TP, Walker BJ, Young SK, Russ C, Nusbaum C, MacCallum I, Jaffe DB. Finished bacterial genomes from shotgun sequence data. Genome Res. 2012;22(11):2270–7. doi:10.1101/gr.141515.112.

11. Pevzner PA, Tang H, Waterman MS. An Eulerian path approach to DNA fragment assembly. Proc Natl Acad Sci. 2001;98(17):9748–53. doi:10.1073/pnas.171285098.

12. Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, Lesin VM, Nikolenko SI, Pham S, Prjibelski AD, Pyshkin AV, Sirotkin AV, Vyahhi N, Tesler G, Alekseyev MA, Pevzner PA. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19(5):455–77. doi:10.1089/cmb.2012.0021.

13. Salmela L, Rivals E. LoRDEC: accurate and efficient long read error correction. Bioinformatics. 2014;30(24):3506–14. doi:10.1093/bioinformatics/btu538.

14. Yang X, Chockalingam SP, Aluru S. A survey of error-correction methods for next-generation sequencing. Brief Bioinform. 2013;14(1):56–66. doi:10.1093/bib/bbs015.

15. Benoit G, Lavenier D, Lemaitre C, Rizk G. Bloocoo, a memory efficient read corrector. In: European Conference on Computational Biology (ECCB); 2014. https://gatb.inria.fr/software/bloocoo/.

16. Wang M, Ye Y, Tang H. A de Bruijn graph approach to the quantification of closely-related genomes in a microbial community. J Comput Biol. 2012;19(6):814–25. doi:10.1089/cmb.2012.0058.

17. Huang L, Popic V, Batzoglou S. Short read alignment with populations of genomes. Bioinformatics. 2013;29(13):361–70. doi:10.1093/bioinformatics/btt215.

18. Dilthey A, Cox C, Iqbal Z, Nelson MR, McVean G. Improved genome inference in the MHC using a population reference graph. Nat Genet. 2015;47(6):682–8. doi:10.1038/ng.3257.

19. Holley G, Peterlongo P. Blastgraph: Intensive approximate pattern matching in sequence graphs and de-bruijn graphs. In: Stringology; 2012. p. 53–63. http://alcovna.genouest.org/blastree/.

20. Karp RM. Reducibility Among Combinatorial Problems. In: 50 Years of Integer Programming 1958-2008. Berlin, Heidelberg: Springer; 2010. p. 219–41. doi:10.1007/978-3-540-68279-0_8. http://link.springer.com/10.1007/978-3-540-68279-0_8.

21. Chikhi R, Limasset A, Jackman S, Simpson JT, Medvedev P. On the representation of de bruijn graphs. In: Research in Computational Molecular Biology. Springer; 2014. p. 35–55. doi:10.1007/978-3-319-05269-4-4.

22. Zerbino DR, Birney E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18(5):821–9. doi:10.1101/gr.074492.107.

23. Chikhi R, Rizk G. Space-efficient and exact de Bruijn graph representation based on a Bloom filter. Algorithm Mol Biol. 2013;8(1):22. doi:10.1186/1748-7188-8-22.

24. Vroland C, Salson M, Touzet H. Lossless seeds for searching short patterns with high error rates. In: Combinatorial Algorithms. Springer; 2014. p. 364–75.
25. Sacomoto GA, Kielbassa J, Chikhi R, Uricaru R, Antoniou P, Sagot MF, Peterlongo P, Lacroix V. Kissplice: de-novo calling alternative splicing events from rna-seq data. BMC Bioinformatics. 2012;13(Suppl 6):5.

26. Ye Y, Tang H. Utilizing de Bruijn graph of metagenome assembly for metatranscriptome analysis. Bioinformatics. 2015btv510. Oxford Univ Press. arXiv preprint arXiv:1504.01304.

Submit your next manuscript to BioMed Central
and we will help you at every step:

* We accept pre-submission inquiries

* Our selector tool helps you to find the most relevant journal

* We provide round the clock customer support

* Convenient online submission

* Thorough peer review

* Inclusion in PubMed and all major indexing services

* Maximum visibility for your research

Submit your manuscript at
www.biomedcentral.com/submit

**BioMed Central**
Chapter 4

Improving an assembly paradigm
In this chapter we are interested in improving the de Bruijn graph assembly. We first describe the limitations of actual assembly methods and propose a high order de Bruijn graph as a scalable and efficient solution. Then we propose a proof of concept method to construct such graph based on successive mapping. In the last part, we present different results of our proof of concept assembler and show its potential to tackle several existing limitations.

# 4.1 Assembly challenges

## 4.1.1 Two assembly paradigms

As we previously saw, OLC and string graph approaches are able to use the full read information. They allow better repeats handling than the de Bruijn graph and produce more contiguous assemblies. From a qualitative point of view, the best assemblers should be those based on string graphs, as assessed by the Assemblathon assembly competition [122] [123]. However, in the introduction, we presented those approaches as intractable because of the need to index the very large reads sets and to compute pairwise alignments. Large genomes are costly to assemble, and alternative approaches have been preferred for scalability reasons. Contrary to string graphs, despite being tractable, the de Bruijn graph do not use full read information and produce fragmented assembly. Moreover, high order de Bruijn graphs are hard to obtain because they require a high sequencing depth.

## 4.1.2 Why read length matters: polyploid assembly

The major practical difference between the de Bruijn graph and the string graph is thus the sequence length information they use. Having full length helps solve repeats, but since read lengths are bounded and particularly short for NGS, the theoretical edge of the string graph over the de Bruijn graph is not major. De Bruijn graph assemblers are not really behind string graph in practice, they can use an order above 60 from a high coverage, when the reads length goes up to a few hundred. Then, the advantages of the string graphs will come in the range of repeats between *k* and the read length. We mentioned in the introduction that *C. elegans* genome presents 1,892,518 repeated 61-mers and 824,254 repeated 301-mers. Even a string graph using 300 bp reads could not solve those repeats. But in some cases the number of short repeats can become extremely high, for instance in diploid genomes. As mentioned before, the variations between two alleles create bubble patterns in assembly graphs. Yet, according to the heterozygosity rate of a given genome, the three following situations may be distinguished.

* Low heterozygosity rate (< 0.1%): most variant are isolated from the others and each appears as a bubble in the graph (Figure 3.3). Most nodes of the graph represent homozygous regions of the underlying genome.

* High heterozygosity rate (> 5%): the divergence between the sequences is so high that the two alleles are mostly separated in the assembly graph (for example with a high *k* as seen in the Figure 4.1). Nodes of the graph mostly represent heterozygous
regions, they yield large bubbles that relatively well represent and differentiate the diverging alleles.

* Medium heterozygosity rate (≈ 1%): present a mix of heterozygous and homozygous regions that may be difficult to distinguish. It results into complex graphs, hard to produce long contigs in practice.

The heterozygosity rates given are rather orders of magnitude than accurate values (Figures 4.5 and 4.6). Furthermore, these rates are based on the use of hundred of bases NGS reads. If we were able to produce longer reads, we could use a higher *k* and a heterozygosity rate of 1% would be easier to manage.

Heterozygosity rates vary a lot across species as seen in Figure 4.2 from [124]. From our current knowledge, one of the lowest rates is present in Lynx individuals that have a rate of the order of 0.01%, where the heterozygous events are likely to be far apart. On the opposite, in species like *C. remanei*, with a rate of 5%, the variations should be quite close to each other on average and therefore a large order of Bruijn graph would be able to differentiate two alleles in the heterozygous regions. Most of Bruijn graph assemblers present a "comfort zone" below the rate of 0.5%, where they are able to produce large contigs by crushing bubbles. Medium rates (around 1%) result in extremely poor assemblies. As mentioned, very high rates (> 5%) genomes are manageable, although we currently count few of them among known assembled species.

### 4.1.3 Taking the best of both worlds

In any case, haploid or polyploid, the faculty to solve repeats is the key to more contiguous assemblies. Even if this problem is dramatically more frequent in polyploid assembly, enhancing this capacity would positively impact any assembly problem. We explained that each presented graph structure shows advantages as well as drawbacks. The question we should ask is "What kind of graph structure would we ideally want ?"

<figure><img src="image_62.png" /><figcaption>Figure 4.1: Example of phased and non-phased variants depending on *k*. With *k* = 4 the graph presents two bubbles, and the assembler has no mean to output phased contigs. With *k* = 6 the variants are phased and the graph is reduced into a single heterozygous region of two contigs.</figcaption></figure>
<figure><img src="image_63.png" /><figcaption>Figure 4.2: Heterozygosity rates across species [124].</figcaption></figure>
<figure><img src="image_64.png" /><figcaption>Figure 4.3: The general idea to improve the de Bruijn graph using read length. The reads information is used here in order to "solve" the gray repeat, *ie* to determine the possible contexts of this repeats into contigs. Here two contexts exist according to reads, the orange-gray-orange contig an the green-gray-green.</figcaption></figure>

Preferably we look for a graph:

* Using the whole read information, in order to get events phased in the reads also phased in the graph (not a de Bruijn graph)

* Easy to compute (not a string graph)

While several contributions improved the string graphs scalability, the de Bruijn graph assemblers main drawbacks was also subject to improvements. The fact that most people used de Bruijn graph assemblers led to new methods to improve the repeat resolution of such graphs, in order to produce less fragmented assemblies [39, 125]. The global idea is represented in the Figure 4.3. We have a graph made of trusted sequences, built from reads kmers. We are interested in adding back the read information in order to simplify some situations. The questions about how to make the best use of such information (how to produce longer contigs according the read information) and how to obtain it (how to map reads on the graph) has not reached a consensus yet.

A simple solution would be to construct a very high order de Bruijn graph with a *k* value near the size of the reads. But such a graph would require a tremendous coverage as shown in Table 4.1.

This raises the question: "How to construct a high order de Bruijn graph without a massive coverage ?" Several assemblers have tried to improve their repeat resolution power by constructing a high order De Bruin graph. Tools like IDBA [39] or SPAdes [36] tackle this question by using multiple kmer sizes approaches. The idea is to use several de Bruijn graph with different orders, to be able to detect all overlaps with the low *k* values, and solve some repeats with a high *k* value. Each of those techniques have shown how to improve assemblies in a consistent way, even if they present drawbacks that may be addressed. In the following we propose a new way to use the multiple kmer idea to
<table><thead><tr><th>Coverage</th><th>k=240</th><th>k=230</th><th>k=220</th></tr></thead><tbody><tr><td>100X</td><td>67</td><td>44</td><td>27</td></tr><tr><td>100X filter</td><td>94</td><td>81</td><td>63</td></tr><tr><td>200X</td><td>45</td><td>18</td><td>6</td></tr><tr><td>200X filter</td><td>81</td><td>50</td><td>24</td></tr><tr><td>500X</td><td>14</td><td>1</td><td>0</td></tr><tr><td>500X filter</td><td>42</td><td>8</td><td>1</td></tr><tr><td>1000X</td><td>2</td><td>0</td><td>0</td></tr><tr><td>1000X filter</td><td>9</td><td>0</td><td>0</td></tr></tbody></table>

Table 4.1: Percentage of missing kmers in a graph created with 250bp reads simulated from a small virus genome of 48,502 base pairs called Lambda phage with 1% error rate. The filtering consisted in removing kmers with an abundance of one. Even a tremendous coverage does not allow to construct directly a high order de Bruijn graph.

build a very high order de Bruijn graph based on the techniques previously proposed in this document.

### "Assembly challenges" core messages:

* Using the whole length of the reads can reduce assembly fragmentation

* This is especially important in polyploid assembly

* A high order de Bruijn graph would be a scalable structure very similar to the string graph

## 4.2 Building very high order de Bruijn graph

In this section we propose a new method using the previously described tools able to construct very high order de Bruijn graphs without the need of a costly coverage.

### 4.2.1 Successive mapping strategy

We have shown previously that the mapping of the reads on the de Bruijn graph produced very accurate corrected reads. Such reads with their lower error rates could be used in order to build a higher order de Bruijn graph. However we can see unitigs as trusted genome sequences that we want to order. Our hypothesis is that a read mapping on a path of unitigs validates their ordering. It means that a mapped read validates the whole path and not only the nucleotides where the read actually maps. As the Masurca assembler [126], we call super-read (SR) the paths of unitigs that correspond to a read mapping. An interesting property of a SR is that its length is superior (or equal) to its associated read because it may be extended by the first and last unitig sequences (Figure 4.4). Unlike unitigs, that can be smaller than reads in complex regions, the SR are guaranteed to be at least the size of the reads, and most of the time larger.
<figure><img src="image_65.png" /><figcaption>Figure 4.4: Super-reads generated by mapping the reads on the de Bruijn graph. For example the orange read maps on the unitigs 1, 3 and 4, its associated super-read is therefore the concatenation of those unitigs. We can see that super-reads may be longer than their associated reads, for example the orange read "AGCACTCATGCATGCGGTAT" is associated to the super-read "ACATGCATGCTAGCACTCATGCATGCGGTATATATATATATG". Several reads mapping on the same path generate the same super-read as the two pink reads, and some super-reads may be included in another super-read. Super-reads that are not included in any other super-read are called maximal super-reads.</figcaption></figure>
<figure><img src="image_66.png" /></figure>

Figure 4.5: How the de Bruijn graph unitigs can enlarge the overlaps. In this example the two reads in blue share a small overlap. They are mapped on the de Bruijn graph and they both map on the gray unitig. Their SR in purple, enlarged with the unitig sequence, share a larger overlap than their associated reads.

To assess the quality of the SR obtained by mapping on a clean de Bruijn graph we conducted the following experiment.

* Simulation of 100X reads from a reference genome

* Read correction with Bloocoo

* de Bruijn graph construction with $k=51$

* Tip removal

* Read mapping on the graph to generate SR

* Duplicates removal from SR

* SR mapping on the reference genome

On the *E. coli* genome 1,642 different SR were obtained, all mapped on the reference genome with a similarity of 99%. All SR but 5 mapped with a similarity of 100% on the reference. On the *C. elegans* genome 357,183 different SR were obtained. All SR but 38 were mapped on the reference with a similarity of 95%. With a higher similarity of 99%, the number of unmapped SR was only of 260 and all SR but 502 mapped perfectly. Those experiments show that SR are very accurate sequence that can be almost considered as genome sequences.

Our idea is to use the sequences from such SR in order to produce higher order de Bruijn graph (Figure 4.9). By constructing a de Bruijn graph from the SR with an augmented *k*, we obtain a simpler de Bruijn graph than the first one constructed, with longer unitigs, less repeats and less spurious edges. As in the classical multiple *k* approaches like IDBA, the interest of starting with a low *k* is to be able to detect small overlaps between reads. Small overlaps are detected with the first value of *k* and kept in the latter graph with higher *k* because the overlaps between the SR are larger than the overlaps between reads. (Figure 4.5). However, if most of the time the overlaps are enlarged by the unitigs graph, it is possible that the augmented overlap is not large enough to be kept in the later steps.

We can raise *k* more gradually in order for the unitigs of the overlaps to be larger, along with the simplification of the graph, to avoid losing overlaps. Furthermore some
<figure><img src="image_67.png" /><figcaption>Figure 4.6: How the graph cleaning at each step helps the mapping of reads. In this example, two reads in blue are mapped on the de Bruijn graph unitigs. They share a small overlap that can not be augmented because the graph is branching due to a sequencing error. The tip removal allows to get a large unitig and two SR in purple that share a large overlap. If a too high value of *k* was used directly to build the de Bruijn graph, the overlaps between the two reads would have been lost.</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>Figure 4.7: This example shows the interest of using several value of *k*. In this situation, the reads A and B share a small green overlap and the reads B and C share a large one with the reads Y and Z (in green and orange respectively). However, these overlaps share a medium repeat in brown. The low *k* assembly is fragmented because of the repeat but is able to detect the small overlap. The high *k* graph is able to avoid the repeat fragmentation but misses the small overlap. The successive construction allows to take advantages of both *k* values.</figcaption></figure>

errors would become tips since the higher the *k*, the more sequencing errors become tips, so the graph can be cleaned in a better way at each step. Then the removal of such errors would also allow the generation of more contiguous unitigs, and would help in mapping reads with a higher *k* as it is shown in Figure 4.6. The Figure 4.7 shows the interest of such an iterative approach. We propose a workflow based on the presented ideas as a proof of concept assembler called BWISE. We give a high level view of the BWISE workflow and a more detailed version is presented in Figure 4.8.

* Read correction

* de Bruijn graph construction

* Graph tipping

* Read mapping on the de Bruijn graph

* Super-Reads selection

* Higher de Bruijn graph construction
<figure><img src="image_69.png" /><figcaption>Figure 4.8: Bwise workflow. Data are represented in blue and software in green. The initial reads are corrected and a first de Bruijn graph is constructed. Then the corrected reads are mapped on the first de Bruijn graph to generate SR that are used to construct a second de Bruijn graph. The corrected reads are once again mapped on the second de Bruijn graph to produce new SR and so on.</figcaption></figure>

First, the reads are corrected to contain as low sequencing error as possible to enable fast mapping of the reads on the graph. Secondly, the unitigs are computed and a tip removal step is performed to further clean the graph. Thirdly, reads are mapped on the de Bruijn graph, generating SR. In a fourth step, the redundant SR are removed from the SR set. Two phenomenons are responsible for redundancy. First, a very high number of reads may map on large unitigs. Many super-reads will therefore just consist in one large unitig. Secondly, some super-reads may be included in other super-reads and are therefore not informative. Such super-reads are discarded and only Maximal Super Reads (SR that are not included in any other) are conserved. Then MSR sequences are used to build a higher order de Bruijn graph.

We show the interest of our iterative approach in order to construct a very high k de Bruijn graph. Let us fix the goal of constructing a high order de Bruijn graph from 250 base pairs reads. We start by building a de Bruijn graph with $k = 51$ from corrected reads, then we use the MSR to obtain kmers of superior order. We create the successive graphs by adding 50 to the order at each step, results about unitig numbers, false positives and false negatives are presented in Table 4.2. We can see that we are able to construct a very high order de Bruijn graph with very few holes and errors, showing the interest of this approach.

### 4.2.2 Beyond read length

In order to further improve our assembler, we integrated the use of paired end reads. PE reads come from very close regions of the genome. Classical PE libraries have a fragment size between 300 and 800bp. In some cases the paired reads can even overlap if the size of the two reads is larger than the fragment size (for example pair of 250 bp reads with an fragment size of 450). The idea to use this longer distance information is to map a pair of reads separately on the graph and to create a single SR if the two SR generated by the read pair overlap. The toy example in Figure 4.9 shows that we can create a long
<table><caption>Table 4.2: High order de Bruijn graph construction for 100X of 250 bp reads with 1% error rate from E.coli containing more than 4.6 million 251-mer.</caption><thead><tr><th>Order</th><th>Unitig number</th><th>FP</th><th>FN</th></tr></thead><tbody><tr><td>k=51</td><td>949</td><td>408</td><td>32</td></tr><tr><td>k=101</td><td>447</td><td>72</td><td>32</td></tr><tr><td>k=151</td><td>329</td><td>25</td><td>32</td></tr><tr><td>k=201</td><td>270</td><td>54</td><td>32</td></tr><tr><td>k=251</td><td>222</td><td>198</td><td>92</td></tr></tbody></table>

<figure><img src="image_70.png" /><figcaption>Figure 4.9: How to construct a high order de Bruijn graph with read mapping. A first de Bruijn graph is constructed with $k = 5$. Reads are mapped on it, producing super-reads. A new de Bruijn graph is constructed from the super-reads with a higher order, here $k = 7$. With paired reads we are able to phase more distant events or to solve larger repeats. In this example the purple reads in the right figure are paired. They produce longer SR than the unpaired reads from the left figure and are able to phase the three variants.</figcaption></figure>

SR from a pair of non overlapping reads. The simpler the graph is, the longer the unitigs are and the more we will be able to connect distant paired reads. A SR generated from a pair of reads can be larger than the fragment size of the pair. Using this fact, we can generate SR longer than the reads and construct a de Bruijn graph with an order above the read length.

We already shown in the previous table a de Bruijn graph with $k = 251$ constructed from reads of length 250. However paired end reads allow both an easier process and the use of orders that can be superior to the size of the reads (Table 4.3). Following this principle, there is actually no upper bound on the distance we may authorize between reads as soon as we have unitigs large enough to connect them. Further work will consist in including longer range information as mate pair reads that present a way larger fragment size. We can also consider adding different types of data such as long reads but such data would require specific steps in the workflow. The use of long range information should solve large repeats and phase distant variants, providing a more continuous assembly, haploid or polyploid.
<table><caption>Table 4.3: High order de Bruijn graph construction for 100X of 250 bp reads from E.Coli containing more than 4.6 million 301-mer.</caption><thead><tr><th>Order</th><th>Unitig number</th><th>FP</th><th>FN</th></tr></thead><tbody><tr><td>k=51</td><td>949</td><td>408</td><td>32</td></tr><tr><td>k=101</td><td>576</td><td>1,208</td><td>32</td></tr><tr><td>k=151</td><td>356</td><td>341</td><td>32</td></tr><tr><td>k=201</td><td>276</td><td>125</td><td>32</td></tr><tr><td>k=251</td><td>224</td><td>210</td><td>108</td></tr><tr><td>k=301</td><td>122</td><td>67</td><td>5,992</td></tr></tbody></table>

### "Building very high order de Bruijn graph" core messages:

* Super Reads created from read mapping can be seen as trusted sequences

* Iterative mapping allows the construction of high order de Bruijn graph

* Paired end reads allow the construction of a de Bruijn graph with *k* above the reads length

## 4.3 Assembly results

### 4.3.1 Haploid genome

In Table 4.4 we present some assembly results on various haploid genomes to assess the efficiency of our proposed solution against other assemblers.

We observe that BWISE is able to produce contigs as long as state of the art assemblers. We also note the scalability of BWISE that successfully performed a human genome assembly with less than 60 GB of RAM in two days. SPAdes and Platanus were

<table><caption>Table 4.4: Simulated of 100x datasets on various genomes. Simulated reads are paired with 250 base pairs and a fragment size of 800.</caption><thead><tr><th>Assembly</th><th>#contigs</th><th>N50</th><th>N80</th></tr></thead><tbody><tr><td><i>E. coli</i> BWISE</td><td>56</td><td>210,994</td><td>126,805</td></tr><tr><td><i>E. coli</i> SPAdes</td><td>71</td><td>178,400</td><td>83,005</td></tr><tr><td><i>E. coli</i> Platanus</td><td>272</td><td>133,252</td><td>28,361</td></tr><tr><td><i>E. coli</i> Minia</td><td>309</td><td>107,878</td><td>40,937</td></tr><tr><td><i>C. elegans</i> BWISE</td><td>2,527</td><td>122,287</td><td>56,547</td></tr><tr><td><i>C. elegans</i> SPAdes</td><td>6,550</td><td>103,128</td><td>43,904</td></tr><tr><td><i>C. elegans</i> Platanus</td><td>21,413</td><td>52,661</td><td>18,330</td></tr><tr><td><i>C. elegans</i> Minia</td><td>35,317</td><td>23,641</td><td>6,906</td></tr><tr><td>Human BWISE</td><td>132,272</td><td>74,148</td><td>29,705</td></tr><tr><td>Human Minia</td><td>884,788</td><td>20,202</td><td>7,381</td></tr></tbody></table>
<table><caption>Table 4.5: Assembly comparison of simulated 100x datasets on *E. coli* genome with simulated heterozygosity. Simulated reads are paired with 250 base pairs and a fragment size of 800.</caption><thead><tr><th>Assembly</th><th>#contigs</th><th>N50</th><th>N80</th></tr></thead><tbody><tr><td><i>E. coli</i> 1% BWISE</td><td>307</td><td>64,217</td><td>29,717</td></tr><tr><td><i>E. coli</i> 1% SPAdes</td><td>5,027</td><td>1,827</td><td>728</td></tr><tr><td><i>E. coli</i> 1% Platanus</td><td>2,887</td><td>6,526</td><td>1,999</td></tr><tr><td><i>E. coli</i> 2% BWISE</td><td>60</td><td>325,047</td><td>181,988</td></tr><tr><td><i>E. coli</i> 2% SPAdes</td><td>2,411</td><td>12,632</td><td>3,450</td></tr><tr><td><i>E. coli</i> 2% Platanus</td><td>2,488</td><td>6,896</td><td>2,028</td></tr><tr><td><i>E. coli</i> 3% BWISE</td><td>8</td><td>2,610,736</td><td>605,526</td></tr><tr><td><i>E. coli</i> 3% SPAdes</td><td>121</td><td>246,674</td><td>116,251</td></tr><tr><td><i>E. coli</i> 3% Platanus</td><td>1,614</td><td>10,019</td><td>3,157</td></tr><tr><td><i>E. coli</i> 4% BWISE</td><td>5</td><td>3,424,913</td><td>2,655,158</td></tr><tr><td><i>E. coli</i> 4% SPAdes</td><td>4</td><td>3,424,851</td><td>3,424,851</td></tr><tr><td><i>E. coli</i> 4% Platanus</td><td>1,071</td><td>13,346</td><td>5,279</td></tr><tr><td><i>E. coli</i> 5% BWISE</td><td>2</td><td>4,639,590</td><td>4,639,574</td></tr><tr><td><i>E. coli</i> 5% SPAdes</td><td>3</td><td>4,639,654</td><td>4,639,654</td></tr><tr><td><i>E. coli</i> 5% Platanus</td><td>711</td><td>22,275</td><td>9,117</td></tr></tbody></table>

not able to perform the human assembly on the used machine with 250 GB of RAM.

### 4.3.2 Diploid genome

As seen before, assembly of highly heterozygous genomes is a case where the number of small repeats can be tremendous and where the use of the whole read length can lead to very significant improvements. In Table 4.5 we present assembly results of several tools, including BWISE, when used on regions with a known uniform heterozygosity. We first simulated uniform heterozygosity on haploid genomes and thereafter simulated reads from the artificial diploid references. We may argue that this does not represent a realistic scenario as the simulated heterozygosity "destroys" the repeats along the genomes. If a large repeat appear in the genome, we may add different variations on the different occurrences of the repeats. Thus, each occurrence of the repeat may be different from the other in the artificial diploid references. Therefore, presented results should only be seen as comparative evaluation. We can make several observations from these results. First we see that high heterozygosity is easy to solve and that tools achieve to produce long phased contigs since most variants were close to each other. Still, it is interesting to see that BWISE performs better than SPAdes on all rates. Especially on medium rates where SPAdes is not able to produce long contigs, results show that BWISE is able to do so.

On wider genomes, Table 4.6 shows that BWISE present similar performances and is still able to produce large contigs. On such dataset SPAdes and Platanus were not able to complete on the used server (2 days of timeout and 250 GB of RAM).
<table><caption>Table 4.6: Assembly comparison of simulated 100x datasets on *C. elegans* genome with simulated heterozygosity. Simulated reads are paired with 250 base pairs and a fragment size of 800. SPAdes exceeded the timeout.</caption><thead><tr><th>Assembly</th><th>#contigs</th><th>N50</th><th>N80</th></tr></thead><tbody><tr><td><i>C. elegans</i> 1% BWISE</td><td>9,505</td><td>50,870</td><td>22,850</td></tr><tr><td><i>C. elegans</i> 2% BWISE</td><td>1,674</td><td>310,249</td><td>136,205</td></tr><tr><td><i>C. elegans</i> 3% BWISE</td><td>429</td><td>1,089,405</td><td>488,453</td></tr><tr><td><i>C. elegans</i> 4% BWISE</td><td>196</td><td>2,626,274</td><td>1,074,422</td></tr></tbody></table>

<table><caption>Table 4.7: Assembly comparison of simulated 100x datasets on *E. coli* genome with low simulated heterozygosity. Simulated reads are paired with 250 base pairs and a fragment size of 800.</caption><thead><tr><th>Assembly</th><th>#contigs</th><th>N50</th><th>N80</th></tr></thead><tbody><tr><td><i>E. coli</i> 0.5% BWISE</td><td>1,370</td><td>10,390</td><td>5,633</td></tr><tr><td><i>E. coli</i> 0.5% SPAdes</td><td>973</td><td>17,487</td><td>6,778</td></tr><tr><td><i>E. coli</i> 0.5% Platanus</td><td>2,224</td><td>11,128</td><td>3,671</td></tr><tr><td><i>E. coli</i> 0.1% BWISE</td><td>3,835</td><td>4,520</td><td>3,233</td></tr><tr><td><i>E. coli</i> 0.1% SPAdes</td><td>55</td><td>210,905</td><td>88,598</td></tr><tr><td><i>E. coli</i> 0.1% Platanus</td><td>926</td><td>40,837</td><td>14,688</td></tr></tbody></table>

### 4.3.3 Future works

The essential limitation of BWISE workflow is noticeable when the heterozygosity rate is low (Table 4.7). Each time variants are pairwise too far away from each other, the assembly is broken since the contigs cannot be extended in a sure way. The results show that with 0.5% heterozygosity, neither BWISE or SPAdes are able to produce long contigs. On 0.1% SPAdes is able to produce long contigs where BWISE produces very short ones. This can be explained as SPAdes crushes bubbles and therefore produces a haploid assembly of 4.6 mega-bases while BWISE tries to produce a diploid assembly of 9.2 mega-bases.

The fact that we are not able to phase all heterozygous regions raise the question of the assembly representation. As it can be seen on Figure 4.10, the first and last variants were not phased by the reads, and BWISE produces the four contigs of solution 1. While being correct, this format may not be optimal for each usage. In some cases, we may be interested to get longer sequences by crushing bubbles. Another way to solve this problem would be to add longer range information in order to phase more distant variants. With more distant linked reads (as mate pair reads) phasing more distant variant should be possible and therefore produce more contiguous assembly.

The proposed solution based on successive mapping provides very interesting result and scale up to large genome. However several aspects have to be studied as the assembly quality, the behavior of BWISE with different coverage and the integration of longer range information. Several operations of the proposed workflow are quite redundant, as the repeated mapping of the whole read set. Thus, optimization of the workflow could
<figure><img src="image_71.png" /><figcaption>Figure 4.10: Potential representations of non phases contigs. The first solution proposed is a bit redundant. To limit the redundancy we can crush the remaining small bubbles as showed in solution 2. To further improve contigs length we could also extend the contigs with the merged bubbles sequences, as it is done in solution 3.</figcaption></figure>
lead to sensible performances improvement. This method is not completely explored
and in depth investigations are required to provide a robust and usable assembler for
haploid and diploid genomes and eventually more complex assembly cases as pool-seq or
meta-genomes.

"Assembly results" core messages:

* Iterative mapping leads to state of the art haploid assembly

* Iterative mapping outperforms state of the art tools in highly heterozygous zones

* Low heterozygosity polyploid assembly is difficult with short reads alone
# Chapter 5

## Conclusions and perspectives
## 5.1 Conclusion

Two main directions were explored during this thesis. The first aim was to propose new structures to handle the scaling of de Bruijn graph usages. The second aim was to propose new methods and techniques to no longer consider the de Bruijn graph only as an assembler internal representation but as a proper reference structure that can be used for various applications. This part reminds and summarizes our contributions and main results toward those objectives.

Several reasons justify the need for the development of efficient data structures in genome assembly. The assembly of human sized genomes still presents important scalability issues, with the need of very important resources. Very large genomes are common across the living organisms, especially plant genomes, and efficient assembly methods for such genomes are necessary to remove the present locks. The second interest of efficient methods, for any genome, is to democratize the treatments of such datasets and make their analysis less costly and more accessible. In order to address the large genome assembly scalability using de Bruijn graphs, we proposed several practical and theoretical results. We described a new category of data structure called Navigational Data Structure that enable de Bruijn graph assembly with lower space requirements than the previously known bounds. Implementations of such structures could lead to extremely memory efficient assemblers. We proposed a proof of concept low memory usage assembler called DGBFM. The core idea of this method is to index the simple paths (unitigs) of the de Bruijn graph instead of its kmers into a memory efficient data structure, here a FM-index. DGBFM present of very low memory usage and is able to assemble a human genome with 1.5 GB of RAM. A take-home message is that considering a de Bruijn graph as a set of unitigs is an efficient practical choice. However the construction of such sets in low memory is not trivial. We therefore proposed a new method based on external memory to construct this set called BCALM. This method was slow and sequential, and the construction of the compacted de Bruijn graph was the bottleneck of DGBFM as in many assemblers. We therefore conceived a efficient and parallel method called BCALM2 to compute de Bruijn graph unitigs. BCALM2, relying on external memory present order of magnitude less memory that known implementation and was also significantly faster. Such a method can lead to a low memory assembly even for genome larger than the human one. In BCALM2, and in several other applications, we encountered the challenge of indexing extremely large numbers of object. The need for memory efficient and fast structure lead us to use minimal perfect hash functions (MPHF). Existing MPHF library presented very high memory usage and running time during construction. Thus we proposed a new scalable MPHF library, designed to construct MPHF from very large key set with moderate resources called BBHASH. The main property of BBHASH is the ability to construct a MPHF without memory overhead, and the construction is also parallel and faster than other available implementations. This implementation could be used virtually in any HPC field. We believe that such an efficient data structure could especially benefit in numerous applications in bioinformatics, where the need to index large number of objects is fundamental.

The de Bruijn graph is mainly seen as an assembly graph, a simple temporary structure used to generate contigs. If the de Bruijn graph is mainly used for assembly, many properties of this structure have been used outside of its first domain of application. For
numerous reasons, the classical representation of a genome into linear sequences is not fit-
ted for all application. To cope with these problems, several new genome representations
based on graph are proposed. If the de Bruijn graph is not an efficient representation for
finished genomes, its properties can be used in advantageous ways in several applications.
With efficient techniques to construct a (compacted) de Bruijn graph, we argued that it
could be a good reference to represent one or multiple genomes. We presented different
examples of the interest of such a representation. The unfinished genomes, where the de
Bruijn graph may contain more information than the contigs sequences. Heterozygous
genomes, as a de Bruijn graph contain and can represent efficiently the different varia-
tions that may be present in such genomes. Very redundant dataset that can be factored
by the efficient representation of the repeated sequences in a de Bruijn graph.

Multiple applications could use such a reference graph directly, as read set comparison
in meta-genomics or mapping read on the de Bruijn graph sequences. Thus, we argued
that we need a read mapping method working directly on the de Bruijn graph structure
to permit its versatile use as a reference. In order to enable the use of such representation,
we provided a tool able to map reads efficiently on a de Bruijn graph called BGREAT.
Despite the fact that we proved that mapping reads on a de Bruijn graph is a NP-
Complete problem, the use of efficient heuristics allows BGREAT to be very efficient
on current datasets. We presented several direct advantages of such a method. Results
confirm that mapping reads directly on the de Bruijn graph instead of mapping on contigs
allows a higher ratio of read mapped, showing the interest of mapping on the de Bruijn
graph over mapping on its sequences. We also proved that the de Bruijn graph is a
very accurate reference with very few sequencing errors and keeping the heterozygosity
information. This facts show the de Bruijn graph as a very sound and easy to construct
structure. Using this property, we presented an out of the box usage with the read
corrector BCOOL. Results of the proof of concept read correction by mapping them on
the de Bruijn graph are very interesting. BCOOL corrected several time more errors than
state of the art correctors on high coverage simulated data. Outside of read correction, we
believe that such a method could be used in very diverse usage of NGS data as assembly,
scaffolding, read compression or variant calling.

Applied to genome assembly, the de Bruijn graph is essentially used for its scalability,
but the splitting of the reads into kmers fragments the assembly. One of the motiva-
tion behind the conception of a tool like BGREAT was to improve the de Bruijn graph
assembly continuity. The idea was to map reads on a de Bruijn graph in order to add
back the reads information into de Bruijn graph assembly. To use this information, we
proposed a new method to build a high order de Bruijn graph using the read length at a
reasonable cost. The idea is to map reads on the de Bruijn graph and to use the mapped
path sequences to build a de Bruijn graph with a higher order. We implemented this
idea into a proof of concept assembler called BWISE. Similarly to the string graphs, the
produced high order de Bruijn graph are able to use almost the whole length of the reads
in order to solve repeats and produce more contiguous assembly than regular de Bruijn
graph. We also showed that such an approach could benefit of paired end reads in order
to produce de Bruijn graph with kmer size above the read length. Using those paired
reads information improve the assembly continuity and the size of the contigs without
a scaffolding step as we work directly on the de Bruijn graph. Preliminary results on
simulated data on haploid genomes exhibit assemblies equivalent to state of the art as-
semblers. However, the most interesting results of this approach was a good assembly
continuity in heterozygous regions were state of the art assemblers produce assembly of
very poor quality. On our simulated data, BWISE is able to produce phased contigs and
to assemble separately the different allele where other de Bruijn graph assembler are not
able to make the distinction. We believe that the next assemblers will have to propose
such heterozygosity awareness and that BWISE is a promising approach to do so.

Those results convey the idea that the de Bruijn graph as a data structure has not
been used yet to the maximum of its potential.

## 5.2 Proposed methods

This thesis resulted in the following open-source softwares:

**BCALM** Proof of concept of de Bruijn graph compaction in low memory. BCALM was used in the DBGFM assembler to compute the unitigs to index and is described in [74].
Available https://github.com/Malfoy/BCALM

**BCALM2** Efficient and parallel de Bruijn graph compaction in low memory. BCALM2 was used for the de Bruijn graph results presented in this document. It is used for the de Bruijn graph compaction before the use of BGREAT2 in BCOOL and BWISE. It is described in [82]. Available https://github.com/GATB/bcalm

**BBHASH** Scalable minimal perfect hash function construction. BBHASH MPHF is used in BCALM2, GATB and SRC softwares. It is described in [86]. Available https://github.com/rizkg/BBHash

**BGREAT** Proof of concept of read mapping on de Bruijn graph. It is described in [116].
Available https://github.com/Malfoy/BGREAT

**BGREAT2** Read mapping on de Bruijn graph. An improved and self contained version of BGREAT. As said before it is used with BCALM2 in BCOOL and BWISE workflows.
The improvements applied are described in chapter 3. Available https://github.com/Malfoy/BGREAT2

**BCOOL** de Bruijn graph based read correction. This read corrector and its preliminary results are presented in chapter 3. Available https://github.com/Malfoy/BCOOL

**BWISE** High order de Bruijn graph Assembly. This assembler and its preliminary results on haploid and diploid genomes are presented in chapter 4. Available https://github.com/Malfoy/BWISE
## 5.3 Future of sequencing

In 2008, the next generation sequencing revolutionized the genomic field. While it is hard to forecast the future of such fast evolving and diverse technologies, we can observe emerging usages of new kinds of data and imagine the potential benefits and challenges they may bring.

### 5.3.1 Third generation sequencing

New sequencing technologies are arising. Called third generation sequencing, those technologies propose order of magnitude longer sequences, but currently with a very high error rate (up to 15% [127]). Both Pacific Bioscience [6] and Oxford Nanopore [7] technologies are improving and may provide even longer reads with a low error rate in the near future. Those technologies have the potential to be a major shift in genome assembly. With such long sequences, most repeats could be solved and the "one contig one chromosome" dream could be achieved. But those sequences still present huge challenges. Tools have to deal with the important error rate. To make it worse, the type of error is mainly insertions and deletions errors which are harder to deal with than the substitutions usually encountered in NGS reads. The costs of such sequencing are currently high compared to NGS, while a high coverage is necessary in order to cope with the error rate. To assemble such sequences, the de Bruijn graph is not usable in a straightforward way. The string graph is therefore recovering its central place in assembly. But the computational cost of such an assembly is going through the roof and many tools have been proposed to make the overlap detection phase extremely efficient. Two approaches are advocated. The hybrid assembly approaches are based on the assumption that long reads could be used coupled with short reads to produce a high quality assembly in a cost efficient way. Either by correcting long reads with short reads [128] [129], or to scaffold short reads contigs with long reads [130] or assemble long and short reads together directly [131]). The pure assembly approaches [132] argues that long reads alone have the potential to produce the best assemblies because of their less biased sequencing protocols and that we should work on the computational cost of the assemblers.

### 5.3.2 Long range short reads sequencing

New short reads based technologies are also making an entrance into bioinformatics. The 3C protocol [133] could provide read pairs from the same chromosome with no upper bound on the insert size and therefore help phasing and solving repeats with no upper bound of size [51]. The 10X [134] technology provides short reads localized on short genome region of size around 10,000 bases. This kind of sequencing has shown great potential for phasing since reads from a pack come from the same allele [135] while it seems a priori harder to do so with long reads given their error rates. But for haploid assembly those technologies also great potential for solving repeats while being accurate sequences. The long reads present the advantages of representing a continuous sequences while 10X, like paired reads give sequences separated by holes.
## 5.4 Future of this work and perspectives

The different methods proposed in this documents could be continued in many ways, with refinements, new applications or integrations of other kinds of information. Some perspective considered are presented here.

While assembly of short reads is often considered as a well known subject, large genomes still require tremendous resources. But even for smaller genomes, faster structures and assembly methods will be required for the throughput of assemblers to follow the sequencer one. An efficient assembler, implementing a NDS based on unitigs indexed on a MPHF, could provide an extremely efficient method to work on very large genomes without current scalability problems or provide assembly with a very high throughput.

We showed that, beyond assembly, the de Bruijn graph is an extremely efficient data structure to represent a genome, or a set of genomes, with interesting properties that could be used in other fields than assembly. A scalable read corrector providing a very high ratio of perfect read would be a very useful tool for applications as assembly or variant calling. Our proposed corrector BCOOL could be tuned to adapt to real data and to different coverage in order to provide almost perfect short reads.

Several other applications of mapping read on a de Bruijn graph can be found as quantification, variant calling or read compression could be proposed in the near future. As the read correction, reads could be mapped on the graph and we could encode their path in the de Bruijn graph instead of their sequences in order to compress them. More direct usages would be quantification in genomic, meta-genomic or RNA-seq by mapping the read on a de Bruijn graph and extracting information from the amount of reads mapped on given components of the de Bruijn graph. Simple variant calling have been shown to be done directly from a de Bruijn graph structure, an interesting perspective would be to use the reads mapping on this graph and use this information to capture more complex cases.

Polyploid assembly is also a widely encountered problem while representing a very small part of the literature. We pointed those problems and presented new solutions and perspectives to address those challenges. BWISE polyploid assembly methods have still to be tuned to fit to the challenges of the real data but the door to the assembly of polyploid genomes is open. The use of paired reads directly on a de Bruijn graph could also inspire new scaffolding techniques that are mainly based on read mapping on contigs. The use of other kind of information in order to produce larger contigs is also a natural prolongation of this work. 10X reads seems a really promising technique to improve the polyploid assembly of BWISE, as the mapping method does not need important changes. The length of 10X DNA fragments could lead to extremely large and phased contigs.

We essentially talked about polyploid assembly but techniques used by BWISE could also be relevant for meta-genomic assembly. The read mapping could help distinguish the different individuals present in a de Bruijn graph as it is used now to separate the different alleles. To do so we have to address the challenges of heterogeneous coverages that can exist in meta-genomic. Almost the same argumentation could be used for transcriptome assembly of RNA sequences, an efficient assembler able to use Paired End data could tackle the combinatorial aspect of RNA assembly.

As BGREAT is able to align short reads on a de Bruijn graph, we are naturally interested to adapt our methods to be able to align long reads on a de Bruijn graph.
This kind of tool could have, as BGREAT, several applications. The alignment of the long reads on the de Bruijn graph could be used to correct their errors, or the alignment information could be used to order contigs or unitigs of the de Bruijn graph. Such a method could also be used in the context of pure long reads assembly as a efficient mean to detect overlap between them. The computational cost of long reads assembly has also to be addressed and efficient data structures will be required as well in order to efficiently assemble the future deluge of such datasets. We believe that with decreasing error rate the de Bruijn graph could be useful even for pure long reads assembly [136] as the de Bruijn graph paths represent easy to compute efficient consensus sequences.
Bibliography

[1] James D Watson and Francis HC Crick. A structure for deoxyribose nucleic acid. *Nature*, 171, 2004.

[2] Eric S Lander, Lauren M Linton, Bruce Birren, Chad Nusbaum, Michael C Zody, Jennifer Baldwin, Keri Devon, Ken Dewar, Michael Doyle, William FitzHugh, et al. Initial sequencing and analysis of the human genome. *Nature*, 409(6822):860–921, 2001.

[3] J Craig Venter, Mark D Adams, Eugene W Myers, Peter W Li, Richard J Mural, Granger G Sutton, Hamilton O Smith, Mark Yandell, Cheryl A Evans, Robert A Holt, et al. The sequence of the human genome. *science*, 291(5507):1304–1351, 2001.

[4] Frederick Sanger, Steven Nicklen, and Alan R Coulson. Dna sequencing with chain-terminating inhibitors. *Proceedings of the National Academy of Sciences*, 74(12):5463–5467, 1977.

[5] David R Bentley, Shankar Balasubramanian, Harold P Swerdlow, Geoffrey P Smith, John Milton, Clive G Brown, Kevin P Hall, Dirk J Evers, Colin L Barnes, Helen R Bignell, et al. Accurate whole human genome sequencing using reversible terminator chemistry. *nature*, 456(7218):53–59, 2008.

[6] John Eid, Adrian Fehr, Jeremy Gray, Khai Luong, John Lyle, Geoff Otto, Paul Peluso, David Rank, Primo Baybayan, Brad Bettman, et al. Real-time dna sequencing from single polymerase molecules. *Science*, 323(5910):133–138, 2009.

[7] Andrew H Laszlo, Ian M Derrington, Brian C Ross, Henry Brinkerhoff, Andrew Adey, Ian C Nova, Jonathan M Craig, Kyle W Langford, Jenny Mae Samson, Riza Daza, et al. Decoding long nanopore sequencing reads of natural dna. *Nature biotechnology*, 32(8):829–833, 2014.

[8] Daniel Aird, Michael G Ross, Wei-Sheng Chen, Maxwell Danielsson, Timothy Fen-
nell, Carsten Russ, David B Jaffe, Chad Nusbaum, and Andreas G nirke. Analyzing
and minimizing pcr amplification bias in illumina sequencing libraries. *Genome
biology*, 12(2):R18, 2011.

[9] Juliane C Dohm, Claudio Lottaz, Tatiana Borodina, and Heinz Himmelbauer. Substantial biases in ultra-short read data sets from high-throughput dna sequencing. *Nucleic acids research*, 36(16):e105–e105, 2008.
[10] Project Encode. The encode (encyclopedia of dna elements) project. *Science*, 306(5696):636-640, 2004.

[11] Sushmita Roy, Jason Ernst, Peter V Kharchenko, Pouya Kheradpour, Nicolas Negre, Matthew L Eaton, Jane M Landolin, Christopher A Bristow, Lijia Ma, Michael F Lin, et al. Identification of functional elements and regulatory circuits by drosophila modencode. *Science*, 330(6012):1787-1797, 2010.

[12] Jennifer Harrow, Adam Frankish, Jose M Gonzalez, Electra Tapanari, Mark Diekhans, Felix Kokocinski, Bronwen L Aken, Daniel Barrell, Amonida Zadissa, Stephen Searle, et al. Gencode: the reference human genome annotation for the encode project. *Genome research*, 22(9):1760-1774, 2012.

[13] Todd J Treangen and Steven L Salzberg. Repetitive dna and next-generation sequencing: computational challenges and solutions. *Nature Reviews Genetics*, 13(1):36-46, 2012.

[14] Jonghwan Kim, Akshay A Bhinge, Xochitl C Morgan, and Vishwanath R Iyer. Mapping dna-protein interactions in large genomes by sequence tag analysis of genomic enrichment. *Nature Methods*, 2(1):47-53, 2005.

[15] David S Johnson, Ali Mortazavi, Richard M Myers, and Barbara Wold. Genome-wide mapping of in vivo protein-dna interactions. *Science*, 316(5830):1497-1502, 2007.

[16] Kari-Jouko Räihä and Esko Ukkonen. The shortest common supersequence problem over binary alphabet is np-complete. *Theoretical Computer Science*, 16(2):187-198, 1981.

[17] Granger G Sutton, Owen White, Mark D Adams, and Anthony R Kerlavage. Tigr assembler: A new tool for assembling large shotgun sequencing projects. *Genome Science and Technology*, 1(1):9-19, 1995.

[18] Xiaoqiu Huang and Anup Madan. Cap3: A dna sequence assembly program. *Genome research*, 9(9):868-877, 1999.

[19] René L Warren, Granger G Sutton, Steven JM Jones, and Robert A Holt. Assembling millions of short dna sequences using ssake. *Bioinformatics*, 23(4):500-501, 2007.

[20] William R Jeck, Josephine A Reinhardt, David A Baltrus, Matthew T Hickenbotham, Vincent Magrini, Elaine R Mardis, Jeffery L Dangl, and Corbin D Jones. Extending assembly of short dna sequences to handle error. *Bioinformatics*, 23(21):2942-2944, 2007.

[21] Eugene W Myers, Granger G Sutton, Art L Delcher, Ian M Dew, Dan P Fasulo, Michael J Flanigan, Saul A Kravitz, Clark M Mobarry, Knut HJ Reinert, Karin A Remington, et al. A whole-genome assembly of drosophila. *Science*, 287(5461):2196-2204, 2000.
[22] Stephen F Altschul, Warren Gish, Webb Miller, Eugene W Myers, and David J Lipman. Basic local alignment search tool. *Journal of molecular biology*, 215(3):403-410, 1990.

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

[24] Alexey Gurevich, Vladislav Saveliev, Nikolay Vyahhi, and Glenn Tesler. Quast: quality assessment tool for genome assemblies. *Bioinformatics*, 29(8):1072-1075, 2013.

[25] Eugene W Myers. The fragment assembly string graph. *Bioinformatics*, 21(suppl 2):ii79-ii85, 2005.

[26] Jared T Simpson and Richard Durbin. Efficient construction of an assembly string graph using the fm-index. *Bioinformatics*, 26(12):i367-i373, 2010.

[27] Jared T Simpson and Richard Durbin. Efficient de novo assembly of large genomes using compressed data structures. *Genome research*, 22(3):549-556, 2012.

[28] Nicolaas Govert De Bruijn. A combinatorial problem. 1946.

[29] Dmitri Loguinov, Anuj Kumar, Vivek Rai, and Sai Ganesh. Graph-theoretic analysis of structured peer-to-peer systems: routing distances and fault resilience. In *Proceedings of the 2003 conference on Applications, technologies, architectures, and protocols for computer communications*, pages 395-406. ACM, 2003.

[30] M Frans Kaashoek and David R Karger. Koorde: A simple degree-optimal distributed hash table. In *International Workshop on Peer-to-Peer Systems*, pages 98-107. Springer, 2003.

[31] Pavel A Pevzner, Haixu Tang, and Michael S Waterman. An eulerian path approach to dna fragment assembly. *Proceedings of the National Academy of Sciences*, 98(17):9748-9753, 2001.

[32] Paul Medvedev, Konstantinos Georgiou, Gene Myers, and Michael Brudno. Computability of models for sequence assembly. In *International Workshop on Algorithms in Bioinformatics*, pages 289-301. Springer Berlin Heidelberg, 2007.

[33] Daniel R Zerbino and Ewan Birney. Velvet: algorithms for de novo short read assembly using de bruijn graphs. *Genome research*, 18(5):821-829, 2008.

[34] Rayan Chikhi and Paul Medvedev. Informed and automated k-mer size selection for genome assembly. *Bioinformatics*, page btt310, 2013.

[35] Ruibang Luo, Binghang Liu, Yinlong Xie, Zhenyu Li, Weihua Huang, Jianying Yuan, Guangzhu He, Yanxiang Chen, Qi Pan, Yunjie Liu, et al. Soapdenovo2: an empirically improved memory-efficient short-read de novo assembler. *Gigascience*, 1(1):18, 2012.
[36] Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A Gurevich, Mikhail Dvorkin, Alexander S Kulikov, Valery M Lesin, Sergey I Nikolenko, Son Pham, Andrey D Prjibelski, et al. Spades: a new genome assembly algorithm and its applications to single-cell sequencing. *Journal of Computational Biology*, 19(5):455-477, 2012.

[37] Mark J Chaisson and Pavel A Pevzner. Short read fragment assembly of bacterial genomes. *Genome research*, 18(2):324-330, 2008.

[38] Jonathan Butler, Iain MacCallum, Michael Kleber, Ilya A Shlyakhter, Matthew K Belmonte, Eric S Lander, Chad Nusbaum, and David B Jaffe. Allpaths: de novo assembly of whole-genome shotgun microreads. *Genome research*, 18(5):810-820, 2008.

[39] Yu Peng, Henry CM Leung, Siu-Ming Yiu, and Francis YL Chin. Idba—a practical iterative de bruijn graph de novo assembler. In *Annual International Conference on Research in Computational Molecular Biology*, pages 426–440. Springer Berlin Heidelberg, 2010.

[40] Jared T Simpson, Kim Wong, Shaun D Jackman, Jacqueline E Schein, Steven JM Jones, and Inanc Birol. Abyss: a parallel assembler for short read sequence data. *Genome research*, 19(6):1117-1123, 2009.

[41] Ruiqiang Li, Hongmei Zhu, Jue Ruan, Wubin Qian, Xiaodong Fang, Zhongbin Shi, Yingrui Li, Shengting Li, Gao Shan, Karsten Kristiansen, et al. De novo assembly of human genomes with massively parallel short read sequencing. *Genome research*, 20(2):265-272, 2010.

[42] Sante Gnerre, Iain MacCallum, Dariusz Przybylski, Filipe J Ribeiro, Joshua N Burton, Bruce J Walker, Ted Sharpe, Giles Hall, Terrance P Shea, Sean Sykes, et al. High-quality draft assemblies of mammalian genomes from massively parallel sequence data. *Proceedings of the National Academy of Sciences*, 108(4):1513-1518, 2011.

[43] Thomas Conway, Jeremy Wazny, Andrew Bromage, Justin Zobel, and Bryan Beresford-Smith. Gossamer—a resource-efficient de novo assembler. *Bioinformatics*, 28(14):1937-1938, 2012.

[44] Rayan Chikhi and Guillaume Rizk. Space-efficient and exact de bruijn graph representation based on a bloom filter. In *International Workshop on Algorithms in Bioinformatics*, pages 236–248. Springer, 2012.

[45] Niranjan Nagarajan and Mihai Pop. Parametric complexity of sequence assembly: theory and applications to next generation sequencing. *Journal of computational biology*, 16(7):897-908, 2009.

[46] Mark Chaisson, Pavel Pevzner, and Haixu Tang. Fragment assembly with short reads. *Bioinformatics*, 20(13):2067-2074, 2004.
[47] Daniel H Huson, Knut Reinert, and Eugene W Myers. The greedy path-merging algorithm for contig scaffolding. *Journal of the ACM (JACM)*, 49(5):603–615, 2002.

[48] Marten Boetzer, Christiaan V Henkel, Hans J Jansen, Derek Butler, and Walter Pirovano. Scaffolding pre-assembled contigs using sspace. *Bioinformatics*, 27(4):578–579, 2011.

[49] Marten Boetzer and Walter Pirovano. Toward almost closed genomes with gapfiller. *Genome biology*, 13(6):1, 2012.

[50] Song Gao, Denis Bertrand, Burton KH Chia, and Niranjan Nagarajan. Opera-lg: Efficient and exact scaffolding of large, repeat-rich eukaryotic genomes with performance guarantees. *Genome biology*, 17(1):1, 2016.

[51] Jean-François Flot, Hervé Marie-Nelly, and Romain Koszul. Contact genomics: scaffolding and phasing (meta) genomes using chromosome 3d physical signatures. *FEBS letters*, 589(20PartA):2966–2974, 2015.

[52] Rei Kajitani, Kouta Toshimoto, Hideki Noguchi, Atsushi Toyoda, Yoshitoshi Ogura, Miki Okuno, Mitsuru Yabana, Masayuki Harada, Eiji Nagayasu, Haruhiko Maruyama, et al. Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads. *Genome research*, 24(8):1384–1395, 2014.

[53] Yu Peng, Henry CM Leung, Siu-Ming Yiu, and Francis YL Chin. Meta-idba: a de novo assembler for metagenomic data. *Bioinformatics*, 27(13):i94–i101, 2011.

[54] Toshiaki Namiki, Tsuyoshi Hachiya, Hideaki Tanaka, and Yasubumi Sakakibara. Metavelvet: an extension of velvet assembler to de novo metagenome assembly from short sequence reads. *Nucleic acids research*, 40(20):e155–e155, 2012.

[55] Jaume Pellicer, Michael F Fay, and Ilia J Leitch. The largest eukaryotic genome of them all? *Botanical Journal of the Linnean Society*, 164(1):10–15, 2010.

[56] Temple F Smith and Michael S Waterman. Identification of common molecular subsequences. *Journal of molecular biology*, 147(1):195–197, 1981.

[57] Serafim Batzoglou, David B Jaffe, Ken Stanley, Jonathan Butler, Sante Gnerre, Evan Mauceli, Bonnie Berger, Jill P Mesirov, and Eric S Lander. Arachne: a whole-genome shotgun assembler. *Genome research*, 12(1):177–189, 2002.

[58] Paolo Ferragina and Giovanni Manzini. Opportunistic data structures with applications. In *Foundations of Computer Science, 2000. Proceedings. 41st Annual Symposium on*, pages 390–398. IEEE, 2000.

[59] Ben Langmead, Cole Trapnell, Mihai Pop, and Steven L Salzberg. Ultrafast and memory-efficient alignment of short dna sequences to the human genome. *Genome biology*, 10(3):1, 2009.

[60] Heng Li and Richard Durbin. Fast and accurate short read alignment with burrows-wheeler transform. *Bioinformatics*, 25(14):1754–1760, 2009.
[61] Inanc Birol, Anthony Raymond, Shaun D Jackman, Stephen Pleasance, Robin Coope, Greg A Taylor, Macaire Man Saint Yuen, Christopher I Keeling, Dana Brand, Benjamin P Vandervalk, et al. Assembling the 20 gb white spruce (picea glauca) genome from whole-genome shotgun sequencing data. *Bioinformatics*, page btt178, 2013.

[62] Shaun D Jackman, Benjamin P Vandervalk, Hamid Mohamadi, Justin Chu, Sarah Yeo, S Austin Hammond, Golnaz Jahesh, Hamza Khan, Lauren Coombe, Rene L Warren, et al. Abyss 2.0: Resource-efficient assembly of large genomes using a bloom filter. *Genome Research*, pages gr-214346, 2017.

[63] Thomas C Conway and Andrew J Bromage. Succinct data structures for assembling large genomes. *Bioinformatics*, 27(4):479-486, 2011.

[64] Daisuke Okanohara and Kunihiko Sadakane. Practical entropy-compressed rank/select dictionary. In *Proceedings of the Meeting on Algorithm Engineering & Experiments*, pages 60–70. Society for Industrial and Applied Mathematics, 2007.

[65] Guillaume Marçais and Carl Kingsford. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. *Bioinformatics*, 27(6):764–770, 2011.

[66] Chris Purcell and Tim Harris. Non-blocking hashtables with open addressing. In *International Symposium on Distributed Computing*, pages 108–121. Springer, 2005.

[67] Guillaume Rizk, Dominique Lavenier, and Rayan Chikhi. Dsk: k-mer counting with very low memory usage. *Bioinformatics*, page btt020, 2013.

[68] Sebastian Deorowicz, Marek Kokot, Szymon Grabowski, and Agnieszka Debudaj-Grabysz. Kmc 2: Fast and resource-frugal k-mer counting. *Bioinformatics*, 31(10):1569–1576, 2015.

[69] Chengxi Ye, Zhanshan Sam Ma, Charles H Cannon, Mihai Pop, and W Yu Douglas. Exploiting sparseness in de novo genome assembly. *BMC bioinformatics*, 13(6):1, 2012.

[70] Sara El-Metwally, Magdi Zakaria, and Taher Hamza. Lightassembler: fast and memory-efficient assembly algorithm for high-throughput sequencing reads. *Bioinformatics*, page btw470, 2016.

[71] Jason Pell, Arend Hintze, Rosangela Canino-Koning, Adina Howe, James M Tiedje, and C Titus Brown. Scaling metagenome sequence assembly with probabilistic de bruijn graphs. *Proceedings of the National Academy of Sciences*, 109(33):13272–13277, 2012.

[72] Rayan Chikhi and Guillaume Rizk. Space-efficient and exact de bruijn graph representation based on a bloom filter. *Algorithms for Molecular Biology*, 8(1):1, 2013.

[73] Kamil Salikhov, Gustavo Sacomoto, and Gregory Kucherov. Using cascading bloom filters to improve the memory usage for de bruijn graphs. *Algorithms for Molecular Biology*, 9(1):1, 2014.
[74] Rayan Chikhi, Antoine Limasset, Shaun Jackman, Jared T Simpson, and Paul Medvedev. On the representation of de bruijn graphs. In *International Conference on Research in Computational Molecular Biology*, pages 35–55. Springer International Publishing, 2014.

[75] Evangelos Georganas, Aydın Buluç, Jarrod Chapman, Steven Hofmeyr, Chaitanya Aluru, Rob Egan, Leonid Oliker, Daniel Rokhsar, and Katherine Yelick. Hipmer: an extreme-scale de novo genome assembler. In *Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*, page 14. ACM, 2015.

[76] Jintao Meng, Bingqiang Wang, Yanjie Wei, Shengzhong Feng, and Pavan Balaji. Swap-assembler: scalable and efficient genome assembly towards thousands of cores. *BMC bioinformatics*, 15(Suppl 9):S2, 2014.

[77] Jintao Meng, Sangmin Seo, Pavan Balaji, Yanjie Wei, Bingqiang Wang, and Shenzhong Feng. Swap-assembler 2: Optimization of de novo genome assembler at extreme scale. In *Parallel Processing (ICPP), 2016 45th International Conference on*, pages 195–204. IEEE, 2016.

[78] Michael Burrows and David J Wheeler. A block-sorting lossless data compression algorithm. 1994.

[79] Michael Roberts, Wayne Hayes, Brian R Hunt, Stephen M Mount, and James A Yorke. Reducing storage requirements for biological sequence comparison. *Bioinformatics*, 20(18):3363–3369, 2004.

[80] Steven L Salzberg, Adam M Phillippy, Aleksey Zimin, Daniela Puiu, Tanja Magoc, Sergey Koren, Todd J Treangen, Michael C Schatz, Arthur L Delcher, Michael Roberts, et al. Gage: A critical evaluation of genome assemblies and assembly algorithms. *Genome research*, 22(3):557–567, 2012.

[81] T. C Conway and A. J Bromage. Succinct data structures for assembling large genomes. *Bioinformatics*, 27(4):479, 2011.

[82] 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.

[83] 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.

[84] Erwan Drezen, Guillaume Rizk, Rayan Chikhi, Charles Deltel, Claire Lemaitre, Pierre Peterlongo, and Dominique Lavenier. Gatb: Genome assembly & analysis tool box. *Bioinformatics*, 30(20):2959–2961, 2014.

[85] Aleksey Zimin, Kristian A Stevens, Marc W Crepeau, Ann Holtz-Morris, Maxim Koriabine, Guillaume Marçais, Daniela Puiu, Michael Roberts, Jill L Wegrzyn,
Pieter J de Jong, et al. Sequencing and assembly of the 22-gb loblolly pine genome. *Genetics*, 196(3):875-890, 2014.

[86] Antoine Limasset, Guillaume Rizk, Rayan Chikhi, and Pierre Peterlongo. Fast and scalable minimal perfect hashing for massive key sets. *arXiv preprint arXiv:1702.03154*, 2017.

[87] Camille Marchet, Antoine Limasset, Lucie Bittner, and Pierre Peterlongo. A resource-frugal probabilistic dictionary and applications in (meta) genomics. *arXiv preprint arXiv:1605.08319*, 2016.

[88] Camille Marchet, Lolita Lecompte, Antoine Limasset, Lucie Bittner, and Pierre Peterlongo. A resource-frugal probabilistic dictionary and applications in bioinformatics. *arXiv preprint arXiv:1703.00667*, 2017.

[89] Lars Feuk, Andrew R Carson, and Stephen W Scherer. Structural variation in the human genome. *Nature Reviews Genetics*, 7(2):85-97, 2006.

[90] Pierre Peterlongo, Nicolas Schnel, Nadia Pisanti, Marie-France Sagot, and Vincent Lacroix. Identifying snps without a reference genome by comparing raw reads. In *International Symposium on String Processing and Information Retrieval*, pages 147-158. Springer Berlin Heidelberg, 2010.

[91] Jacob F Degner, John C Marioni, Athma A Pai, Joseph K Pickrell, Everlyne Nkadori, Yoav Gilad, and Jonathan K Pritchard. Effect of read-mapping biases on detecting allele-specific expression from rna-sequencing data. *Bioinformatics*, 25(24):3207-3212, 2009.

[92] Débora YC Brandt, Vitor RC Aguiar, Bárbara D Bitarello, Kelly Nunes, Jérôme Goudet, and Diogo Meyer. Mapping bias overestimates reference allele frequencies at the hla genes in the 1000 genomes project phase i data. *G3: Genes/Genomes/Genetics*, 5(5):931-941, 2015.

[93] Peter H Sudmant, Tobias Rausch, Eugene J Gardner, Robert E Handsaker, Alexej Abyzov, John Huddleston, Yan Zhang, Kai Ye, Goo Jun, Markus Hsi-Yang Fritz, et al. An integrated map of structural variation in 2,504 human genomes. *Nature*, 526(7571):75-81, 2015.

[94] 1000 Genomes Project Consortium et al. A global reference for human genetic variation. *Nature*, 526(7571):68-74, 2015.

[95] Agnieszka Danek, Sebastian Deorowicz, and Szymon Grabowski. Indexes of large genome collections on a pc. *PloS one*, 9(10):e109384, 2014.

[96] Guillaume Holley, Roland Wittler, and Jens Stoye. Bloom filter trie-a data structure for pan-genome storage. In *International Workshop on Algorithms in Bioinformatics*, pages 217-230. Springer, 2015.

[97] Benedict Paten, Adam M Novak, Jordan M Eizenga, and Erik Garrison. Genome graphs and the evolution of genome inference. *bioRxiv*, page 101816, 2017.
[98] Adam M Novak, Glenn Hickey, Erik Garrison, Sean Blum, Abram Connelly, Alexander Dilthey, Jordan Eizenga, MA Saleh Elmohamed, Sally Guthrie, André Kahles, et al. Genome graphs. bioRxiv, page 101378, 2017.

[99] Alexander Dilthey, Charles Cox, Zamin Iqbal, Matthew R Nelson, and Gil McVean. Improved genome inference in the mhc using a population reference graph. *Nature genetics*, 47(6):682–688, 2015.

[100] Ngan Nguyen, Glenn Hickey, Daniel R Zerbino, Brian Raney, Dent Earl, Joel Armstrong, W James Kent, David Haussler, and Benedict Paten. Building a pan-genome reference for a population. *Journal of Computational Biology*, 22(5):387–401, 2015.

[101] Richard Durbin. Efficient haplotype matching and storage using the positional burrows-wheeler transform (pbwt). *Bioinformatics*, 30(9):1266–1272, 2014.

[102] Adam M Novak, Erik Garrison, and Benedict Paten. A graph extension of the positional burrows-wheeler transform and its applications. In *International Workshop on Algorithms in Bioinformatics*, pages 246–256. Springer, 2016.

[103] Mingjie Wang, Yuzhen Ye, and Haixu Tang. A de bruijn graph approach to the quantification of closely-related genomes in a microbial community. *Journal of Computational Biology*, 19(6):814–825, 2012.

[104] Manuel Garber, Manfred G Grabherr, Mitchell Guttman, and Cole Trapnell. Computational methods for transcriptome annotation and quantification using rna-seq. *Nature methods*, 8(6):469–477, 2011.

[105] Bo Liu, Hongzhe Guo, Michael Brudno, and Yadong Wang. debga: read alignment with de bruijn graph-based seed and extension. *Bioinformatics*, page btw371, 2016.

[106] Gaëtan Benoit, Dominique Lavenier, Claire Lemaitre, and Guillaume Rizk. Bloocoo, a memory efficient read corrector. In *European Conference on Computational Biology (ECCB)*, 2014.

[107] Ilya Minkin, Anand Patel, Mikhail Kolmogorov, Nikolay Vyahhi, and Son Pham. Sibelia: a scalable and comprehensive synteny block generation tool for closely related microbial genomes. In *International Workshop on Algorithms in Bioinformatics*, pages 215–229. Springer, 2013.

[108] Shoshana Marcus, Hayan Lee, and Michael C Schatz. Splitmem: a graphical algorithm for pan-genome analysis with suffix skips. *Bioinformatics*, 30(24):3476–3483, 2014.

[109] Ilia Minkin, Son Pham, and Paul Medvedev. Twopaco: An efficient algorithm to build the compacted de bruijn graph from many complete genomes. *Bioinformatics*, page btw609, 2016.

[110] Timo Beller and Enno Ohlebusch. Efficient construction of a compressed de bruijn graph for pan-genome analysis. In *Annual Symposium on Combinatorial Pattern Matching*, pages 40–51. Springer, 2015.
[111] Hayan Lee and Michael C Schatz. Genomic dark matter: the reliability of short read mapping illustrated by the genome mappability score. *Bioinformatics*, 28(16):2097–2105, 2012.

[112] Zamin Iqbal, Mario Caccamo, Isaac Turner, Paul Flicek, and Gil McVean. De novo assembly and genotyping of variants using colored de bruijn graphs. *Nature genetics*, 44(2):226–232, 2012.

[113] Gustavo AT Sacomoto, Janice Kielbassa, Rayan Chikhi, Raluca Uricaru, Pavlos Antoniou, Marie-France Sagot, Pierre Peterlongo, and Vincent Lacroix. K is s plice: de-novo calling alternative splicing events from rna-seq data. *BMC bioinformatics*, 13(6):S5, 2012.

[114] Son K Pham and Pavel A Pevzner. Drimm-synteny: decomposing genomes into evolutionary conserved segments. *Bioinformatics*, 26(20):2509–2516, 2010.

[115] Ben Langmead and Steven L Salzberg. Fast gapped-read alignment with bowtie 2. *Nature methods*, 9(4):357–359, 2012.

[116] Antoine Limasset, Bastien Cazaux, Eric Rivals, and Pierre Peterlongo. Read mapping on de bruijn graphs. *BMC bioinformatics*, 17(1):237, 2016.

[117] Xiao Yang, Sriram P Chockalingam, and Srinivas Aluru. A survey of error-correction methods for next-generation sequencing. *Briefings in bioinformatics*, 14(1):56–66, 2013.

[118] Yongchao Liu, Jan Schröder, and Bertil Schmidt. Musket: a multistage k-mer spectrum-based error corrector for illumina sequence data. *Bioinformatics*, 29(3):308–315, 2013.

[119] Heng Li. Bfc: correcting illumina sequencing errors. *Bioinformatics*, page btv290, 2015.

[120] Li Song, Liliana Florea, and Ben Langmead. Lighter: fast and memory-efficient sequencing error correction without counting. *Genome biology*, 15(11):509, 2014.

[121] Hamid Mohamadi, Hamza Khan, and Inanc Birol. ntcard: A streaming algorithm for cardinality estimation in genomics data. *Bioinformatics*, page btw832, 2017.

[122] Dent Earl, Keith Bradnam, John St John, Aaron Darling, Dawei Lin, Joseph Fass, Hung On Ken Yu, Vince Buffalo, Daniel R Zerbino, Mark Diekhans, et al. Assemblathon 1: a competitive assessment of de novo short read assembly methods. *Genome research*, 21(12):2224–2241, 2011.

[123] Keith R Bradnam, Joseph N Fass, Anton Alexandrov, Paul Baranay, Michael Bechner, Inanç Birol, Sébastien Boisvert, Jarrod A Chapman, Guillaume Chauspuis, Rayan Chikhi, et al. Assemblathon 2: evaluating de novo methods of genome assembly in three vertebrate species. *GigaScience*, 2(1):10, 2013.
[124] Ellen M Leffler, Kevin Bullaughey, Daniel R Matute, Wynn K Meyer, Laure Segurel, Aarti Venkat, Peter Andolfatto, and Molly Przeworski. Revisiting an old riddle: what determines genetic diversity levels within species? *PLoS Biol*, 10(9):e1001388, 2012.

[125] Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A. Gurevich, Mikhail Dvorkin, Alexander S. Kulikov, Valery M. Lesin, Sergey I. Nikolenko, Son K. Pham, Andrey D. Prjibelski, Alex Pyshkin, Alexander Sirotkin, Nikolay Vyahhi, Glenn Tesler, Max A. Alekseyev, and Pavel A. Pevzner. SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. *Journal of Computational Biology*, 19(5):455-477, 2012.

[126] Aleksey V Zimin, Guillaume Marçais, Daniela Puiu, Michael Roberts, Steven L Salzberg, and James A Yorke. The masurca genome assembler. *Bioinformatics*, 29(21):2669-2677, 2013.

[127] Mark J Chaisson and Glenn Tesler. Mapping single molecule sequencing reads using basic local alignment with successive refinement (blasr): application and theory. *BMC bioinformatics*, 13(1):238, 2012.

[128] Kin Fai Au, Jason G Underwood, Lawrence Lee, and Wing Hung Wong. Improving pacbio long read accuracy by short read alignment. *PloS one*, 7(10):e46679, 2012.

[129] Mohammed-Amin Madoui, Stefan Engelen, Corinne Cruaud, Caroline Belser, Laurie Bertrand, Adriana Alberti, Arnaud Lemainque, Patrick Wincker, and Jean-Marc Aury. Genome assembly using nanopore-guided long and error-free dna reads. *BMC genomics*, 16(1):327, 2015.

[130] Marten Boetzer and Walter Pirovano. Sspace-longread: scaffolding bacterial draft genomes using long read sequence information. *BMC bioinformatics*, 15(1):211, 2014.

[131] Viraj Deshpande, Eric DK Fung, Son Pham, and Vineet Bafna. Cerulean: A hybrid assembly using high throughput short and long reads. In *International Workshop on Algorithms in Bioinformatics*, pages 349-363. Springer Berlin Heidelberg, 2013.

[132] Chen-Shan Chin, David H Alexander, Patrick Marks, Aaron A Klammer, James Drake, Cheryl Heiner, Alicia Clum, Alex Copeland, John Huddleston, Evan E Eichler, et al. Nonhybrid, finished microbial genome assemblies from long-read smrt sequencing data. *Nature methods*, 10(6):563-569, 2013.

[133] Job Dekker, Karsten Rippe, Martijn Dekker, and Nancy Kleckner. Capturing chromosome conformation. *science*, 295(5558):1306-1311, 2002.

[134] Jacob O Kitzman. Haplotypes drop by drop. *Nature biotechnology*, 34(3):296-298, 2016.

[135] Neil I Weisenfeld, Vijay Kumar, Preyas Shah, Deanna Church, and David B Jaffe. Direct determination of diploid genome sequences. *bioRxiv*, page 070425, 2016.
[136] Leena Salmela, Riku Walve, Eric Rivals, and Esko Ukkonen. Accurate self-correction of errors in long reads using de bruijn graphs. *Bioinformatics*, page btw321, 2016.
<figure><img src="image_72.png" /></figure>