Christophe Vroland

Mikael Salson

Sébastien Bini

Hélène Touzet
email: helene.touzet@univ-lille1.fr

Mikaël Salson

Approximate search of short patterns with high error rates using the 010 lossless seeds

Keywords: approximate string matching, lossless approximate seeds, Levenshtein distance, computational biology

come

Introduction

We consider the approximate string matching problem, where an error between a pattern P and a text T is measured by the Levenshtein distance: an error is either a substitution, an insertion, or a deletion. The first approach to solve this problem was to perform dynamic programming with an O(|T | × |P |) algorithm [START_REF] Wagner | The string-to-string correction problem[END_REF]. This kind of approach has been improved by using bit-parallelism or finite automaton (for instance [START_REF] Wu | Fast text searching: allowing errors[END_REF]). However it is still not practical to search large databases of sequences or a full genome.

Rather than traversing all the sequences, it is possible to start by filtering candidate regions which show a high-enough similarity between portions of the text and of the pattern. Such portions are called seeds. Their occurrences must then be extended to check if the pattern occurs within k errors. This seed-andextend approach is now widely used to solve the approximate pattern matching problem. The seed definition can be very diverse and depending on the type of application (size of the alphabet, rate of errors, . . . ) the preferred seeds may differ. The most widely known seed is the contiguous seed, which consists in a substring that must be shared between the pattern and the text. For instance it is used in the Blast software for homology search in bioinformatics [START_REF] Altschul | Basic local alignment search tool[END_REF]. Depending on the error rate and on the length of the contiguous seed, those seeds may be lossy, meaning that the query algorithm is a heuristic that can lead to false negatives. The classical pigeonhole principle avoids this pitfall by splitting the pattern in exactly k + 1 parts, which are searched independently as instances of contiguous seeds. This is an example of lossless seed, that never discards an actual occurrence. A recent example of this method, using a modified Burrows-Wheeler transform is shown in [START_REF] Petri | Efficient indexing algorithms for approximate pattern matching in text[END_REF]. However the more errors we allow, the shorter the parts will be and therefore the more potential occurrences we may have. Thus the filtration efficiency becomes lower with higher values of k. This limitation can be partially overcome by the usage of spaced seeds [START_REF] Ma | Patternhunter -faster and more sensitive homology search[END_REF] or subset seeds [START_REF] Kucherov | A unifying framework for seed sensitivity and its application to subset seeds[END_REF] that offer a better selectivity/specificity trade off than contiguous seeds [START_REF] Brown | A tutorial of recent developments in the seeding of local alignment[END_REF][START_REF] Keich | On spaced seeds for similarity search[END_REF]. One main drawback is that those seeds tolerate substitutions only. They do not improve the sensitivity over contiguous seeds with insertions and deletions. Attempts to generalize spaced seeds to the full Levenshtein distance were presented in [START_REF] Burkhardt | One-Gapped q-Gram Filters for Levenshtein Distance[END_REF] for example, with a strong restriction on the total number of errors.

Other approaches are hybrid approaches that combine seed filtering and exact search. The pattern is split in non-overlapping parts that can be searched with a given number of errors. Each method differs in the number of pieces in the patterns and on the maximal number of errors in each piece. For instance, Navarro and Baeza-Yates [START_REF] Navarro | A Hybrid Indexing Method for Approximate String Matching[END_REF] designed a hybrid method which consists of splitting the pattern P in j = (|P | + k)/ log σ |T | parts, where σ is the alphabet size, and searching these parts with k j errors. This approach has also been used with a LZ-index or an FM-index [START_REF] Russo | Approximate string matching with compressed indexes[END_REF]. In [START_REF] Kärkkäinen | Faster filters for approximate string matching[END_REF], Kärkkäinen and Na introduded a filtering method based on suffix filters, where the pattern is cut in k + 1 factors and the set of seeds is the set of all suffixes with strong matches. In general, hybrid approaches are more flexible than contiguous or spaced seeds as they tolerate any kind of errors, without restriction on the positions. However, their filtration efficiency tends to decrease as the number of errors increases, or they involve factorizations where some parts have to be searched with multiple errors.

Despite this large body of research on seeding techniques for the approximate string matching problem, we believe that there is still room for improvement for small patterns with high error rates (10-20%). In this article we will present a new hybrid method where the pattern will be split in k + 2 non-overlapping parts, some of them being searched without error, while others are searched with exactly one error. This approach results in the definition of a new kind of lossless seeds, called 01 * 0 seeds, that can accomodate a large number of insertions and deletions 1 . In Section 2, we will first introduce the theoretical framework of these new seeds. Once the type of seeds has been chosen, there are different ways of using them. The ability to index them or to search them efficiently in a text index is an important criterion. To illustrate this, we will present two real-life applications of those seeds in a bioinformatics context. In Section 3, we will present an algorithm for the approximate pattern matching problem in large genomic sequences. In this application, 01 * 0 seeds are used in conjunction with a compressed full-text index, namely the FM-index [START_REF] Ferragina | Indexing compressed text[END_REF]. In Section 4, we will address the complementary problem of comparing a large collection of patterns against a relatively short DNA sequence such as a gene. For doing so, we will explain how to index the seeds themselves. Finally, in Section 5, we will discuss the application of the 01 * 0 seeds to the k-mismatch problem and compare them to spaced seeds.

2. The 01 * 0 seeds

Definition

Let Σ be a finite alphabet of size σ. Given two strings U and V of Σ * , define lev(U, V ) to be the Levenshtein distance between U and V . This is the minimum number of operations needed to transform U into V , where the only allowed operations are substitution of a single character and deletion or insertion of a single character. Each such operation is also called an error. From now on, we assume that a given natural number k corresponds to a maximum number of errors.

Let P be a pattern of length m over Σ. Using the pigeonhole principle, it is well-known that if P is partitioned into k + 1 parts, then every string U , such that lev(P, U ) ≤ k, contains at least one of these parts as a substring. The parts do not need to be of the same length. Similarly, if P is partitioned into k + 2 parts, denoted P 1 , . . . , P k+2 , then U should contain at least two disjoint parts of P . The following lemma allows to push the analysis further. It is indeed possible to request that these two parts be separated by parts with exactly one error.

Lemma 1. Let U be a string of Σ * such that lev(P, U ) ≤ k. Then there exist i, j, 1 ≤ i < j ≤ k + 2, and U 1 , . . . , U j-i-1 of Σ * such that

1. P i U 1 . . . U j-i-1 P j is a substring of U , and 2. When j > i + 1, for each , 1 ≤ ≤ j -i -1, lev(P i+ -1 , U ) = 1.
Example 1. Assume k = 3. Given the pattern P = AACGTGAGGTAGGTTCCATG of length 20, we partition it into five parts, of equal length: P 1 = AACG, P 2 = TGAG, P 3 = GTAG, P 4 = GTTC, and P 5 = CATG. Consider three strings whose Levenshtein distance with P is 3: AACGGAGGTAAGTTCTCATG, AACGTAGGCAAGTTCCATG The pattern AACGTGAGGTAGGTTCCATG is on the top of each alignment. The two parts with no error (written P i and P j in the Lemma) are highlighted in dark grey, and the parts with 1 error surrounded by P i and P j are highlighted in light grey. and ATCGTGACGTAGGGTCCATG. For each string, we show in Figure 1 the parts that fulfil the conditions of Lemma 1.

P 1 P 2 P 3 P 4 P 5 A A C G | | | | A A C G 0 T G A G | | | -G A G 1 G T A G | | | G T A A 1 G T T C | | | | G T T C 0 -C A T G | | | | T C A T G P 1 P 2 P 3 P 4 P 5 A A C G | | | | A A C G T G A G | | | T -A G G T A G | | G C A A G T T C | | | | G T T C 0 C A T G | | | | C A T G 0 P 1 P 2 P 3 P 4 P 5 A A C G | | | A T C G T G A G | | | T G A C G T A G | | | | G T A G 0 G T T C | | | G G T C 1 C A T G | | | | C A T G 0
We rephrase Lemma 1 as a pure counting problem and establish its proof. Lemma 2. Let k be a natural number. Assume you have k + 2 containers numbered from 1 to k + 2, and y tokens with 0 ≤ y ≤ k. Then there exists two containers i and j, 1 ≤ i < j ≤ k + 2, such that 1. containers i and j are empty, and 2. for each , i < < j, the container contains exactly one token.

Proof. Let µ be the number of containers with at least two tokens, called the multiple containers. In total there are at least µ + 2 empty containers. The 2 comes from the fact that we have k tokens and k + 2 containers. Therefore there are at least µ + 1 pairs of empty containers i and j such that each container between i and j is non empty. Each of the µ multiple containers can lie between one distinct pair of empty containers, breaking the property for µ pairs of empty containers. Stated otherwise, we have at most µ pairs of empty containers with one intervening multiple container. This leaves us with µ + 1µ = 1 pair satisfying the property that the intervening containers (if any) have exactly one token.

As a consequence of Lemma 1, we can design a seeding framework for lossless filtering for the approximate pattern matching problem with k errors. To this end, we introduce some terminology that will be used in the remainder of the paper.

Definition 1. Let P = P 1 . . . P k+2 be a pattern divided into k + 2 parts. Then the 01 * 0 seed for P and k is the regular expression ∪ k+1 i=1 ∪ k+2 j=i+1 P i lev 1 (P i+1 ) . . . lev 1 (P j-1 ) P j

where lev 1 (u) denotes the set of strings whose Levenshtein distance with u is 1.

A subseed is the regular expression associated to a given pair (i, j):

P i lev 1 (P i+1 ) . . . lev 1 (P j-1 ) P j .
An instance of a seed, or of a subseed, is a string u of Σ * which is recognized by the seed, or the subseed. Given a text T on Σ * , an occurrence of the seed for P is a substring of T which is an instance of the seed. Therefore, an occurrence is characterized by its start position and its end position in T .

Filtration efficiency

The filtration efficiency is the primary criterion used to evaluate the performance of a seed. The thorough analysis of the average case for the 01 * 0 seeds, such as what is done in [START_REF] Baeza-Yates | Fast and practical approximate string matching[END_REF] for example, is beyond the scope of this paper. Having a formula for the expected numbers of occurences of subseed instances for a given pattern is also a non trivial task, since this formula would be highly dependent of the structure of the pattern and of its internal repeats. We leave it as an open question for future work. Here we present a simple lemma, that gives an upper bound for the probability of an occurrence of a subseed.

Lemma 3.

If the text T is a sequence of independent and identically distributed characters of Σ, then the probability to have an occurrence of the subseed P i . . . P j at a given position of the text is smaller than

1 σ pi × L(p i+1 ) × . . . × L(p j-1 ) × 1 σ pj
where p i , . . . , p j are the lengths of the parts P i , . . . , P j and L(p) = p(σ-1)

σ p + p σ p-1 + σ(p+1) σ p+1 .
Proof. By Definition 1, P i lev 1 (P i+1 ) . . . lev 1 (P j-1 ) P j gives the set of all possible instances of the subseed. The probability to observe P i and P j is 1 σ p i and 1 σ p j respectively. For each part P of length p in {P i+1 , . . . , P j-1 }, p(σ-1) σ p + p σ p-1 + σ(p+1) σ p+1 bounds the probability to observe lev 1 (P ). Indeed, p(1-σ) σ p is an upper bound for the probability of the set of words with one mismatch from P , p σ p-1 is an upper bound for the probability of the set of words with one deletion from P , and σ(p+1) σ p+1 is an upper bound for the probability of the set of words with one insertion from P . This lemma can be used to obtain an upper bound on the average number of occurrences found in a random text, by summing up all the probabilities of all subseeds. For example, for k = 3, a pattern of length m = 20 divided in five equal parts of length 4, a text of size 10 8 over the DNA aphabet {A, C, G, T } and the result is 6,747 occurrences.

Complementary to this lemma, we generated an independent and identically distributed random sequence of length 10 8 over the DNA alphabet as well as 100 patterns of length 20. We then searched for the 01 * 0 seeds for k = 3. For each pattern, we counted the total number of occurrences of the seeds in the text, including overlapping occurrences. The distribution is plotted in Table 1. The average number of observed occurrences per pattern is 6,665 (98.7% of our estimated upper bound).

For the purpose of comparison, we also estimated the filtration efficiency of several other types of seeds, for the same text and the same collection of patterns. First, we considered the classical pigeonhole principle, which guarantees lossless filtration. The pattern is split in k + 1 = 4 non-overlapping parts of length 5. In Table 1, we call it the k + 1 pigeonhole seed. The average number of occurrences is 385,651. We then analysed a k + 2 version of the pigeonhole principle: the pattern is divided in k + 2 = 5 non-overlapping parts of length 4, and two parts must match exactly. In a first case, we require that any two exact parts are separated by a distance compatible with at most k errors: the number of positions between two exact parts P i and P j (i < j) should range from

p i+1 + • • • + p j-1 -k to p i+1 + • • • + p j-1 + k.
In a second case, we consider a weakened form of Lemma 1: all parts in-between the two exact parts should contain exactly one error. As a consequence, the number of positions between P i and P j ranges from p i+1 +• • •+p j-1 -(j -i-1) to p i+1 +• • •+p j-1 +(j -i-1). We refer to this seed as the bounded k + 2 pigeonhole seed. It can be seen as an intermediate seed between the k + 2 pigeonhole seed and the 01 * 0 seeds and allows us to better understand the behavior of 01 * 0 seeds. In those two cases, the average number of occurrences is respectively of 86,145 and 44,746 (see Table 1, k + 2 pigeonhole and bounded k + 2 pigeonhole). These plots show that moving from k + 1 parts to k + 2 parts allows to improve the filtration efficiency by a factor greater than 4. The additional constraint on distances beween parts further decreases the number of occurrences by almost half. Finally, the 01 * 0 seed offers a new gain of a factor 6.7. Overall, this is more than 57 times better than the filtering rate achieved with the k + 1 pigeonhole principle. As a last case, we chose to divide the pattern in only three parts, of lengths 6, 7, and 7. Contrary to the four preceding seeds, this seed is lossy since it allows for some false negatives. The average number of occurrences is 36,207 (see Table 1, lossy seed). Interestingly, this value is close to the bounded k + 2 pigeonhole seed, but the former maintains full sensitivity. It is far from being as selective as 01 * 0 seed.

These empirical measurements show that the 01 * 0 seed is significantly more selective than all other seeds presented in Table 1. It provides orders of magnitude more efficient filtering. Of course, this higher selectivity comes at the price of some additional work to locate seeds in the text. Each seed is composed 1: Distribution of the number of occurrences for five different seeds on a dataset constituted by 100 patterns of size m = 20 in a random sequence of length 10 8 on the DNA alphabet. The 01 * 0 seed and the three pigeonhole seeds are designed for k = 3 errors. So for the 01 * 0 seed and the k + 2 pigeonhole seeds, the pattern is divided in 5 parts of length 4.

01 * 0 seed k + 1 k + 2
For the k + 1 pigeonhole seed, the pattern is divided in 4 parts of length 5. For the lossy seed, each pattern is split in 3 parts of lengths 6, 7, and 7, respectively. There are on average 26.85 occurrences of the whole pattern within 3 errors.

of (k+1)(k+2) 2 subseeds, and some parts of the subseed should be searched with errors. This task, however, is facilitated by the regular structure of the 01 * 0 seeds. Firstly, each part of the subseed contains exactly 0 or 1 error. Secondly, distinct subseeds may share a common exact prefix or a common exact suffix as well as parts with one error. Those two properties confer a practical advantage over suffix filters introduced in [START_REF] Kärkkäinen | Faster filters for approximate string matching[END_REF], for example.

In the two following sections, we will provide two examples of algorithms based on 01 * 0 seeds, that show the practicability of these seeds.

Application to the pattern matching problem

Let us consider the problem of finding matches of the pattern P with at most k errors in the text T . We assume a small pattern (several dozens of letters) and a large text (millions or billions of letters) over a small alphabet (e.g. DNA) that is known in advance. This problem arises naturally in several applications in computational biology, such as identifying regulatory signals along the genome, predicting targets of non-protein coding small RNAs or analysing spacers in CRISPR for potential transfers from viruses or plasmids, to mention a few. More generally, introducing some errors would improve the sensitivity in the presence of sequencing errors or variants.

We devise an efficient filtration algorithm for this problem based on the 01 * 0 seeds. We will first justify our choice of using an FM-index for T . Then we will explain how seeds are searched in the FM-index, and how the property of the 01 * 0 seeds is used during the verification step to restrict the number of errors searched for.

Choice of index

We are in the situation where the text is known in advance, and we may have hundreds of short sequences to be queried in the text. This situation is particularly suitable for text indexes. Using a text index, it is possible to reduce time consumption at the expense of space consumption. Two main families of indexes are used: q-gram indexes and full-text indexes. The former allow to efficiently recover occurrences of a fixed-length word, while the latter allow to search for any pattern of any length. A third family of indexes consists of indexes specifically designed for approximate search [START_REF] Chávez | A Metric Index for Approximate String Matching[END_REF][START_REF] Maaß | Text indexing with errors[END_REF][START_REF] Chan | A linear size index for approximate pattern matching[END_REF][START_REF] Belazzougui | Improved space-time tradeoffs for approximate full-text indexing with one edit error[END_REF]. However, these indexes are not compressed indexes (i.e. whose space consumption is proportional to the empirical entropy of the text) and, to the best of our knowledge, no implementation of the proposed solutions exists.

Compressed full-text indexes appear to be the indexes of choice, because they limit space consumption and allow for patterns of arbitrary sizes. Among compressed indexes, FM-indexes [START_REF] Ferragina | Indexing compressed text[END_REF] have an optimal time complexity for counting the occurrences of a pattern, while pattern search is more complex and counting is more time consuming with LZ-indexes [START_REF] Ferragina | Compressed text indexes: From theory to practice[END_REF]. The FM-index accompanied by the Burrows-Wheeler transform has been successfully used in a variety of bioinformatics tools [START_REF] Petri | Efficient indexing algorithms for approximate pattern matching in text[END_REF][START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF][START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF]. We will show in the remainder of this section that 01 * 0 seeds are appropriate for a search in the FM-index.

Seed filtration

Given a pattern P , we enumerate all possible subseeds for the pattern. Each subseed for P is characterized by two parts P i and P j , 1 ≤ i < j ≤ k + 2, that occur exactly in the text. According to Lemma 1, all the intervening parts between P i and P j must be searched with exactly one error. We recall that in the FM-index, patterns are searched backwards, therefore, we first start by searching any part P , with 1 < ≤ k + 2, assuming it is P j . This is an exact search in the index. Then the parts preceding P are searched with at most one error (by backtracking as in BWA for instance [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF]). When a part is found exactly, we know that P i has been reached. Starting with P , we can have several parts that fulfil our requirements; on reaching different parts P i1 , . . . , P iq each of them matching exactly at different locations in the text. All the possible solutions are searched. If P cannot be found exactly or if a part cannot be found with at most one error, this P is skipped and we move on the next one. Therefore, at most we will have considered the (k+1)(k+2) 2 possible pairs (i, j). All those pairs are searched within the FM-index.

Example 2. Let us continue with Example 1, also shown in Figure 1: k = 3 and P = AACG TGAG GTAG GTTC CATG, which is partitioned into 5 parts of equal length. Assume that the text T is the concatenation of the three strings at distance 3 from P : T =AACGGAGGTAAGTTCTCATGAACGTAGGCAAGTTCCATGATCGTGACGTAGGGTCCATG.

-The algorithm first tries j = 5. P 5 = CATG is found with no error in the FMindex. So, it has some exact occurrences in the text. Therefore, we continue to go through the FM-index to extend P 5 to the left and find all possible values for i. We find i = 4 (P 4 occurs exactly), i = 3 (P 4 occurs with one error and P 3 exactly) and i = 1 (P 4 , P 3 and P 2 occur with one error and P 1 exactly). This gives three different seed instances, leading to three seed occurrences.

P 1 P 2 P 3 P 4 P 5 P 4 P 5 P 3 P 4 P 5 AACG TGAG GTAC GTTC-CATG GTTC CATG GTAG GTTC CATG | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | AACG -GAG GTAA GTTCT CATG AACGTAGGCAA GTTC CATG ATCGTGAC GTAG GGTC CATG
-With j = 4, GTTC occurs exactly in the FM-index, and corresponds to two occurrences in T . By extending P 4 to the left, we keep just one instance since the second one cannot be extended to P 3 = GTAG with at most one error.

P 1 P 2 P 3 P 4 P 4 AACG TGAG GTAC GTTC GTTC | | | | | | | | | | | | | | | | | | AACG -GAG GTAA GTTC TCATGAACGTAGGCAA GTTC CATGATCGTGACGTAGGGTCCATG
Note that in this particular case, the first occurrence of P in T is covered by two overlapping 01 * 0 subseeds, characterized by i = 1 and j = 5, and i = 1 and j = 4, respectively. This redundancy is solved with the extension and verification step, which is described in the next subsection.

-With j = 3, we have two occurrences of GTAG in the text. The first one cannot be extended to the left with P 2 = TGAG. As for the second occurrence, P 2 is found with one error, but P 1 = AACG does not exactly match. So, the occurrence is discarded.

P 3 P 2 P 3 GTAG TGAG GTAG | | | | | | | | | | | AACGGAGGTAAGTTCTCATGAAC GTAG GCAAGTTCCATGATCG TGAC GTAG GGTCCATG -With j = 2,
there is no exact occurrence of the part TGAG in the text.

At this point, all the seed instances occurring in the text are identified. We then proceed to the extension and verification step.

Elongation and verification

To perform the elongation of an instance of the seed, we first need to have a deeper look at the error distribution along the pattern. We know that the subseed instance has a Levenshtein distance of ji -1 with P i . . . P j , which makes ji -1 errors. The following lemma gives an upper bound for the number of errors in P 1 . . . P i-1 and in P j+1 . . . P k+2 . Lemma 4. Let U be a string of Σ * such that lev(P, U ) = y ≤ k. Then there exists a 01 * 0 subseed P i . . . P j such that the prefix P 1 . . . P i-1 contains exactly i -1 -(ky) errors and the suffix P j+1 . . . P k+2 contains exactly k + 2j errors.

Proof. The hypothesis of the Lemma can be rephrased as follows. Assume you have k + 2 containers numbered from 1 to k + 2, and y tokens with 0 ≤ y ≤ k. Let δ = ky. The proof is by recurrence on k.

We start by noting that for k = 1, the solutions are (1, 2) (token in container 3), (2, 3) (token in container 1) or (1, 3) (token in container 2), when y = 1. When y = 0 the unique solution is [START_REF] Wu | Fast text searching: allowing errors[END_REF][START_REF] Altschul | Basic local alignment search tool[END_REF].

Assume now that k > 1. We choose i and j, such that the pair (i, j) fulfils the conditions of Lemma 2, and i is the container furthest to the left satisfying this condition. Let η be such that i -1δ + η is the number of tokens present in containers 1 to i -1. There are three possibilities for η.

If η = 0: the tokens from containers 1 to i -1 sum up to i -1δ and the tokens from j + 1 to k + 2 sum up to k + 2j. So the pair of containers (i, j) fulfils the criteria of our lemma.

If η > 0: there are (i -1δ + η) + (ji -1) = jδ -2 + η tokens in containers 1 to j. This leaves y -(jδ -2 + η) tokens in containers j + 1 to k + 2. In other words we have kj + 2η tokens in the kj + 3 containers from j to k + 2 (as container j is empty). By recurrence, we know that there is a pair (j , j ), with j ≤ j < j ≤ k + 2, such that containers j and j are empty and there are jj -(η -1) tokens in containers j to j and there are k + 2j tokens in containers j + 1 to k + 2. Finally from container 1 to j we have (jδ -2 + η) + (jj -(η -1)) = j -1δ tokens. Therefore we found a pair (j , j ) that fulfils our conditions.

If η < 0: since container i is empty there are also i -1δ + η tokens until container i. By recurrence, we know that we can find a pair (i 1 , i 2 ) which fulfils Lemma 2. That contradicts our initial hypothesis, as i 1 < i 2 ≤ i.

As a consequence of this Lemma, seed instances are first extended to the left, to find P 1 . . . P i-1 with at most i -1 errors. To gain more efficiency, this extension is directly carried out in the FM-index. Indeed, the retrieval of the positions of occurrences is the most time consuming part in an FM-index (in O(log 1+ n) per occurrence [START_REF] Ferragina | Compressed representations of sequences and full-text indexes[END_REF]). Once this extension has been performed, all occurrences of the prefix P 1 . . . P j are retrieved. Then the extension to the right is performed in the text using a banded dynamic programming algorithm. The starting point of the extension is the ending position of the occurrence of P 1 . . . P j in the text. Following Lemma 4, P j+1 . . . P k+2 must be searched with at most k-j+2 errors in the text. Therefore, the bandwidth is 2×(k-j+2)+1 in the dynamic programming algorithm. Note that the extension to the right could also have been performed in the index using a bidirectional Burrows-Wheeler transform [START_REF] Schnattinger | Bidirectional search in a string with wavelet trees[END_REF][START_REF] Belazzougui | Versatile succinct representations of the bidirectional Burrows-Wheeler transform[END_REF]. That would, however, increase the memory footprint and provide only a moderate speed up, since many false positive seed instances have been removed at this step.

Implementation

This algorithm was implemented in a software called Bwolo. Bwolo is written in C++, with the help of SeqAn library and the FM-Index implemented within [START_REF] Döring | SeqAn an efficient, generic C++ library for sequence analysis[END_REF]. It can be downloaded from http://bioinfo.lifl.fr/olo. In this implementation, patterns are divided into parts whose length differs by at most one character.

Experimental results

In this section, we present some experimental results in order to measure the performance of our algorithm. We compare Bwolo to a selection of tools that were chosen for their complementarity. Widely utilized in bioinformatics, Exonerate is a generic tool for pairwise sequence alignment, which uses exact sparse dynamic programming to perform the search [START_REF] Slater | Automated generation of heuristics for biological sequence comparison[END_REF]. We use it as a standard for an on-line exact algorithm for our problem. RazerS3 is a read mapping program based on counting q-grams [START_REF] Weese | RazerS 3: Faster, fully sensitive read mapping[END_REF]. It performs the verification via an implementation of the improved Myers bit-vector algorithm proposed by Hyyr [START_REF] Hyyrö | A Bit-vector Algorithm for Computing Levenshtein and Damerau Edit Distances[END_REF]. RazerS3 works without a precomputed index for the text. Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF], like our tool, indexes the text with an FM-index and searches for exact contiguous seeds. It then uses backtracking for handling errors and dynamic programming to build the full alignment. Readaligner [START_REF] Mäkinen | Unified view of backward backtracking in short read mapping[END_REF] relies on suffix filters [START_REF] Kärkkäinen | Faster filters for approximate string matching[END_REF] searched within an FM-index too. Lastly, we used an in-house implementation for approximate search in an FM-index written with the SeqAn library. It is based on a breadth-first search method with no prior filtration. Unfortunately, we were not able to include hybrid methods described in [START_REF] Russo | Approximate string matching with compressed indexes[END_REF] in our benchmark, since the implementation is not available.

All these tools were configured to be fully sensitive and output all occurrences of the pattern. Exonerate was launched using the options --exhaustive -m affine:bestfit --bestn 1 --score -3 --gapextend -1 --gapopen -1 --showalignment yes --verbose 0, RazerS3 with --filter pigeonhole --percent-identity

[Id] --recognition-rate 100 -f such that [Id] = 100×(1-k m ), Bowtie2 with -a -L [Seeds] -i C,[Seeds],0 such that [Seeds] = m
k+1 and Readaligner with options --all -i 3. Moreover, for each tool the score system is based on the unit score, which computes the Levenshtein distance.

The tests were run on a single thread of a server equipped with two Intel(R) Xeon(R) CPU E5-2420 and 205GB of RAM. The CPU time and the memory consumption were measured using the GNU time command.

Randomly generated sequences

This first test uses independent and identically distributed sequences on the DNA alphabet. The size n of the sequences ranges from 10 4 to 10 9 . We also generated 100 patterns of length 20 at random and measured the computation time of each tool for k = 2 and k = 3. Results are shown in Figure 2.

Bwolo is the fastest tool for long sequences, from 10 6 nucleotides, apart for k = 2 where Readaligner is the fastest when the sequence is longer than 2 × 10 8 . As expected, the added-value of Bwolo is even more obvious when k = 3. The two tools with no filtration, Exonerate and the exact search in the FM-index, are slow. Bowtie2 operates slowly compared to all the other tools, especially with larger values of n. This confirms that Bowtie2's heuristics, which has been designed for long patterns (at least 50 nt) and few errors, is not well adapted to shorter patterns with higher error rate. Unfortunately, there is not yet a specialized tool for this type of problem. In our benchmark, Bowtie2 is obliged to use a seed with low filtering power that lets too many occurrences happen. This dramatically increases the verification effort due to the cost of retrieving text positions from the FM-index. Interestingly enough, RazerS3, which uses the same seed, functions well on this data. This is consistent with the fact that a linear method can, in certain conditions for large k and n, be more efficient than a method based on a text index [START_REF] Navarro | A Hybrid Indexing Method for Approximate String Matching[END_REF]. However, Bwolo is still five times faster than RazerS3 for sequences of length 10 9 . Indeed, the number of seed occurrences is an order of magnitude less with Bwolo, which offsets the additional time needed to query the FM-index in the verification step. With the configuration we used, Readaligner relies on suffix filters to align the sequences. Those filters have very good performances when the number of errors is not too high (k = 2 in our case). With suffix filters the most time-consuming case occurs when all parts contain exactly one error but the last one. With k = 3, the last part is 5nt long. Searching such a small part leads to many false positives.

For k = 2, we can observe that there are smaller differences in the CPU time between FM-index and Bwolo on larger texts. The former takes 18.4 s on the 1GB sequence while the latter takes 13.8 s. This small difference is actually misleading. Unserializing the index (which is the same in both cases) takes 12 s on that same sequence. Ignoring that step, which is constant whatever the number of sequences to search is, leads to a three-fold speedup using the 01 * 0 seeds compared to the breadth-first approach. With a higher error rate (k = 3) we have a seventy-five-fold speedup on the 1GB sequence. For the sake of comprehensiveness, we should mention that Readaligner takes 1s to unserialize its index and RazerS3 takes 8s to load the 1GB sequence. Taking that into account shows why Readaligner becomes quicker with long sequences and a small error rate: this is due to the loading of the index which is much more efficient than for RazerS3 and Bwolo.

All tools have a reasonable memory consumption, independent of the value of k, which grows linearly with the size of the text. For example, it is 13 MB for Readaligner, 27 MB for Bwolo, 99 MB for Bowtie2, 25 MB for RazerS3, and 31 MB for Exonerate for n = 10 7 . The memory consumption of Bwolo and Bowtie2 is dominated by the size of the FM-index. It is larger for Bowtie2 because it also deals with the inverted text and uses a different implementation. It is quite surprising that RazerS3 and Exonerate have a memory peak in the same order. It may be possible that they load both all the text and keep all results in memory.

Reads from the Human genome

The second test uses an external dataset made of short sequences, taken from [START_REF] Schbath | Mapping Reads on a Genomic Sequence: An Algorithmic Overview and a Practical Comparative Analysis[END_REF]. In this article, the H 3 dataset contains 10 millions of sequencing reads of length 40 that have been generated from the Human genome (assembly 37.1 from the NCBI, 25 chromosomes for 2.7 Gbp) with exactly three mismatches. Compared to the previous test, it allows us to evaluate the performance of the software with longer patterns, hence longer seeds. The maximum number k of errors is 3 (including indels, not only substitutions). We ran Bwolo, RazerS3, and Bowtie2 on the full set of reads (10 7 reads). Since we were not able to obtain results with both Bowtie2 and Readaligner on the full dataset within a reasonable amount of time, we also used a selection of 10,000 reads.

Application to the approximate dictionary matching problem

In this section, we consider the approximate dictionary matching problem: finding all approximate occurrences of a set P of string patterns in a query text Q. More specifically, we assume that P is a database of short patterns (about 20 letters each) that is known in advance and can be preprocessed. We want to query this database with a sequence Q of a few hundred letters. We will show how the 01 * 0 seeds can be used to solve this problem by indexing P.

Indexing pairs of exacts parts in 01 * 0 seeds

Let us assume that we have a collection P of patterns and that all patterns have the same length m. We will discuss later, in Subsection 4.4, the case where the patterns have different lengths. The query is a sequence Q whose length n is greater than m.

In accordance with the principle of the 01 * 0 seeds, every pattern P of P is split into k + 2 non overlapping parts, P 1 , . . . , P k+2 . For the sake of simplicity we will assume for now that k + 2 divides m and that all parts have the same length p. The first step of the algorithm is to search for all possible pairs of exact parts of all patterns of P that are present in the text Q. For this purpose, we work subseed by subseed and for each pair i and j, 1 ≤ i < j ≤ k + 2 corresponding to a given subseed, we will use an array A i,j , indexed by strings of length 2p, corresponding to pairs of parts. Given a string S of length 2p, A i,j [S] is the set of all patterns whose part i concatenated to part j equals S: A i,j [S] = {P ∈ P; P i P j = S} There are as many arrays as possible subseeds, which makes k(k + 1) arrays, each of size σ 2p . This is of course practical only for small values of p and σ, which is the case in our application.

Example 3. We show an example of construction of the A i,j arrays in Figure 3. The set P consists of three patterns of length 8: P 1 = ATACCACT, P 2 = TAACATCT, P 3 = ACCATTAT. The number of errors is k = 2. So each pattern is split in four parts of length 2. We have 6 arrays, A 1,2 , A 1,3 , A 1,4 , A 2,3 , A 2,4 and A 3,4 . We give the content of each of these 6 arrays. For example, A 1,2 has three elements. The first one is 3 indexed by ACCA, meaning that the pattern P 3 is such that P 3 1 = AC and P 3 2 = CA. For the array A 2,4 , we have A 2,4 [ACCT] = {1,2}, because P 1 and P 2 share the same second part, which is AC, and the same fourth part, which is CT.

Querying the A i,j arrays

The two exact parts of the 01 * 0 seeds are identified using the A i,j arrays. We still have to determine what is the range of distances between those two parts on Q to make it possible that P i and P j are the exact parts of a subseed of a 01 * 0 seed. For that, we will use the weakened form of Lemma 1, introduced in Section 2.2 for the bounded k + 2 pigeonhole seed. Let assume that for a given pattern P , the parts P i and P j match on the query sequence Q. There are ji -1 parts between the exact parts P i and P j in the pattern. Each part has length p, and should occur with exactly one error in Q. So the total number of AT AC CA CT TA AC AT CT AC CA TT AT

P 1 1 2 3 4 P 2 1 2 3 4 P 3 1 2 3 4 ACCA 3 ATAC 1 TAAC 2 A 1,2 ACTT 3 ATCA 1 TAAT 2 A 1,3 ACAT 3 ATCT 1 TACT 2 A 1,4 ACAT 2 ACCA 1 CATT 3 A 2,3 ACCT 1,2 CAAT 3 A 2,4
ATCT 2 CACT 1 TTAT 3 positions between P i and P j on Q should be comprised between (ij -1)(p -1) and (ij -1)(p + 1). Therefore at each position of the text Q, we must search for: where occ is the number of occurrences. Indeed, queries in A i,j arrays can be done in constant time (plus the time needed to retrieve the occurrences). This step is illustrated in the example below.

A 3,4 Q = A A T C A C A A T C T T C A A
1. two consecutive parts in A 1,2 , A 2,3 , . . . , A k+1,k+2 ,
Example 3 (continued). We now have a text to query: Q = AATCACAATCTTCAA. We show how the text is scanned in the bottom part of Figure 3. At a given point, the parts queried will be the ones in a frame. While maintaining the first part identical (AC), the second part queried in A 2,4 will be AT, TC and CT as we must account for potential deletion or insertion in the inner part. The case shown in the figure illustrates the case where an insertion is assumed.

Aligning candidates against the query

Once the candidates have been retrieved, they must be aligned against Q to check if there are at most k errors between the pattern and Q. Let assume that P is a pattern of P that has been selected by the preceding filtering step of the algorithm. By construction, there exist i and j, 1 ≤ i < j ≤ k + 2, q, 0 ≤ q ≤ n -1 and δ, (ji -1)(p -1) ≤ δ ≤ (ji -1)(p + 1), such that P i = Q[q . . . q + p -1] and P j = Q[q + δ . . . q + δ + p -1]. First we need to check whether the two parts P i and P j could be the exact parts of a 01 * 0 subseed. For that, the alignment of P i+1 . . . P j-1 against Q[q + p . . . q + δ -1] must verify that each part has exactly one error to the corresponding substring in Q. Then to extend this subseed, we must align the two remaining substrings P 1 . . . P i-1 and P j+1 . . . P k+2 to Q by dynamic programming. Similarly to what is done in Bwolo, we can bound the number of errors for each subtring and speed up the computation accordingly. When aligning P 1 . . . P i-1 against Q, we know by Lemma 4 that we can assume that there are at most i -1 errors.

Example 3 (continued). We continue with the example of Figure 3. Recall that at this step of the algorithm, the framed subwords correspond to parts P 1 2 and P 1 4 in P 1 , and to parts P 2 2 and P 2 4 in P 2 . The next step is to align AAT (the inner sequence between the two considered parts in Q) against P 1 3 = CA and P 2 3 = AT with exactly one error. As lev(CA, AAT) > 1, P 1 is not a candidate anymore. Regarding P 2 , since lev(AT, AAT) = 1, we have found an instance of a subseed. It means that P 2 is still a valid candidate, that must be fully aligned against Q. As we are considering parts 2 and 4, we must align the prefix P 2 1 on Q, before position 5, with exactly one error (if part 1 had no error we would already have a hit).

Dealing with variable-length database sequences

For the sake of simplicity, we started by assuming that patterns in P all have the same length m and that k + 2 divides this length. We will now explain how to relax those assumptions.

First we will consider the case where a pattern P of P is strictly longer than m. The first k + 1 parts of P will all have the same length p. The length of the last part will be p = |P | -(k + 1)p. Therefore the first k + 1 parts of P are treated in the same way as the parts of the other sequences. This cannot be true for the last part whose size can vary. This part of P is still indexed on its first p letters. However the A i,k+2 arrays (with 1 ≤ i ≤ k + 2) have to store slightly more information compared to the other arrays. In addition to the indices, the trailing pp letters of the last parts are also stored:

A i,k+2 [S] = {(P, w); P i P k+2 [1 . . . p] = S and w = P k+2 [p + 1 . . . p ]}
Therefore when querying any array A i,k+2 with some subtring S = Q[q . . . q + p -1]Q[q . . . q + p -1], all the entries (P , w) ∈ A i,k+2 [S] are retrieved, such that w = Q[q + p . . . q + p + |w| -1]. Rather than storing entries of A j,k+2 as plain lists, they are stored in a trie. In this trie, a path labelled w from the root to a leaf, gives all the patterns stored in A j,k+2 that have w as trailing letters.

Second we consider the case where a pattern P in P is shorter than m. We will assume that only the last part is shorter than the other ones. Let p = |P k+2 |, for each part i, 1 ≤ i ≤ k + 1, we will have P ∈ A i,k+2 [P i P k+2 w], for all w ∈ Σ p-p . If other parts also have to be smaller, we process them similarly.

Implementation

This algorithm for the multi-pattern matching problem has been implemented in C++ in a software called Piccolo, which is available at http:// bioinfo.lifl.fr/olo/. Piccolo uses the same partition as Bwolo. For the A i,j arrays, the index data structure is simply a lookup table, like in Blastn.

Experimental results

We measure the performance of Piccolo on a case study from computational biology: the classification of plant miRNAs (microRNAs) in families. Plant miRNAs are short non-coding RNAs that are 21-24 nucleotides long, and that are produced from larger precursor sequences that can be hundreds of nucleotides long. Many miRNA genes are conserved across plant species, and form gene families. So an important step for the miRNA annotation is to compare a precursor sequence to a known miRNAs, in order to identify homologs. For that, the protocol is to align the sequence of the precursor to a database of known miRNAs, such as miRBase [START_REF] Kozomara | miRBase: annotating high confidence mi-croRNAs using deep sequencing data[END_REF]. What is found commonly in the literature is to require a full-length alignment of a miRNA with up to 3 errors [START_REF] Chorostecki | Identification of new microRNA-regulated genes by conserved targeting in plant species[END_REF].

To our knowledge, there exist no dedicated tools for such a task, either with a precomputed index or not. One possibility is to perform dynamic programming on each miRNA sequence of the database against the precursor sequence. This is slow in practice, and one has to resort to heuristic methods, in particular Blastn. In this case Blastn provides quick results at the cost of some false positive and false negative results.

We will therefore test our program against Exonerate (as in Section 3) and Blastn. The data was retrieved from miRBase 21. It consists of 124 miRNA precursors from Amborella trichopoda (minimum length: 80, maximum length: 2354, average length: 284, median length: 228) and 4,713 plant miRNAs (excluding those from Amborella trichopoda, minimum length: 17, maximum length: 26, average length: 22, median length: 21). The programs were launched on one thread of an Intel i7-4600 dual core processor. Exonerate was run with the same options as in Section 3. For fair comparison we used a slightly modified version that searches matches on the forward strand only. Blastn was launched in its standalone version but with the smallest seed available in the online version (as most users usually do), which spans seven characters. Blastn does not allow arbitrary costs for the gaps; we chose the costs that were closest to those for Exonerate. The options were -strand plus -outfmt 6 -task blastn-short -min raw gapped score 53 -penalty -5 -reward 4 -gapopen 5 -gapextend 5.

The computation times, including the indexation time for Blastn and Piccolo, are presented in Table 3. Unsurprisingly the dynamic programming approach of Exonerate is much slower. It takes more than 400 times longer than Blastn and Piccolo. The very well optimized Blastn is not the fastest, as Piccolo is slightly faster. Moreover, Blastn has higher memory footprint. This shows that even a naive implementation of the A i,j arrays based on a lookup table yields good performances on the space consumption side. We deliberately time memory alignments false negatives false positives opted for this index data structure to enhance the running time. In presence of a larger database of patterns, the algorithm could also be implemented with hashtables or a compressed full-text index. Exonerate and Piccolo, which are exact methods, both find 1,428 alignments. Since Blastn is a heuristic, it has false positive and false negative results: Blastn misses 146 alignments (almost 10%), and outputs 225 local alignments that contain more than three errors when extended to the full length miRNA. Such false results can bias the type of miRNAs found and therefore could hinder the understanding in their evolution. In conclusion, Piccolo is as fast as Blastn, but with the advantage of being fully sensitive.

Application to the k-mismatch problem

In this last section, we make a short digression and study the application of the 01 * 0 seeds to the k-mismatch problem: given a pattern P , a threshold k, find all strings U whose Hamming distance between P and U is less than k. 01 * 0 seeds are not specifically designed to solve this problem. The reason why we perform this analysis is because it allows us to compare the filtering power of 01 * 0 seeds to spaced seeds. A spaced seed consists of an -mer in which some predefined positions can match any character symbol in the alphabet. Exact match positions are usually denoted #, and positions which are not required to match -. Spaced seeds are known to offer a good sensitivity-selectivity trade off when dealing with mismatches [START_REF] Brown | A tutorial of recent developments in the seeding of local alignment[END_REF][START_REF] Keich | On spaced seeds for similarity search[END_REF]. In particular, spaced seeds are better than contiguous seeds. Note also that the k-mismatch problem has attracted attention very recently in [START_REF] Kucherov | Approximate string matching using a bidirectional index[END_REF], and that our 01 * 0 seeds when restrained to substitutions could be viewed as a special case of the framework introduced in their paper.

Given a pattern P of length m partitioned in k + 2 parts, P 1 , . . . , P k+2 , we define S(P 1 , . . . , P k+2 ) to be the set of words of Σ m that are recognized by the 01 * 0 seed based on this partition. Our goal is to compute the cardinality of S(P 1 , . . . , P k+2 ). Since we have only mismatches, this value depends only on the distribution of p 1 , . . . , p k+2 , the lengths of P 1 , . . . , P k+2 , and on the size σ of the alphabet. of 01 ⇤ 0 seeds to spaced seeds. Spaced seeds consist of a `-mer in which some predefined positions can match any character symbol in the alphabet. Exact match positions are usually denoted #, and positions which are not required to match -. Spaced seeds are known to o↵er a good sensitivity-selectivity trade o↵ when dealing with mismatches [START_REF] Brown | A tutorial of recent developments in the seeding of local alignment[END_REF][START_REF] Keich | On spaced seeds for similarity search[END_REF]. In particular, spaced seeds are better than contiguous seeds. Note also that the k-mismatch problem has attracted attention very recently in [START_REF] Kozomara | miRBase: annotating high confidence mi-croRNAs using deep sequencing data[END_REF], and that our 01 ⇤ 0 seeds when restrained to substitutions could be viewed as a special case of the framework introduced in their paper. Given a pattern P of length m partitioned in k + 2 parts, P 1 , . . . , P k+2 , we define S(P 1 , . . . , P k+2 ) to be the set of words of ⌃ m that are recognized by the 01 ⇤ 0 seed based on this partition. Our goal is to compute the cardinal of S(P 1 , . . . , P k+2 ). Since we have only mismatches, this value depends only of the distribution of p 1 , . . . , p k+2 , the lengths of P 1 , . . . , P k+2 , and on the size of the alphabet.

Lemma 5. The cardinal of S(P 1 , . . . , P k+2 ) is k+2 where k+2 is defined by the linearly recursive sequence.

0 @ ↵ 0 0 0 1 A = 0 @ 1 0 0 1 A 0 @ ↵ i i i 1 A = 0 @ p i 1 p i ( 1)p i 1 0 1 ( 1)p i 0 0 1 p i 1 A 0 @ ↵ i 1 i 1 i 1 1 A Proof. We convert the problem in a new alphabet, {0, 1, 2}. For each i, 1  i  k + 2, define H i : ⌃ p i ! {0, 1, 2} by
H i (U ) = 0, if ham(P i , U) = 0 H i (U ) = 1, if ham(P i , U) = 1 H i (U ) = 2, otherwise (ham(P i , U) 2)
where ham denotes the Hamming distance. In this context, the set S(P

  α 0 β 0 δ 0   =   1 0 0     α i β i δ i   =   σ pi -1 σ pi -(σ -1)p i -1 0 1 (σ -1)p i 0 0 1 σ pi     α i-1 β i-1 δ i-1   Proof. We convert the problem in a new alphabet, {0, 1, 2}. For each i, 1 ≤ i ≤ k + 2, define H i : Σ pi → {0, 1, 2} by
H i (U ) = 0, if ham(P i , U ) = 0 H i (U ) = 1, if ham(P i , U ) = 1 H i (U ) = 2, otherwise (ham(P i , U ) ≥ 2)
where ham denotes the Hamming distance. In this context, the set S(P

(U 1 ) • • • H i (U i ).
From the DFA, we infer the following relations. where for each i, 1 ≤ i ≤ k + 2, for each x in {0, 1, 2}, H i (x) denotes the number of strings U of Σ pi that are mapped to x by H i . In other words, this is the cardinality of H -1 i (x). From the definition of H i , it is straightfoward to see that

  α 0 β 0 δ 0   =   1 0 0     α i β i δ i   =   H i (1) + H i (2) H i (2) 0 H i (0) H i (1) 0 0 H i (0) H i (0) + H i (1) + H i (2)     α i-1 β i-1 δ i-1   m = 20 k + 1 k + 2 spaced-
m = 20 k = 1 k = 2 k = 3 k = 4 k + 1
H i (0) = 1 H i (1) = (σ -1)p i H i (2) = σ pi -(σ -1)p i -1
This gives the expected result.

Lemma 5 allows us to characterize exactly the number of words caught by the 01 * 0 seed. In Table 4, we report a comprehensive comparison with spaced seeds and k + 1, k + 2 pigeonhole seeds for the special case of a pattern of length m = 20 on the DNA alphabet, and values of k ranging from 1 to 4. seeds were designed with Iedera [START_REF] Kucherov | A unifying framework for seed sensitivity and its application to subset seeds[END_REF][START_REF] Kucherov | Subset seed automaton[END_REF]: in each case, we selected the lossless spaced seed with maximum weight (the number of #), and then maximum length (the number of # plus the number of -). The two k + 2 pigeonhole seeds considered so far are identical here, since only mismatches are allowed: Two exact parts P i and P j should be separated by exactly p i+1 + • • • + p j-1 positions. In the bottom table, we report the number of words of length m = 20 recognized by the corresponding seed. For the k + 1 pigeonhole seed, we require that one part occurs exactly in the string, at the same position. For the k + 2 pigeonhole seed, we require that two parts occur exactly in the string, at the same positions. For the spaced seed, we require that there exists at least one seed shared by the pattern and the string. Finally, for the 01 * 0 seed, we use Lemma 5. This table shows that the 01 * 0 seed is always the best seed, and that it improves the filtering power of spaced seeds by a factor of at least 20 as soon as there are at least 2 errors in the pattern.

Conclusion

We have introduced a new type of seeds, named 01 * 0 seeds, that are especially well-suited to deal with patterns containing a high rate of errors. We have shown that these seeds are easy to use with two complementary applications: the pattern matching problem, where the seeds are efficiently searched in a compressed full-text index, and the multi-pattern matching problem where the seeds are indexed in a classical lookup table and the text is processed on-line. Each of these two algorithms were implemented. The programs are open source and are available at http://bioinfo.lifl.fr/olo. We also made some theoretical comparisons to contiguous and spaced seeds and we showed that 01 * 0 seeds have a better filtering power than their counterparts under the k-mismatch problem.

We believe that 01 * 0 seeds constitute a promising alternative to existing approaches. They should prove useful for a variety of applications where the sensitivity of the search is critical. We have illustrated this through a selection of examples coming from computational biology. Another application in this field concerns the processing of third generation sequencing data, such as produced by PacBio or Oxford Nanopore sequencers [START_REF] Thompson | The properties and applications of single-molecule DNA sequencing[END_REF]. Those emerging technologies are known to deliver longer sequencing reads at the price of a higher error rate. So exploiting the power of these sequencers requires the development of new algorithms for the assembly problem as well as for the mapping problem that are able to accomodate error-prone sequences, which is a key point of 01 * 0 seeds.

The 01 * 0 seeds also raise new computational problems. The first one is a formal analysis of the selectivity of those seeds. In this article, we have only provided an empirical study on randomly generated sequences. Another question of interest is the indexing strategy of 01 * 0 seeds, that could still be improved as we must take into account the potential shifts introduced by insertions or deletions in the 1 * parts of the 01 * 0 seeds. A future research would be to devise a more efficient indexing strategy that could help in lowering the query complexity.

Figure 1 :

 1 Figure 1: Application of Lemma 1 for sequences of Example 1.The pattern AACGTGAGGTAGGTTCCATG is on the top of each alignment. The two parts with no error (written P i and P j in the Lemma) are highlighted in dark grey, and the parts with 1 error surrounded by P i and P j are highlighted in light grey.

Figure 2 :

 2 Figure 2: Running time for 100 randomly sequences. Both axes are in logarithmic scale. Bwolo is our algorithm. FM-index refers to the breadth-first search implementation in FM-index.

15 A 2 , 4 [Figure 3 :

 15243 Figure 3: Example of indexing and querying strategy on a set P of three patterns of length 8: P 1 = ATACCACT, P 2 = TAACATCT, P 3 = ACCATTAT. The number of errors is k = 2. Each pattern is split in four parts of length 2.

2 .

 2 two parts separated by p -1, p or p + 1 positions in A 1,3 , . . . , A k,k+2 , 3. . . . 4. two parts separated by k(p -1), . . . , k(p + 1) positions in A 1,k+2 . So, for each d in [0..k] and for each i ∈ [1..k-d+1], the array A i,i+d+1 needs 2d+1 queries. Since there are k + 1d such arrays for a given value of d, this gives a total of k d=0 (2d + 1)(k + 1d) = O(k 3 ) queries. Finally identifying candidates in P which have compatible pairs with Q can be done in time O(|Q|k 3 + occ),

Figure 5 :

 5 Figure 5: Deterministic Finite Automaton for the regular expression {0, 1, 2} ⇤ 01 ⇤ 0{0, 1, 2} ⇤ . ↵ is the initial state and is the unique final state.

  3 #----#----#----# 4,4,3,3,3,3

Table

					bounded k + 2	lossy
			pigeonhole pigeonhole	pigeonhole	seed
	average	6,598	385,651	86,145	44,747 36,207
	std. dev.	106	4,221	2,078	967	419

Table 2

 2

		index construction	10, 000 reads	10 7 reads
		time	mem.	time	mem.	time	mem.
		(min.)	(GB.)	(min.)	(GB.)	(min.)	(GB.)
	Bwolo	127	9.6	1.6	6.5	925	9.1
	RazerS3	0	0	8.4	6.5	7, 790	152
	Readaligner	451	5.5	41	3.6	NA	NA
	Bowtie2	176	5.4	2, 603	8.3	NA	NA

shows the results. As in the previous test, Bwolo achieves the best performances. However, the difference between Bwolo, Readaligner and RazerS3 is even more striking than in the previous test. This is due to the time needed to load the index. It was negligible on this dataset, but it constituted an important part of the search time with a much smaller dataset as in the previous test.

Table 2 :

 2 Running time and space consumption on the Human genome benchmark. NA: non available.

Table 3 :

 3 Computation time and memory usage for Exonerate, Blastn and Piccolo for a dataset of 4,713 plant miRNAs searched among each of 124 miRNA precursors. The column alignments gives the total number of alignments found by the tool and satisfying the specification of the problem: up to 3 errors. For Blastn, we also report the number of false negatives (alignments that have not been found) and false positives (alignments that do not satisfy the specification).

  1 , . . . , P k+2 ) is exactly the set of words U of ⌃ m such that H 1 (U 1 ) . . . H k+2 (U k+2 ) is accepted by the regular expression {0, 1, 2} ⇤ 01 ⇤ 0{0, 1, 2} ⇤ , where U 1 , . . . , U k+2 is the unique partition of U in parts of length p 1 , . . . , p k+2 respectively. The DFA 18 Figure 4: Deterministic Finite Automaton for the regular expression {0, 1, 2} * 01 * 0{0, 1, 2} * . α is the initial state δ is the unique final Lemma 5. The cardinality of S(P 1 , . . . , P k+2 ) is δ k+2 where δ k+2 is defined by the linearly recursive sequence.

  p1+

  ••+pi , such that the DFA starting at the initial state ends in the state α (respectively β, δ) after reading the string H 1

Table 4 :

 4 The top table displays the form of seed used for each value of k. For the k + 1 pigeonhole, k + 2 pigeonhole and 01 * 0 seeds, each line gives the length of each non-overlapping part. For the space seed, # are for fixed positions andare for variable positions. In the bottom table, we report the number of words of length m = 20 recognized by the corresponding seed.

Part of this work has already been published as an extended abstract in the proceedings of the conference Iwoca

[START_REF] Vroland | Lossless seeds for searching short patterns with high error rates[END_REF].


Solène Comby-Dassonneville

Fabien Volpi
email: fabien.volpi@simap.grenoble-inp.fr

Guillaume Parry

Didier Pellerin

Marc Verdier

Resistive-nanoindentation: contact area monitoring by real-time electrical contact resistance measurement

come

Introduction

In the past three decades, efforts have been made to couple instrumented indentation with resistive measurements [START_REF] Nili | In situ nanoindentation: Probing nanoscale multifunctionality[END_REF]. This development was driven by several motivations, such as the local monitoring of phase transformation [START_REF] Clarke | Amorphization and Conductivity of Silicon and Germanium Induced by Indentation[END_REF][START_REF] Pharr | Electrical resistance of metallic contacts on silicon and germanium during indentation[END_REF][START_REF] Mann | Size-dependent phase transformation during point loading of silicon[END_REF][START_REF] Mann | Contact resistance and phase transformation during nanoindentation of silicon[END_REF][START_REF] Bradby | In situ electrical characterisation of phase transformation in Si during indentation[END_REF][START_REF] Ruffell | An in situ electrical measurement technique via a conducting diamond tip for nanoindentation in silicon[END_REF][START_REF] Ruffell | Identification of nanoindentation-induced phase changes in silicon by in situ electrical characterization[END_REF], the study of native oxide fracture [START_REF] Pethica | Contact of characterised metal surface at very low loads: deformation and adhesion[END_REF][START_REF] Stauffer | Plastic response of the native oxide on Cr and Al thin films from in situ conductive nanoindentation[END_REF][START_REF] Nguyen | Electrical contact resistance for monitoring nanoindentationinduced delamination[END_REF], the investigation of MEMS operation at small scales [START_REF] Duvivier | Time dependence investigation of the electrical resistance of Au/Au thin film micro contacts[END_REF][START_REF] Arrazat | Discrete analysis of gold surface asperities deformation under spherical nano-indentation towards electrical contact resistance calculation[END_REF] and the contact area computation during nanoindentation test [START_REF] Fang | Continuous electrical in situ contact area measurement during instrumented indentation[END_REF][START_REF] Sprouster | Quantitative electromechanical characterization of materials using conductive ceramic tips[END_REF]. The later point is of particular interest for the determination of the contact area which is essential to compute both the sample Young's modulus and the hardness from nanoindentation tests:


The sample Young's modulus is determined using the Sneddon's relation [START_REF] Sneddon | The relation between load and penetration in the axisymmetric boussinesq problem for a punch of arbitrary profile[END_REF]:

(1)

with the contact stiffness (measured continuously during nanoindentation test) and the reduced modulus, expressed as:

(2) with , , , the tip and sample Poisson's ratios and moduli, respectively.



The sample hardness is computed using Eq. ( 3), with the load on sample.

(

) 3
The contact area (defined as the projected area of the contact interface between the tip and the sample) is however complex to monitor, even for the simplest case of homogeneous semi-infinite specimens. This is because contact area depends both on the tip geometry and on the contact depth ( being the length of the tip effectively in contact with the sample (Supplementary Fig. S1)). The tip geometry is either determined by AFM imaging or by using a calibration sample with well-known mechanical properties [START_REF] Oliver | Measurements of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF]. As shown in Supplementary Fig. S1, the contact depth depends on sample rheology: the sample can either sink-in or pile-up around the tip during nanoindentation. The contact depth then strongly differs from the total penetration depth (which is the depth reached by the tip from the initial specimen surface).

The standard methods widely used to extract contact area require analytical models based on material rheology assumptions, which use penetration depth and contact stiffness values to compute the contact depth [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF][START_REF] Loubet | Mechanical Properties and Deformation Behavior of Materials Having Ultra-Fine Microstructures : Nanoexperiments with a surface force apparatus[END_REF]. On the contrary, the direct monitoring of contact area by resistive-nanoindentation brings supplementary inputs to the quantitative analysis of indentation tests. This is expected to be an experimental alternative to standard analytical models.

A first attempt to extract mechanical contact area from electrical measurements was reported by Fang et al [START_REF] Fang | Continuous electrical in situ contact area measurement during instrumented indentation[END_REF] for micro-indentation. The authors successfully related non-linear current-bias curves performed during micro-indentation to real-time contact area. However, this approach is inadequate for nanoindentation. More promising results were obtained by Sprouster et al [START_REF] Sprouster | Quantitative electromechanical characterization of materials using conductive ceramic tips[END_REF], with resistive nanoindentation performed on a gold bulk sample with a vanadium carbide tip. The measured resistance during indentation was successfully fitted using an analytical model based on Maxwell expression for the spreading resistance (Eq. ( 4), with the material resistivity and the contact radius). With this method, the authors managed to compute their specimen resistivity but the ability to monitor the contact area during nanoindentation was not discussed.

(4)

In the present work, a complete method is proposed to compute contact area using resistivenanoindentation of noble metals. First, the different contributions to the measured resistance are discussed, then the computation procedure is detailed and applied step-by-step: (1) measurement of the tip shape, [START_REF] Clarke | Amorphization and Conductivity of Silicon and Germanium Induced by Indentation[END_REF] set-up calibration on a model specimen (gold bulk single crystal) and (3) extraction of the contact area on a specimen with complex rheology (gold film on sapphire substrate), and finally, the data generated by this method are compared to experimental AFM measures of contact areas.

Experimental details

Nanoindentation testing is performed with a home-made device, which includes a commercial actuator from Nanomechanics Inc. Resistance measurements are conducted with a ResiScope apparatus from Scientec [START_REF] Houzé | Imaging the local electrical properties of metal surfaces by atomic force microscopy with conducting probes[END_REF]. Nanoindentation tests are performed with a Berkovich boron-doped diamond (BDD) tip with resistivity in the range [0.2-2] Ω.cm. Resistive-nanoindentation tests were performed on a gold bulk single crystal (100) for the calibration step, and the whole method was validated on a 200 nm-thick polycrystalline gold film deposited on a sapphire substrate. Gold samples were chosen as model specimens, due to the absence of native oxide. Nonetheless this method can be applied to any oxide-free highly-conductive material. Indentation of a composite geometry such as an elasto-plastic film (gold) on an elastic substrate (sapphire) is a model case of complex rheology which is depth-dependent and where no analytical modeling exists. AFM imaging of residual imprints is performed after each test for validation purpose. Post-mortem contact areas are computed using Gwyddion software [START_REF]Gwyddion software[END_REF]: a map of local plane inclination angles is first extracted from the topological map of these imprints [START_REF] Charleux | A method for measuring the contact area in instrumented indentation testing by tip scanning probe microscopy imaging[END_REF]. Then manually-adjusted edge detection is performed on the obtained map (standard deviation on contact area lower than 3%). AFM images are performed after the final elastic unloading (viscous relaxation is negligible at room temperature): this leads to a small modification of the indentation imprint, as compared to the contact area under loading. According to Charleux et al. [START_REF] Charleux | A method for measuring the contact area in instrumented indentation testing by tip scanning probe microscopy imaging[END_REF], the induced error on contact area is lower than 10%.

Results and discussion

Resistive-nanoindentation curves

The resistance-versus-depth curves obtained during the indentation of the gold bulk sample on different spots are given in Fig. 1(a). This figure shows the high reproducibility of these experiments. As expected, the resistance continuously decreases as the mechanical contact area increases [START_REF] Sprouster | Quantitative electromechanical characterization of materials using conductive ceramic tips[END_REF]. <<color online>>

Identification of resistance contributions

Before starting the procedure for the contact area computation, the different contributions to the measured resistance have to be identified. As long as the electronic transport between the BDD tip and the gold specimen remains in diffusive regime, the measured resistance is the sum of the following contributions (from top to bottom in Fig. 2):

i.

The electrical resistance of the indentation tip ( );

ii.

An interface resistance ( ), coming from a possible interfacial layer;

iii.

The spreading resistance within the sample ( ), due to the current line constriction at the contact; iv.

The resistance from the gold sample ( );

v.

An additional series resistance from the experimental set-up ( ).

FIG. 2. Contributions to the measured resistance. Current lines are schematized by red-dotted lines.

The two resistance contributions coming from the gold sample ( and ) depend linearly on its resistivity, similarly to the tip resistance ( ) that depends linearly on the BDD tip resistivity. Since gold resistivity is at least five orders of magnitude lower than the BDD tip resistivity, and can be safely neglected. can also be neglected, due to the absence of oxide or of other contaminant layer at the tip-to-sample interface. Consequently, the measured resistance is simply the sum of the tip resistance and the additional constant series resistance:

(5)

The expression for can be adapted from a model proposed by Schneegans et al for the computation of a broken AFM tip [START_REF] Schneegans | Capacitance measurements on small parallel plate capacitors using nanoscale impedance microscopy[END_REF]: the broken apex of the AFM tip is, in our case, the part of the indentation tip which has penetrated the sample (Supplementary Fig. S2(a)). This consideration leads to Eq. 6 (with the tip resistivity and a constant related to the tip geometry):

However the experimental determination of is affected by the unavoidable tip defect (rounded apex of height (Supplementary Fig. S2(b)), thus Eq. 6 turns into Eq. 7:

(7)

This expression relies on two strong assumptions: 1/ the current distribution within the tip is homogeneous and parallel to the tip axis and 2/ the part of the tip inside the sample affects only linearly the tip resistance.

Regarding the first assumption, it has been shown in the literature that the error arising from the use of a homogeneous and parallel distribution of current lines within the tip is less than 10% [START_REF] Romano ; R | The conical resistor conundrum: a potential solution[END_REF]. The second assumption is more questionable. Since the BDD tip is five orders of magnitude more resistive than gold, current lines are expected to be localized at the periphery of the contact instead of being homogeneously distributed (Supplementary Fig. S3). This point has been confirmed by numerical computation (not shown). In literature, Nakamura et al [START_REF] Nakamura | Constriction resistance of conducting spots by the boundary element method[END_REF] have shown that the electrical resistance of ring contacts is equivalent to the resistance of the corresponding full-surface contact divided by a shape factor . This shape factor depends on the ratio between the ring thickness and the contact radius. For now, we first assume that this shape factor is constant during indentation (this assumption will be confirmed later).

Finally, the measured resistance is given by Eq. 8, where A and B indicate two constants that depend only on the experimental set-up (series resistance, tip geometry and resistivity) but not on the specimen. These constants have to be determined experimentally during the calibration step (below), and will be further used for the effective nanoindentation measurements.

(8)

Procedure for contact area computation

Now that these preliminary considerations have been stated, the step-by-step procedure for contact area computation can be described:

1. First the BDD tip geometry is determined. This step generates the tip "shape function" that relates the projected contact area to the contact depth.

2. An electrical calibration is carried out in order to determine the A and B constants (Eq. 8). As A and B depend only on the experimental set-up, a one-to-one analytical correspondence can be established between electrical resistance and contact depth (independently of the specimen).

3. The characterisation of a conductive oxide-free specimen can thus be performed. Using the shape function (step 1), the contact area is finally determined from the contact depth values (step 2) for this specimen to be characterized.

Step 1. The determination of the tip geometry allows us to relate the total penetration depth to the corresponding tip cross section . For that purpose, two complementary methods have been used: direct AFM imaging of the tip, and the calibration method described by Oliver and Pharr (where a model specimen with known isotropic elastic properties, like fused SiO 2 , is measured) [START_REF] Oliver | Measurements of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF]. The obtained shape function can be simplified, leading to Eq. 9. The error due to this simplification is less than 5% for penetration depth larger than 45 nm.

(

) 9
Step 2. In this work, electrical calibration has been performed on a gold bulk sample. Several resistive nanoindentation tests were run, up to different contact depths (circa 100, 300, 600 and 800 nm). Postmortem AFM imaging shows that residual imprints exhibit sink-in, probably due to polishing-induced surface hardening (Fig. 3(a)). Due to sink-in evidence, the contact depth can be reliably estimated from

Oliver and Pharr model [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF]. In order to validate the use of this model, the contact area was also computed from contact depth values, thanks to the tip "shape function" determined at step 1. The obtained contact area values were then compared to post-mortem AFM imaging, as shown in Fig. 3 At this point, it is possible to plot the measured resistance versus (Fig. 3(b)) so that Eq. 8 can be fitted. The linear correlation (coefficient of 0.99) validates the proposed analytical expression of the measured resistance (Eq. 8). This linearity also validates the assumption of a constant shape factor during indentation. A and B constants are then obtained by simple linear fitting. Using these fitted values, both experimental and analytical resistance versus depth curves can be plotted (Fig. 3(c)). Remarkably, a very good agreement is found between the experimental results and the analytical model for contact depths larger than 20 nm. Below these depths, the measured resistance is much higher than the model prediction. This deviation is most likely due to a change in electronic transport regimes from diffusive to ballistic since in ballistic regime, the electrical resistance evolves as [START_REF] Sharvin | A possible method for studying Fermi surfaces[END_REF]. This trend is in good agreement with that observed experimentally at small depths (Fig. 3(c)).

Alternatively, this deviation could also be explained by a multi-contact mode [START_REF] Greenwood | Constriction resistance and the real area of contact[END_REF], due to a potential thin interfacial layer at the contact.

Step 3. Once the entire analytical expression is fully set from the calibration sample, the specimen to be characterized can be processed: contact depth can then be directly determined from the resistance measurements. In the present work, this has been applied to a 200 nm polycrystalline gold thin film on sapphire substrate. This sample has a complex rheology, as the film thickness acts as an additional length scale. The deformation field during nanoindentation is no more self-similar as in the bulk case. The standard practice to analyze such a soft-film-on-stiff-substrate structure is usually to determine the contact area by numerically solving a Fredholm integral [START_REF] Perriot | Elastic contact to a coated half-space -Effective elastic modulus and real penetration[END_REF] (or using finite elements simulations) coupled with at least one post-mortem measurement of the remaining imprints at a given depth (by SEM or AFM).

Experimental validation

Resistive-nanoindentation tests with various final penetration depths were performed on this gold thin film. For each test, the contact depths were extracted continuously from resistance measurements (using Eq. 8) and the corresponding contact areas were computed from the tip shape function (using Eq. 9). After each test, post-mortem AFM measurements were conducted for validation. Fig. 4 reports the contact areas extracted from this procedure as well as post-mortem imaging data. An excellent agreement is observed, thus validating the whole methodology. It can be seen that the deformed film at large relative penetration depth (75%) presents a large upheaval. At this depth the contact area extracted from the Oliver-Pharr's method [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF] (based on sink-in assumption) would be ~50% lower. The present procedure should apply to any specimen that satisfies the two following conditions: (1) the specimen is much more conductive than the indentation tip and (2) the specimen is oxide-free. The first condition is easily met for bulk metals, while the relevant magnitude is the sheet resistance for thin films. In the case of oxidized specimens,

Conclusion

This paper reports an efficient method to extract the contact area during nanoindentation from resistance measurements. This method relies on three steps (tip shape measurement, set-up calibration, application to the specimen under test), and applies to oxide-free conductive materials. In the presence of an oxide layer, an additional non-linear contact resistance is expected, thus requiring the adaptation of the approach presented in this paper. Finally, the method is applied to nanoindentation tests performed on a sample with a complex rheology. The excellent agreement with the experimental area determination validates the proposed method. To our knowledge, this is the first time in literature that the contact area is extracted from electrical measurements during nanoindentation. It also paves the way for new calibration procedure of both tip geometry and frame stiffness.

Figure captions

Fig. 1 (FIG. 1 .

 11 FIG. 1. Resistance evolution during indentation of gold bulk single crystal. (a) Different tests

  FIG. 3. Results from resistive-nanoindentation on bulk gold sample. (a) Post-mortem AFM measures

FIG. 4 .

 4 FIG. 4. Analysis of the resistive-nanoindentation tests on the 200 nm gold thin film. (a) Post-mortem

FIG. 1 .FIG. 2 .FIG. 3 .FIG. 4 .

 1234 FIG. 1. Resistance evolution during indentation of gold bulk single crystal. (a) Different tests performed at

Acknowledgments

This work has been performed with the financial support of the Centre of Excellence of Multifunctional Architectured Materials "CEMAM" n° ANR-10-LABX-44-01. The CEMAM program is funded by the French Agence Nationale de la Recherche (ANR).

The authors thank the technical team of SIMaP lab (Bruno Mallery, Stéphane Massucci and Nadine Vidal) for its support, as well as members of CSI/Scientec company (Les Ulis, France) for their support for the functionalisation of the nanoindentation head: Louis Pacheco, Alan Lecoguiec and Sylvain Poulet.

Supplementary Material Available


Atal Anil Kumar
email: atal-anil.kumar@univ-lorraine.fr

Jean-Franc ¸ois
email: jean-francois.antoine@univ-lorraine.fr

Patrick Zattarin
email: patrick.zattarin@univ-lorraine.fr

Gabriel Abba
email: gabriel.abba@univ-lorraine.fr

Keywords: Static Equilibrium, non-linear optimization, cable-driven parallel robot

This paper presents the influence of payload and platform dimensions on the static equilibrium workspace of an under-constrained cable-driven robot with four cables taking into account the forces and the moments due to the forces acting on the moving platform. The problem is formulated as a non-linear optimization problem with maintaining static equilibrium as the objective function. The simulations are done in MATLAB. The maximum force on the cables, the payload acting on the platform and the dimensions of the moving platform are varied and their corresponding effects on the static equilibrium is studied. The obtained results are analyzed to finalize the design of the collaborative cable-driven robot to be installed in existing production lines for the agile handling of parts in a manufacturing industry.

I. INTRODUCTION

Cable-Driven Parallel Robots (CDPRs) belong to the class of parallel manipulators in which the rigid links of the robot is replaced by a set of cables driving a platform with the help of a winch consisting of a tensioning motor and spool or a linear actuator moving a pulley system [START_REF] Merlet | A portable, modular parallel wire crane for rescue operations[END_REF]. The position and orientation of the moving platform are controlled by the coordinated retraction and extension of the cables. CDPRs are very well suited for a number of applications like large-scale material handling and manufacturing process, building construction, rescue operations, aerial cameras, high-speed assembly and pick and place operations and so on [START_REF] Gosselin | Cable-driven parallel mechanisms: state of the art and perspectives[END_REF] because of a number of desirable features such as such as simple structure, large load capacity, high payload-to-weight ratio, and low inertia, compared to serial-link and other parallel type robots [START_REF] Verhoeven | Analysis of the workspace of tendon-based Stewart platforms[END_REF].

Several classifications of CDPRs are available in the literature

We would like to thank the Robotix Academy for funding this work as a part of the project funded by INTERREG V-A Grand Region program. [START_REF] Verhoeven | Analysis of the workspace of tendon-based Stewart platforms[END_REF]. According to the classification suggested by Cone [START_REF] Cone | Skycam-an aerial robotic camera system[END_REF] and Verhoeven [START_REF] Verhoeven | Analysis of the workspace of tendon-based Stewart platforms[END_REF] there are mainly three types of CDPRs namely over-constrained, fully-constrained and under-constrained. A cable robot with m-cables and n-degrees of freedom is said to be completely restrained or fully-constrained if it has one cable more than the number of degrees of freedom, i.e. m = n + 1. In such a type of arrangement, all degrees of freedom can be controlled through the cables. An over-constrained CDPR has the condition m ≥ n+1. An incompletely restrained or underconstrained cable robot is one in which the number of cables is less than or equal to the number of degrees of freedom, i.e. m ≤ n. Such systems have one feasible solution for cable tensions and mostly rely on gravity for keeping the cables taut.

In order to decrease complexity, cost, set-up time, likelihood of cable interference etc., it is advantageous to use CDPRs with a limited number of cables in several industrial applications in which the task to be performed requires a limited number of controlled freedoms or a limitation of dexterity is acceptable [START_REF] Abbasnejad | Direct Geometrico-static problem of underconstrained cable-driven parallel robots with n cables[END_REF]. Several studies are available on underconstrained CDPRs [START_REF] Fattah | On the design of cable-suspended planar parallel robots[END_REF]- [START_REF] Merlet | Wire-driven parallel robot: Open issues[END_REF]. Some important results on the stability analysis and kinematics of under-constrained parallel robots have been discussed by Merlet and Daney [START_REF] Merlet | Kinematic analysis of a spatial four-wire driven parallel crane without constraining mechanism[END_REF], Carricato and Merlet [START_REF] Carricato | Stability analysis of underconstrained cable-driven parallel robots[END_REF].

The major objective of the project is to design a collaborative cable-driven parallel robot to be installed in the existing lines of a manufacturing industry where available space is reduced by already installed machines. Hence, a cable-driven parallel robot with 4 cables is proposed as a solution. The number of cables is restricted to 4 for the ease of installation with simple anchor points in the environment and also to reduce the complexity of cable collision with the equipment installed on the shop floor. It is desired to have the actuators on the moving platform, unlike other designs in which the Influence o f p ayload a nd p latform d imensions on the static workspace of a 4-cable driven parallel robot motors are placed in the base or fixed structure. This paper presents the preliminary studies carried out to help the user in fixing the dimensions of the moving platform and the limits on the cable tension and the platform orientation by analyzing the available static equilibrium workspace for various configurations. The results from this study will help the designer in selecting the appropriate motor and its corresponding elements, the arrangement of cables according to the cable tensions, and the weight of the platform according to the payload. It will also give the an idea about the various regions where the robot performance is satisfactory and acceptable. Earlier works with such CDPRs take into account only the forces acting on the moving platform (point-mass). However, this work takes into account both the forces and the moments due to the forces acting on the moving platform with a given dimension. The results obtained from this study indicate that it is possible to achieve satisfactory results with 4 cables and the results obtained can be used in finalizing the parameters for the design of the CDPR.

II. KINEMATIC AND GEOMETRIC MODELING

The unilateral constraints imposed by the cables prevents the application of the modeling and analysis methods developed for conventional rigid link parallel manipulators to the cabledriven robots. The equations used in the modeling of CDPR are presented in this section [START_REF] Gosselin | Cable-driven parallel mechanisms: state of the art and perspectives[END_REF]. A general cable-driven parallel robot with m-cables is shown in fig. 1. The following assumptions are made to reduce the complexity of computation 1) The mass of the cables is negligible and the cables are non-elastic. 2) The moving platform is assumed to be a rigid body defined by its mass and inertia matrix. 3) The cable is assumed to be taut between points and is considered to be a straight segment.

As seen in fig. 2, a fixed reference frame also known as the base frame is attached to the base of a CDPR. A moving reference frame is attached to the moving platform and contains the point P to be positioned by the mechanism. The orientation of the moving frame with respect to the base frame is given by the rotation matrix Q. Also, from fig. 2, a i and b i are respectively defined as the vector connecting point O to point A i and the vector connecting point P of the The loop closure equations are first developed for each chain of the CDPR to identify the kinematic relationship between the cable space and body space for a CDPR. The length of the cable can be calculated by using the vector loop equation given by

- → d i = - → a i - - → b i -- → p (1)
where, -→ d i is the vector along cable i pointing from B i to A i and (i = 1, 2, ., .m). The length of the cable is calculated using

ρ i = - → d i 2 (2)
Also using eq. 2, the unit vector can be calculated by

di = ( - → a i - - → b i -- → p )/ρ i (3)
This basic model has been established in the early works on CDPRs and is purely based on the geometric analysis. This equation( eq. 3), however, does not take into account the ability of the mechanism to maintain the tensions in the cables. In order to address this issue, the concept of wrench matrix was introduced where statics of CDPR is considered.

When a tension t i is applied, the cable i exerts at point B i , a pure force t i di on the mobile platform. This pure force generates a moment -→ b i t i di at the reference point P of the mobile platform and the wrench (force/moment pair) applied at P by the i th cable is t i w i with the wrench w i defined as

w i = di - → b i t i di (4)
If w p is the total wrench applied at point P by the m-cables then the relationship between the tensions in the cables and the wrench is written in matrix form as

W - → t = w p (5)
where, -→ t is the vector of cable tensions and W is the 6 * m pose dependent wrench matrix.

Most of the analytical results proposed in the literature are based on this equation. By properly defining the vector w p , different solutions can be obtained to understand the workspace of CDPRs.

III. WORKSPACE ANALYSIS

The workspace analysis of cable driven robots is done taking into account the statics of the manipulator along with the change in cable lengths. There are many definitions available for workspace. The workspace of a CDPR is the set of position and orientation in which [START_REF] Tavolieri | A workspace analysis of a overconstrained cable-based parallel manipulator by using the interval analysis[END_REF] 1) The end-effector is controllable, 2) The tensions in the cables are positive, 3) The force values lie between a minimum and a maximum in order to avoid the cables break or become slack, 4) The end-effector is far from singularities, and 5) Cable wrapping is avoided.

Different types of workspaces have been defined and addressed in the literature. The current work however presents the static equilibrium workspace of the cable robot to be designed. The static equilibrium workspace (SEW) is defined as the set of postures (position and orientation of the endeffector) that the end-effector can attain statically (only taking gravity into account). It is important to have a good knowledge about SEW as it will help in the appropriate determination of the design parameters (motor speed, platform dimensions etc.).

In order to calculate the SEW of the CDPR, the static equilibrium condition of the cable robot is used to find the force on each cable. In other words, the equations given below which indicate the sum of forces ( F i , i = x, y, z) acting on the moving platform along the x, y, z direction and the sum of moments ( M i , i = x, y, z) acting on the moving platform along the x, y, z direction respectively, needs to be satisfied to achieve static equilibrium.

F x = 0, F y = 0, F z = mg (6) M x = 0, M y = 0, M z = 0 (7)
The only external force acting on the platform is the gravitational force (along with the mass of the platform and payload acting at the center of mass of the moving platform) and all other external forces and moments will not be discussed in this paper. The equations described above are formulated as a nonlinear optimization problem and is solved using fsolve in MATLAB with static equilibrium condition as the objective function.

IV. RESULTS AND DISCUSSION

The results of the simulation carried out in MATLAB are presented in this section. The SEW of the CDPR inside a room of dimension 5m * 5m * 3m is presented in the figures below. The constraints on the angle of inclination of the moving platform about x -axis and y -axis were fixed at ±30 • . The center of mass (CoM) of the moving platform was considered to be at a height of 0.2m below the cable attachment points while the CoM of the payload was considered at a height of 0.4m below the CoM of the platform. The step size was fixed along x and y axis to 0.1.

The search for possible workspace value was fixed reduced to the region of (0.3 < x < 4.7, 0.3 < y < 4.7, 0 < z < 2.5). The points were checked to see if they were belong to in the SEW by fixing various limits for the cable tension and payload of the platform as explained below. The dimensions of the platform were first varied and the corresponding percentage of SEW is presented in table I. The maximum tension on the cables for this study was set to be at 500N while the minimum tension was 1N . Since the final design includes the actuating motors on the moving platform, a dimension of 0.5m * 0.5m is considered to be a good compromise. Since a 4 cable-driven spatial parallel robot has 4 degrees of freedom (the position x, y, z and one angle) that can be controlled by the cables, it is important to have an idea about how the platform orientation behaves as the simulation is done for static equilibrium condition. The variation of cable tensions and the corresponding changes in the value of the angles about x-axis and y-axis when the moving platform is moving in a plane at the center of the room (z=1.5m) is shown in the figure below for a payload of 30kg (20 kg mass of the platform+10kg payload of the material) [START_REF] Kumar | Workspace Analysis of a 4 Cable-Driven Spatial Parallel Robot[END_REF]. It can be seen that the change in angles about x-axis and y-axis is symmetrical under the consideration that the moving platform has equal length and width. The plot of cable tensions also exhibits symmetrical behavior with respect to x and y axis.

Following this, the weight of the platform was varied from 20kg, 25kg and 30kg respectively. The simulation was carried out to observe the SEW for two maximum cable tensions namely 500N and 300N respectively. The simulation results are presented in the figures below. It can be seen from the figures that for a platform weight of 20kg, the available SEW is higher than the ones obtained for the other platform weight (25kg and 30kg). The main objective is to give further information to designer in selecting the components for the design. Following this, the CoM of the platform was varied from 0.2mto0.4m and the corresponding results obtained are shown in table II . It can be seen from the table that the SEW of the CDPR increases as the CoM of the platform is moved away from the cable attachment points. However, for the design to be done, the CoM is chosen to be 0.2m from the cable attachment points as the practical application of the other points will restrict the use of the CDPR for various heights. V. CONCLUSION The work presents the preliminary studies carried out in the design of a 4 cable-driven parallel robot. The static equilibrium workspace of the robot is calculated by varying the platform dimensions, payload and the corresponding tension limits for the cables. It can be seen that it is possible to have a sufficiently large workspace for such robots by limiting the orientation of the platform within a certain range and by making the appropriate choice of the dimension and payload. Future work will include the development of the dynamic model for the robot and implementation of the control law. The end final product of the current work is to have a prototype of the robot considered, in the agile handling of parts in a manufacturing environment.

Fig. 1 .

 1 Fig. 1. A simple arrangement of a CDPR

Fig. 2 .

 2 Fig. 2. Kinematic model of a CDPR

Fig. 3 .Fig. 4 .

 34 Fig. 3. Tension distribution (N) in 4 cables for z=1.5m for a payload of 30kg (20kg mass of platform + 10kg mass of the payload material)

Fig. 5 .

 5 Fig. 5. Variation of angle (in degrees) about y-axis for z=1.5m

Fig. 6 .

 6 Fig. 6. SEW of the CDPR for a payload of 30kg (20kg platform+ 10kg payload) acting on the platform for a maximum force of 500N (a) 3-D representation and (b) XZ view of the workspace

TABLE I %

 I OF SEW FOR PLATFORM WITH DIFFERENT DIMENSIONS

	Platform dimensions(m)	0.5 * 0.5	0.75 * 0.75 0.6 * 0.4	0.8 * 0.6
	SEW%	58	48	58	46

TABLE II %

 II OF SEW FOR DIFFERENT CENTRE OF MASS FOR A PLATFORM OF DIMENSION 0.5m * 0.5m

	Centre of Mass (CoM) in m	0.2	0.25 0.3	0.4
	SEW%	58	65	72	73

nd Robotix-Academy Conference for Industrial Robotics (RACIR) 2018, University of Luxembourg, Luxembourg, 04-05 June 2018