Camille Marchet

Antoine Limasset

Lucie Bittner

Pierre Peterlongo
email: pierre.peterlongo@inria.fr

A resource-frugal probabilistic dictionary

Keywords: Bioinformatics, Sequences comparison, Genomics, Metagenomics, Data structures, Minimal Perfect Hash Functions, Indexing

published or not. The documents may come

Introduction

A genome or a chromosome can be seen as a word of millions characters long, written in a four letters (or bases) alphabet. Modern molecular genome biology relies on sequencing, where the information contained in a genome is chopped into small sequences (around one hundred bases), called reads. By providing millions of short genomic reads along with reasonable sequencing costs, high-throughput sequencing technologies [START_REF] Schuster | Next-generation sequencing transforms todays biology[END_REF] (HTS) introduced an era where data generation is no longer a bottleneck while data analysis is, as this amount of sequences needs to be pulled together in a coherent way. Thanks to HTS improvements, it is possible to sequence hundreds of single genomes and RNA molecules, giving insight to diversity and expression of the genes. HTS even allow to go beyond the study of an individual by sequencing different species/organisms from the same environment at once, going from genomics to metagenomics. This massive sequencing represents a breakthrough: for instance one now can access and directly investigate the majority of the microbial world, which cannot be grown in the lab [START_REF] Hug | A new view of the tree of life[END_REF]. However, because of the diversity and complexity of microbial communities, such experiments produce tremendous volumes of data, which represent a challenge for bioinformaticians to deal with. The fragmented nature of genomic information, shredded in reads, craves algorithms to organize and make sense of the data.

A fundamental algorithmic need is to be able to index read sets for a fast information retrieval. In particular, given the amount of data an experiment can produce, methods that scale up to large data sets are needed. In this paper we propose a novel indexation method, called the quasi-dictionary, a probabilistic data structure based on Minimal Perfect Hash Functions (MPHF). This technique provides a way to associate any kind of data to any piece of sequence from a read set, scaling to very large (billions of elements) data sets, with a low and controlled false positive rate.

A number of studies have focused on optimizing non-probabilistic text indexation, using for instance FM-index [START_REF] Ferragina | Indexing Compressed Text[END_REF], or hash tables. However, except the Bloomier filter [START_REF] Charles | Bloomier filters: A second look[END_REF], to the best of our knowledge, no probabilistic dictionary has yet been proposed for which the false positive or wrong answer rates are mastered and limited. The quasi-dictionary mimics the Bloomier filter solution as it enables to associate a value to each element from a set, and to obtain the value of an element with a mastered false positive probability if the element was not indexed. Existing published results in [START_REF] Charles | Bloomier filters: A second look[END_REF] indicate that the Bloomier filter and the quasi-dictionary have similar execution times, while our results tend to show that the quasi-dictionary uses approximately ten times less memory. Moreover, there are no available/free Bloomier filter yet implemented.

We propose two applications that use quasi-dictionary for indexing k-mers, enabling to scale up large (meta)genomic instances. As suggested by their names (short read connector counter and short read connector linker, as presented below), these applications have the ability to connect any read to either its estimated abundance in any read set or to a list of reads in any read set. A key point of these applications is to estimate read similarity using k-mers diversity only. This alignment-free approach is widely used and is a good estimation of similarity measure [START_REF] Dubinkina | Assessment of k-mer spectrum applicability for metagenomic dissimilarity analysis[END_REF].

Our first application, called short read connector counter (SRC counter), consists in estimating the number of occurrences of a read (i.e. its abundance) in a read set. This is a central point in high-throughput sequencing studies. Abundance is first very commonly used as indicator value for reads trimming: i.e. reads with relatively low abundance value are considered as amplification errors and/or sequencing errors, and these rare reads are generally removed before thorough analyses [START_REF] Kunin | Wrinkles in the rare biosphere: pyrosequencing errors can lead to artificial inflation of diversity estimates[END_REF][START_REF] Schirmer | Insight into biases and sequencing errors for amplicon sequencing with the illumina miseq platform[END_REF]. The abundance of reads is then interpreted as a quantitative or semi-quantitative metric: i.e. reads abundance is used as a measure of genic or taxon abundance, themselves very commonly used for comparisons of community similarity [START_REF] Amend | Quantifying microbial communities with 454 pyrosequencing: does read abundance count?[END_REF][START_REF] Kembel | Incorporating 16s gene copy number information improves estimates of microbial diversity and abundance[END_REF].

The second proposed application in this work, called short read connector linker (SRC linker), consists in providing a list of similar reads between read sets. We define the read set similarity problem as follows. Given a read set bank and a read set query, provide a similarity measure between each pair of reads b i × q j , with b i a read from the bank set and q j a read from the query set. Note that the bank and the query sets may refer to the same data set. Computing read similarity intra-read set or inter-read sets can be performed by a general purpose tool, such as those computing similarities using dynamic programming, and using heuristic tools such as BLAST [START_REF] Altschul | Basic local alignment search tool[END_REF]. However, comparing all versus all reads requires a quadratic number of read comparisons, leading to prohibitive computation time, as this is shown in our proposed results. There exist tools dedicated to the computation of distances between read sets [START_REF] Benoit | Multiple Comparative Metagenomics using Multiset k-mer Counting[END_REF][START_REF] Maillet | COMMET: comparing and combining multiple metagenomic datasets[END_REF][START_REF] Maillet | Compareads: comparing huge metagenomic experiments[END_REF], but none of them can provide similarity between each pair of reads b i × q j . Otherwise, some tools such as starcode [START_REF] Zorita | Starcode: sequence clustering based on all-pairs search[END_REF] are optimized for pairwise sequence comparisons with mainly the aim of clustering DNA barcodes. As shown in results, such tools also suffer from quadratic computation time complexity and thus do not scale up data sets composed of numerous reads.

Availability and license Our proposed tools SRC counter and SRC linker were developed using the GATB library [START_REF] Drezen | GATB: Genome Assembly & Analysis Tool Box[END_REF]. They may be used as stand alone tools or as libraries. They are licensed under the GNU Affero General Public License version 3 and can be downloaded from http://github.com/GATB/short read connector. Also licensed under the GNU Affero General Public License version 3, the quasi-dictionary can be downloaded from http://github.com/pierrepeterlongo/quasi dictionary.

Methods

We first recall basic notations: a k-mer is a word of length k on an alphabet Σ. Given a read set1 R, a k-mer is said solid in R with respect to a threshold t if its number of occurrences in R is bigger or equal to t. Let |w| denote the length of a word w ∈ Σ * and |R| denote the number of elements contained in R.

Quasi-dictionary index

In the following, we present our indexing solution. Based on this solution, two applications are proposed sections 1.3 and 1.4.

The index we propose associates each solid k-mer from a read set R to a unique value in [0, N -1], with N being the total number of solid k-mers in R. Ideally, when querying a non indexed k-mer (i.e. a non solid k-mer or a k-mer absent from R) the index returns no value. In our proposal, a non indexed k-mer may be associated to a value in [0, N -1] with a probability p > 0. This is why we refer to our index as the quasi-dictionary, since it is a probabilistic index. However, note that querying any indexed k-mer provides a unique and deterministic answer.

We define the quasi-dictionary as follows :

Given a static set composed of N distinct elements, a quasi-dictionary is composed of two structures: a minimal perfect hash function MPHF (see for instance [START_REF] Belazzougui | Cache-oblivious peeling of random hypergraphs[END_REF]) and a table of fingerprints F ingerP rints.

The MPHF for S is a function such that:

∀e ∈ S, MPHF(e) = i ∈ [0, N -1] ∀e 1 , e 2 ∈ S, (MPHF(e 1 ) = MPHF(e 2 )) ⇔ (e 1 = e 2 )
The fingerprint table for S is composed of N elements. It assigns to each element form S an integer value in [0, 2 f -1], with f the size of the fingerprint in bits. This table is used to verify the membership of an element e to the indexed set of elements using the MPHF. When S represents a set of k-mers, the fingerprint table uses f ≤ 2k since a k-mer can be coded as a word of 2k bits. The false positive rate is then

2 (2k-f ) -1 2 2k ≈ 1 2 f .
1.2 Indexing solid k-mers using a quasi-dictionary Algorithm 1: create quasidictionary

Data: Read set R, k ∈ N, t ∈ N, f ∈ N Result:
A quasi-dictionary QD 1 k-mer set K = get solid kmers(R, k, t) ; 2 QD.M P HF = create M P HF (K) ; 3 foreach k-mer w in K do Algorithm 1 presents the construction of the quasi-dictionary. The set of solid k-mers (algorithm 1, line 1) is obtained using the DSK [START_REF] Rizk | DSK: K-mer counting with very low memory usage[END_REF] method. The MPHF (algorithm 1, line 2) is computed using the MPHF library2 .

The fingerprint of a word w (algorithm 1, line 5) is obtained thanks to a hashing function create f ingerprint :

Σ |w| → [0, 2 f -1],
with f ≤ 2k. In practice we chose to use a xor-shift [START_REF] Marsaglia | Xorshift rngs[END_REF] for its efficiency in terms of throughput and hash distribution. The querying of a quasi-dictionary with a word w is straightforward, as presented in Algorithm 2. The index of w is retrieved using the MPHF. Then the fingerprint stored for this index is compared to the fingerprint of w. If they differ, then w is not indexed and the -1 value is returned. If they are equal, the value index ≥ 0 is returned. Note that two distinct words have the same fingerprint with a probability ≈ 1 2 f . It follows that there is a probability ≈ 1 2 f that the quasi-dictionary returns a false positive value despite the fingerprint checking, i.e. an index = -1 for a non indexed word. On the other hand, the index returned for an indexed word is the correct one. In practice we use f = 12 that limits the false positive rate to ≈ 0.02%. Note that our implementation authorizes any value f ≤ 64.

DNA strands DNA molecules are composed of two strands, each one being the reverse complement3 of the other. As current sequencers usually do not provide the strand of each sequenced read, each indexed or queried k-mer should be considered in the Quasi dictionary structure composed of a MPHF and a table of fingerprints. The fingerprints are obtained by hashing the corresponding k-mers. Thanks to the MPHF a unique index is associated to each k-mer where the fingerprint is stored. The fingerprint is associated to the value we want to link to the k-mer. When a k-mer is queried the fingerprint obtained is checked against the fingerprint stored to limit false positives.

forward or in the reverse complement strand. This is why, in the proposed implementations, we index and query only the canonical representation of each k-mer, which is the lexicographically smaller word between a k-mer and its reverse complement.

Time and memory complexities

Our MPHF implementation has the following characteristics. The structure can be constructed in O(N ) time and uses ≈ 4 bits by elements. We could use parameters limiting memory fingerprint to less than 3 bits per element, but we chose parameters to greatly speed up MPHF construction and query. The fingerprint table is constructed in O(N ) time, as the create f ingerprint function runs in O(1). This table uses exactly N × f bits. Thus the overall quasidictionary size, with f = 12 is ≈ 16 bits per element. Note that this value does not take into account the size of the values associated to each indexed element.

The querying of an element is performed in constant time and does not increase memory complexity.

Approximating the number of occurrences of a read in a read set

Algorithm 3: SRC counter: Quasi-dictionary used for counting k-mers

Data: Read set B, read set Q, k ∈ N, t ∈ N, f ∈ N Result: For each read from Q, its k-mer similarity with set B 1 quasi-dictionary QD = create quasidictionary(B, k, t, f ) ; 2 create a table count composed of N integers a ; 3 foreach Solid k-mer w from B do 4
count[query quasidictionary(w)] = number of occurrences of w in B; 5 foreach read q in Q do 6 create an empty vector count q; 7 foreach k-mer w in q do 8 if query quasidictionary(w) ≥ 0 then 9 add count[query quasidictionary(w)] to count q ; 10 Output the q identifier, and (mean, median, min and max values of count q); a with N the number of solid k-mers from B As presented in Algorithm 3, we propose a first straightforward application using the quasi-dictionary. This application is called SRC counter for short read connector counter. It approximates the number of occurrences of reads in a read set.

Two potentially equal read sets B and Q are considered. The indexation phase works as follows. Each solid k-mer of B is indexed using a quasi-dictionary. A thirdparty table named count stores the counts of indexed k-mers. Elements of this table are accessed via the quasi-dictionary index value of indexed items (Algorithm 3 lines 4 and 9). The number of occurrences of each solid k-mer from B (line 4) is obtained from DSK output, used during the quasi-dictionary creation (line 1). Then starts the query phase. Once the count table is created, for each read q from set Q, the count of all its k-mers indexed in the quasi-dictionary are recovered and stored in a vector (lines 8 and 9). Finally, collected counts from k-mers from q are used to output an estimation of its abundance in read set B. The abundance is approximated using the mean number of occurrences of k-mers from q, to supplement we output the median, the min and the max number of occurrences of k-mers from q. In rare occasions, false positives of the method can lead to an over-estimation of the count.

This algorithm is extremely simple. In addition to the quasi-dictionary creation time and memory complexities, it has a constant memory overhead (8 bits by element in our implementation) and it has an additional O( Q∈Q |Q|) time complexity.

Identifying similar reads between read sets or inside a read set

Algorithm 4: SRC linker: Quasi-dictionary used for identifying read similarities Output the id of q and each b tg id associate to its count from targets table;

Data: Read set B, read set Q, k ∈ N, t ∈ N, f ∈ N Result: For each read from Q, its k-mer similarity with each read from set B 1 quasi-dictionary QD = create quasidictionary(B, k, t, f ) ; 2 create a table ids composed of N vectors of integers a ; 3 foreach read b in B do
a with N the number of solid k-mers from B b In practice only tg id whose count value is higher or equal to a user defined threshold are output

Our second proposal, called SRC linker for short read connector linker, compares reads from two potentially identical read sets B and Q. For each read q from Q, a similarity measure with reads from B is provided.

The similarity measure we propose for a couple of reads q × b is the number of positions on q that is covered by at least a k-mer that also occur on b. Note that this measure is not symmetrical as one does not verify that the k-mers do not overlap on b.

The indexation phase of SRC linker works as follows. A quasi-dictionary is created and a third-party table ids of size N is created. Each element of this table stores for a solid k-mer w from B a vector containing the identifiers of reads from B in which w occurs. See lines 2 to 7 of Algorithm 4.

The query phase (lines 8 to the end of Algorithm 4) is straightforward. In practice, for each targeted read b j in B we remind the ending position of the last shared k-mer on q with b j denoted by next f ree position in the Algorithm 4. Given a new shared k-mer, the number of positions that was not already covered by another shared k-mer is added to the similarity measure (line 17 of Algorithm 4).

Once all k-mers of a read q are treated, the identifier of q is output and for each read b j from B its identifier is output together with the number of shared k-mers with q. In practice, in order to avoid quadratic output size and to focus on similar reads, only reads sharing a number of k-mers higher or equal to a user defined threshold are output.

In addition to the quasi-dictionary data structure creation, considering a fixed read size, Algorithm 4 has O(N × m) memory complexity and a O(N + Q∈Q |Q| × m) time complexity, with m the average number of distinct reads from B in which a kmer from Q occurs. In the worst case m = N , for instance with

B = Q = A |read| N .
In practice, in our tests as well as for real sets composed of hundred of million reads, m is limited to ≈ 2.22.

Storing read identifiers on disk Storing the read identifiers as proposed in Algorithm 4 presents important drawbacks as it requires a large amount of RAM. In order to get rid of this limitation we propose a disk version of this algorithm, in which the table ids is stored on disk. As shown in Algorithm 5 (see Appendix), the algorithmic solution is not straightforward as one needs to know for each indexed k-mer w its number of occurrences in the read set B plus the number of occurrences of k-mers = w from B (false positives) that have the same quasi-dictionary index.

This disk based solution enables to scale up very large instances with frugal RAM needs, at the price of a longer computation time, as show in results.

Results

This section presents results about the fundamental quasi-dictionary data structure and about potential applications derived from its usage. To this end, we use a metagenomic Tara Oceans [START_REF] Karsenti | A holistic approach to marine Eco-systems biology[END_REF] read set ERR599284 composed of 189,207,003 reads of average size 97 nucleotides. From this read set, we created six sub-sets by selecting first 10K, 100K, 1M, 10M, 50M and 100M reads (with K meaning thousand and M meaning million).

Tests were performed on a linux 20-CPU nodes running at 2.60 GHz with an overall of 252 GBytes memory.

Paper Submitted to PSC

SRC counter tests and performances

We first provide SRC counter results enabling to evaluate the gain of our proposed data structure when compared to a classical hash table. Secondly we provide results that enable to estimate the impact of false positives on results. 1. Wallclock time and memory used by the SRC counter algorithm for creating and for querying the quasi-dictionary using the default fingerprint size f = 8 (denoted by "QD") and the C++ unordered map, denoted by "Hash". Column "k-mer count time" indicates the time DSK spent counting k-mers. Tests were performed using k = 31 and t = 1 (all k-mers are solid). The query read set was always the 10M set. We additionally provide memory results using the quasi-dictionary with a fingerprint size f = 62 (denoted by "QD62"). Construction and query time for QD62 are not shown as they are almost identical to the QD ones. On the full data set, using a classical hash table, the memory exceeded the maximal authorized machine limits (252 GB).

SRC counter performances compared to standard hash table index

We tested the SRC counter performances by indexing iteratively the six read subsets plus the full ERR59928 set, each time querying reads from set 10M. We compared our solution performances with a classical indexation scheme done using the C++ unordered map hash table. Results are presented in Table 1. These results show that the quasi-dictionary is much faster to compute than this hash table solution, in particular because of parallelisation. Moreover, the quasi-dictionary memory footprint is ≈ 13 times smaller on large enough instances (10 million indexed reads or more). These results show that the hash table is not a viable solution scaling up current read sets composed of several billions k-mers. Results also highlight the fact that the query is fast and only slightly depends on the number of indexed elements.

Importantly, using a fingerprint large enough (here f = 62 for k-mers of length k = 31), we can force the quasi-dictionary to avoid false positives. As expected, the quasi-dictionary data structure size increases with the size of f but interestingly, on this example and as shown in Table 1, the size of the quasi-dictionary with f = 62 remains in average 4 times smaller than the size of the hash-table on large problem instances. Keeping in mind that the quasi-dictionary is faster to construct and to query, the usage of this data structure avoiding false positives presents only advantages compared to the hash table usage for indexing a static set. However, one should recall that this is true because we are using an alphabet of size four, so any 31-mer on the alphabet {A, C, G, T } can be assigned to a unique value in [0, 2 62 -1] and vice versa. With larger alphabets such as the amino-acids or the Latin ones, the usage of a hash table is recommended if false positives are not tolerated.

Approximating false positives impact

We propose an experiment to assess the impact on result quality when using a probabilistic data structure instead of a deterministic one for estimating read abundances. We used the read set 100M both for the indexation and for the querying, thus providing an estimation of the abundance of each read in its own read set. We made the indexation using k = 31, c = 2 and f = 8. Note that, with c = 2 only k-mers seen twice or more in the set are solid and thus are indexed. In this example only 756,804,245 k-mers are solid among the [START_REF] Belazzougui | Cache-oblivious peeling of random hypergraphs[END_REF]191,190,377 distinct k-mers present in the read set. This means that during the query, 85.4% of queried k-mers are not indexed. This enables to measure the impact of the quasi-dictionary false positives. We applied the count algorithm as described in Algorithm 3, and the tuned version using a hash table instead of a quasi-dictionary. We analyzed the count output composed of the average number of occurrences of k-mers of each read in the 100M read set.

Because of the quasi-dictionary false positives, results obtained using this structure are an over-estimation of the real result. Thus, we computed for each read the observed difference in the counts between results obtained using the quasi-dictionary implementation and the hash table implementation. The max over-approximation is 26.9, and the mean observed over-approximation is 7.27 × 10 -3 with a 3.59 × 10 -3 standard deviation. Thus, as the average estimated abundance of each read which is ≈ 2.22, the average count over-estimation represents ≈ 0.033% of this value. Such divergences are negligible.

Identifying similar reads

We set a benchmark of our method with comparisons to state of the art tools that can be used in current pipelines for the read similarity identification presented in this paper. We compared our tool with the classical method BLAST [START_REF] Altschul | Basic local alignment search tool[END_REF] (version 2.3.0), with default parameters. BLAST is able to index big data sets, and consumes a reasonable quantity of memory, but the throughput of the tool is relatively low and only small data sets were treated within the timeout (10h, wallclock time). We also included two broadly used mappers in the comparison. We used Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] (version 2.2.7), and BWA [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] (version 0.7.10) in any alignment mode (-a mode in Bowtie2, -N for BWA) in order to output all alignment found instead of the best ones only. Both tools are not well suited to index large set of short sequences nor to find all alignments and therefore use considerably more resources than their standard usage.

We also compared SRC linker to starcode (1.0), that clusters DNA sequences by finding all sequences pairs below a Levenshtein distance metric. One should notice that benchmark comparisons with tools as starcode is unfair as such tool provides much more precise distance information between pair of reads than SRC linker and performs additional task as clustering. However, our benchmark highlights the fact that such approaches suffer from intractable number of read comparisons, as demonstrated by presented results.

We focused on a practical use case for which our method could be used, namely retrieving similarities in a read set against itself. We used default SRC linker parameters (k = 31, f = 12, c = 2). Because of the limitations of the methods we used for the benchmark, reported in Table 2, we could compare against all methods only up to 1M reads. BWA performed better than the two other tools in terms of memory, being able to scale up to 10M reads, while Bowtie2 and BLAST could only reach 1M reads comparison. On this modest size of read set, we show that we are already ahead both in terms of memory and time. However the gap between our approach and others increases with the amount of data to process. Dealing with the full Tara 2. CPU time and memory consumption for indexing and querying a data set versus itself. We set a timeout of 10h. BLAST crashed for 10M data set, Bowtie2 reached the timeout we set with more than 200h (CPU) for 10M reads. BWA performs best among the mappers, reaching the timeout for 100M reads (more than 200h (CPU) on this data set). On the 100M data set, starcode reached the timeout. Only SRC linker finished on all data sets. On the full data set, it lasted an order of magnitude comparable to what BWA performed on only 10M.

data set reveals the specificity of our approach (Table 2) that requires low resources in comparison to others and is able to deal with bigger data sets. 3. Multithreading and disk performances. The full read set was used to detail the performances of the RAM and Disk algorithm on a large data set. We used default parameters k = 31, f = 12, c = 2. Times are wallclock times.

Finally, we highlight that we provide a parallelised tool (10× speedup for the index and 17× speedup for query for RAM algorithm as shown in Table 3) on the contrary to classical methods that are partly-parallelised as only the alignment step is well suited for parallelisation. The disk version does not fully benefit from multiple cores since the bottleneck is disk access. The main interest of this technique is a highly reduced memory usage at the price of an order of magnitude lower throughput, as presented Table 3.

Discussions and conclusion

In this contribution, we propose a new indexation scheme based on a Minimal Perfect Hash Function (MPHF) together with a fingerprint value associated to each indexed element. Our proposal is a probabilistic data structure that has similar features than Bloomier filters, with smaller memory fingerprint. This solution is resource-frugal (we have shown experiments on sets containing more than eight billion elements indexed in ≈ 3 hours and using less than 25GB RAM) and opens the way to new (meta)genomic applications. As proofs of concept, we proposed two novel applications: SRC counter and SRC linker. The first estimates the abundance of a sequence in a read set. The second detects similarities between pair of reads inter or intra-read sets. These applications are a start for broader uses and purposes.

Two main limitations of our proposal due to the nature of the data structure can be pointed out. Firstly, compared to standard hash tables, our indexing data structure presents an important drawback: the exact set of keys to index has to be defined during the data structure creation and it has to be static. This may be a limitation for non fixed set of keys. Moreover, our data structure can generate false positives during query. Even with the proposed false positive ratio limited to ≈ 10 -2 % with defaults parameters, this may be incompatible with some applications. However we can force our tools to avoid false positives by using as a fingerprint the key itself. Interestingly, this still provides better time and memory performances than using a standard hash table in the DNA k-mer indexing context, with k = 31, which is a very common value for read comparisons [START_REF] Benoit | Multiple Comparative Metagenomics using Multiset k-mer Counting[END_REF]. Secondly, one should notice that our indexation proposal saves space regarding the association between an element and a specific array offset (if the element was indexed). However, our proposal does not limit the space needed for storing the value associated to each indexed element. Thus, with respect to classical hash tables, the memory gain is limited in problem instances in which large values are associated to each key. Indeed, in this case, the memory footprint is mainly due to the value over the indexing scheme. In order to benefit from our proposal even in such cases we proposed an application example in which the values are stored on disk. However, our approach is namely designed for problems where a huge number of elements to index are at stake, along with a small quantity of information to match with.

We could improve our technique to recognize key from the original set, using a technique from the hashing field [START_REF] Kirsch | Less hashing, same performance: Building a better Bloom filter[END_REF] or from the set representation field [START_REF] Belazzougui | Compressed static functions with applications[END_REF]. In such framework, a set can be represented with less memory than the sum of the memory required by the keys. We could thus hope being able to represent a non-probabilistic dictionary without storing keys. Otherwise, we could use the hashing information to achieve a smaller false positive rate with the same or a reduced memory usage. The main challenge will be to keep fast query operation for such complex data structure.

The results we provided show that alignment-based approaches do not scale when it comes to find similar reads in data sets composed of millions of sequences. The fact that HTS data count rarely less than millions reads justifies our approach based on k-mer similarity. Moreover our approach is more straightforward and requires less parameters and heuristics than mapping approaches, that can sometimes turn them into blackboxes. However, such an approach remains less precise than mapping, since the k-mer order is not taken into account and is less sensitive because of the fixed size of k. An important future work will be to evaluate the differences between matches of our pseudo-alignment and matches of well-known and widely used tool as BLAST.

Our tools property of enabling the test of a read set against itself opens the door to applications such as read clustering. Latest sequencing technologies, called Third Generation Sequencers (TGS), provide longer reads [START_REF] Sharon | A single-molecule long-read survey of the human transcriptome[END_REF][START_REF] Tilgner | Defining a personal, allelespecific, and single-molecule long-read transcriptome[END_REF] (more than a thousand bases instead of a few hundreds for HTS). With previous HTS short reads, de novo approaches to reconstruct DNA or RNA molecules were using assembly [START_REF] Grabherr | Full-length transcriptome assembly from RNA-Seq data without a reference genome[END_REF][START_REF] Robertson | De novo assembly and analysis of RNA-seq data[END_REF], based on de Bruijn graphs. For RNA, these TGS long reads mean a change of paradigm as assembly is no more necessary, as one read is long enough to represent one full-length molecule. The important matter becomes to segregate families of RNA molecules within a read set, a purpose our approach could be designed for.

Furthermore, the methods we provide have straightforward applications examples in biology, such as the building of sequences similarity networks (SSN) [START_REF] Atkinson | Using sequence similarity networks for visualization of relationships across diverse protein superfamilies[END_REF] using SRC linker. SSN are extremely useful for biologists because, in addition to allowing a user-friendly visualization of the genetic diversity from huge HTS data sets, they can be studied analytically and statistically using graph topology metrics. SSN have recently been adapted to address an increasing number of biological questions investigating both patterns and processes: e.g. population structuring [START_REF] Forster | Testing ecological theories with sequence similarity networks: marine ciliates exhibit similar geographic dispersal patterns as multicellular organisms[END_REF][START_REF] Fondi | Every gene is everywhere but the environment selects: Global geo-localization of gene sharing in environmental samples through network analysis[END_REF]; genomes heterogeneity [START_REF] Boon | Studying genome heterogeneity within the arbuscular mycorrhizal fungal cytoplasm[END_REF]; microbial complexity and evolution [START_REF] Corel | Bapteste: Network-Thinking: Graphs to Analyze Microbial Complexity and Evolution[END_REF]; microbiome adaptation [START_REF] Bapteste | Evolution of genetic diversity using networks: the human gut microbiome as a case study[END_REF][START_REF] Völkel | Read networks and k-laminar graphs[END_REF] or to explore the microbial dark matter [START_REF] Lopez | Highly divergent ancient gene families in metagenomic samples are compatible with additional divisions of life[END_REF]. In metagenomic microbial studies, SSN offer an alternative to classical and potentially biased methods, and thus facilitate large-scale analyses and hypotheses generation, while notably including unknown/dark matter sequences in the global analysis [START_REF] Forster | Testing ecological theories with sequence similarity networks: marine ciliates exhibit similar geographic dispersal patterns as multicellular organisms[END_REF][START_REF] Lopez | Highly divergent ancient gene families in metagenomic samples are compatible with additional divisions of life[END_REF]. Currently SSN are built upon general purposes tools such as BLAST. They thus hardly scale up large data sets. A future work will consist in checking the feasibility of applying SRC linker for constructing SSN and, in case of success, to use it on large SSN problem instances on which other classical tools cannot be applied.

Algorithm 5: SRC linker Disk: Quasi-dictionary used for identifying read similarities Data: Read set B, read set Q, k ∈ N, t ∈ N, f ∈ N Result: For each read from Q, its k-mer similarity with each read from set B 1 quasi-dictionary QD = create quasidictionary(B, k, t, f ) ; 2 create a table ids composed of N integers a all valued to 0; Output the id of q and each b tg id associate to its count from targets table ;   a with N the number of solid k-mers from B b In practice only tg id whose count value is higher or equal to a user defined threshold are output

4 index=

 4 QD.M P HF (w); 5 QD.F ingerP rints[index] = create f ingerprint(w, f ); 6 return QD;

Algorithm 2 :

 2 query quasidictionary Data: Quasi-dictionary QD, word w Result: A unique value in [0, N -1] (with N the number of indexed elements) or -1 if w detected as non indexed 1 index = QD.M P HF (w); 2 if index ≥ 0 and QD.F ingerP rints[index] = create f ingerprint(w) then 3 return index; 4 return -1;

Figure 1 .

 1 Figure1. Quasi dictionary structure composed of a MPHF and a table of fingerprints. The fingerprints are obtained by hashing the corresponding k-mers. Thanks to the MPHF a unique index is associated to each k-mer where the fingerprint is stored. The fingerprint is associated to the value we want to link to the k-mer. When a k-mer is queried the fingerprint obtained is checked against the fingerprint stored to limit false positives.

4

  foreach k-mer w in b do 5 index = query quasidictionary(w); 6 if index ≥ 0 then 7 add identifier of b to vector ids[index] ; 8 foreach read q in Q do 9 create a hash table associating targets (target read id) to couple(next f ree position, count); 10 foreach i in [0, |q| -k] do 11 w = k-mer occurring position i in q; 12 index = query quasidictionary(w); 13 if index ≥ 0 then 14 foreach tg id in vector ids[index] do 15 if targets[tg id] is empty then 16 targets[tg id].next f ree position = 0 17 targets[tg id].count+ = max (k, i + k -targets[tg id].next f ree position) targets[tg id].next f ree position = i + k 18

3 foreach read b in B do 4 foreach k-mer w in b do 5 index 6 if index ≥ 0 then 7 add 1

 45671 = query quasidictionary(w); to ids[index]; 8 foreach Solid k-mer w from B do 9 index = query quasidictionary(w); 10 if index ≥ 0 then 11 count = ids[index]; 12 ids[index] = T emporary F ile.position; 13 write count + 1 '0' on T emporary F ile; 14 foreach read b in B do 15 foreach k-mer w in b do 16 index = query quasidictionary(w); 17 if index ≥ 0 then 18 position = ids[index]; 19 T emporary F ile.goto(position); 20 write id of b in place of the first 0 found; 21 foreach read q in Q do 22 create a hash table targets (target read id) → couple(next f ree position, count); 23 foreach i in [0, |q| -k] do 24 w = k-mer occurring position i in q; 25 index = query quasidictionary(w); 26 if index ≥ 0 then 27 position = ids[index]; 28 T emporary F ile.goto(position); 29 read from T emporary F ile and put in vector V all integer until a 0 is found; 30 foreach tg id in vector V do 31 if targets[tg id] is empty then 32 targets[tg id].next f ree position = 0 33 targets[tg id].count+ = max (k, i + k -targets[tg id].next f ree position) targets[tg id].next f ree position = i + k 34

Note that formally R should be denoted as a "collection" instead of a "set" as a read may appear twice or more in R. However, to make the reading easier, we use in this manuscript the term "set" usually employed for describing HTS outputs.

https://github.com/rizkg/BooPHF, commit number 852cda2

The reverse complement of a DNA sequence is the palindrome of the sequence, in which A and T are swapped and C and G are swapped. For instance the reverse complement of ACCG is CGGT .

http://www.ebi.ac.uk/ena/data/view/ERR599280

Acknowledgments

This work was funded by French ANR-12-BS02-0008 Colib'read project. We thank the GenOuest BioInformatics Platform that provided the computing resources necessary for benchmarking. We warmly thank Guillaume Rizk and Rayan Chikhi for their work on the MPHF and for their feedback on the preliminary version of this manuscript.