Nicolas Philippe

Mikael Salson

Thierry Lecroq

Martine Léonard

Thérèse Commes

Eric Rivals
email: rivals@lirmm.fr

Mikaël Salson

Querying large read collections in main memory: a versatile data structure

published or not. The documents may come

Querying large read collections in main memory: a versatile data structure Philippe et al.

Background

Next-generation sequencing technologies are presently being used to answer key biological questions at the scale of the entire genome and with unprecedented depth. Whether determining genetic or genomic variations, cataloging transcripts and assessing their expression levels, identifying DNA-protein interactions or chromatin modifications, surveying the species diversity in an environmental sample, all these tasks are now tackled with High Throughput Sequencing (HTS) and require different, but computer intensive bioinformatic analyses. Typically, a recent RNA sequencing experiment (RNA-Seq) produces about 8 million reads of 75 base pairs each [START_REF] Maher | Transcriptome sequencing to detect gene fusions in cancer[END_REF], but both the yield and read length will increase [START_REF] Blow | Transcriptomics: The digital generation[END_REF].

Mapping the reads against a reference genome provides the genomic positions of mapped reads. For instance with RNA-Seq reads, these positions allow to know whether a gene is expressed in the studied condition. The set of mapped positions represents only part of the information needed to analyze the reads, and it can be obtained only if a genome is available. Indeed, other important information are contained in the read collection itself. For instance, to determine the frequency of haplotypes at a SNP position, one needs to align the reads related to this position. These can be obtained by considering for some length k,t h ek-mers overlapping the SNP and searching for the reads sharing this k-mer. This procedure is applicable even in the absence of a reference genome, and similar ones can be designed to search for a binding motif in ChIP-Seq reads, to determine with RNA-Seq data whether different regions of a messenger RNA sequence are susceptible to be differentially expressed, etc.

For tasks like assembly or read clustering, one needs to determine reads overlapping each other or that align partly one to another. Numerous works on similarity search algorithms have developed seed-and-extend strategies and shown that it can be performed efficiently by searching common k-mers between two sequences [START_REF] Altschul | Basic local alignment search tool[END_REF][START_REF] Burkhardt | q-gram Based Database Searching Using a Suffix Array (QUASAR[END_REF].

Surely, now and even more in the near future, we will need efficient indexing data structures to store and query large collections of reads in main memory. Up to now, a lot of computational research has been devoted to read mapping, and the most efficient tools owe their efficiency to the use of involved genome indexing data structures, like the Burrows-Wheeler transform [START_REF] Ferragina | Opportunistic data structures with applications[END_REF]. On the other hand, the question of read indexing remains quite unexplored, although the improvements in sequencing throughput suggest that such structures will b e c o m eac o m p u l s o r yp a r to ff u t u r er e a da n a l y s i sp r ograms. A sign supporting this view: even mapping programs now start to index both the genome and the kmers of the reads to boost efficiency and accuracy [START_REF] Hach | mrsFAST: a cache-oblivious algorithm for short-read mapping[END_REF].

Numerous works have presented data structures to index a single text, like the well known Suffix Tree (ST) or the Suffix Array (SA) [START_REF] Weiner | Linear Pattern Matching Algorithms[END_REF][START_REF] Manber | Suffix Arrays: A New Method for On-Line String Searches[END_REF]. These enable the so-called locate query, that is to locate all occurrences of a pattern P either from its sequence or from a position j of occurrence in the text, as well as count query to obtain the number of occurrences of P. These structures can be adapted to index a set of texts, where each text differ from each other; the structures are then called generalized Suffix Tree [START_REF] Gusfield | Algorithms on Strings[END_REF], or generalized Suffix A r r a y( g S A )[ 1 0 ] .T h i si sd o n eb yc o n c a t e n a t i n ga l l texts and adding a separator symbol that does not belong to the alphabet (e.g., a $ for the DNA alphabet) after each text [START_REF] Gusfield | Algorithms on Strings[END_REF], or directly [START_REF] Shi | Suffix Arrays for Multiple Strings: A Method for On-Line Multiple String Searches[END_REF]. Then it requires to store the length of each text in an additional array to c o r r e c t l ya n s w e rl o c a t eq u e r i e s .S u c ha l g o r i t h m sh a v e not been adapted to collections of texts, where two texts may be equal in sequence but differ in their identifier. The reads obtained from sequencers form a collection, not a set.

When the total text is too large, compressed indexes reduce the memory needed by storing not all, but only a certain proportion of the text positions. Compression is obtained by sampling the positions to be stored, while non sampled positions need to be recomputed at run time. This enables the user to control the balance between amount of memory and query time. Hence, compression has an impact on the time needed to compute a query. Ferragina et al. report in a large practical evaluation of compressed text indexes, that the query time of all tested compressed indexes are between 100 and 1,000 times slower than with a plain SA for an index that is 5 times smaller [START_REF] Ferragina | Compressed text indexes: From theory to practice[END_REF]. The FM-index [START_REF] Ferragina | Opportunistic data structures with applications[END_REF] is used to index all chromosomes in mapping applications [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF]. However, the scalability of neither plain nor compressed indexes to collections of millions of texts has not been investigated so far. We thus address the question of indexing large collections of reads with an uncompressed index and compare its performance to a generalized suffix array and a hash table. Our structure aims to save space compared to those indexes while globally retaining queries as fast. Thus we avoid the pitfall of compressed indexes which are less space consuming but slower by orders of magnitude.

I nt h i sw o r k ,w ep r o p o s ean e wd a t as t r u c t u r et o index reads, an algorithm to build the structure, and procedures to query it. Our structure, named Gk arrays, is kept in main memory once built and repeatedly accessed to answer different kinds of queries like "given a k-mer, get the reads containing this k-mer (once/at least once)".O n ec a na s kb o t hf o rt h ek-mer positions or simply for the reads containing it, which can prove useful in different applications. We focus on cases where millions of queries need to be computed; clearly, memory usage will be the key issue. An alternative solution is to adapt some uncompressed indexing structures designed for long texts (suffix tree or suffix array [START_REF] Gusfield | Algorithms on Strings[END_REF][START_REF] Homann | mkESA: enhanced suffix array construction tool[END_REF]). We compare Gk arrays to such an alternative and show experimentally that they process queries fast, while requiring much less memory ( b e t w e e n2 / 3a n d1 / 3o fas u f f i xa r r a ys o l u t i o n ) .W e also perform experimental comparisons against a method using hash table: it shows that while the hash table method can answer quickly to queries it does not scale to large collections of reads.

If in biology the term k-mer is preferred, computer scientists rather use the equivalent words of k-factor or k-substring; we will stick to the term k-mer. The Gk arrays allow to answer queries related to an input kmer; let us call these k-mer queries. Before entering the algorithms description, we list below the applications of k-mer queries in the analysis of High Throughput Sequencing data. The Results section will first present our data structure, its construction algorithm and the procedures to answer k-mer queries, then detail the experimental comparisons.

Finally, we discuss the advantages of our structure and conclude with future developments.

Note that this study does not tackle the question of read mapping, it focuses on read indexing.

Queries and Applications

Let us give an informal presentation of the problem. We a r eg i v e nac o l l e c t i o no fq reads of length m and a length of substring k such that k ≤ m.

Suppose one is given a string f of length k;o n ed o e s not know whether it appears in some of the reads or not (i.e., whether f is a substring of some read). In the Algorithm section, we describe a data structure in which all substrings of length k of the reads are ordered lexicographically. Hence, one can search for f using a dichotomic search in O(k log((mk +1 ) q))) worst case time in this structure (the dichotomic search is the standard procedure in this context [START_REF] Manber | Suffix Arrays: A New Method for On-Line String Searches[END_REF][START_REF] Gusfield | Algorithms on Strings[END_REF]), and determine whether at least one read contains f as a substring and at which position. If not, the answers to the queries below, which are all related to a sub-collection of reads containing f, are trivially the empty set or zero. Otherwise, one knows that f occurs in some read r of the collection at position j, and wishes to get some information on the other reads where f occurs. One wants to answer the following questions:

Q1: In which reads does f occur? Q2: In how many reads does f occur? Q3: What are the occurrence positions of f in the reads?

Q4: What is the number of occurrences of f in the reads?

Q5: In which reads does f occur only once? Q6: In how many reads does f occur only once? Q7: What are the occurrence positions of f in each read where f occurs only once?

Q8: What is the number of occurrences of f in the reads where it occurs only once?

We state several remarks about the queries before dwelling on applications.

1. The queries go by pairs: the first one computes a set of positions or read indices, while the second computes the cardinality of that set. 2. Note the clear semantic difference between Q1/ Q2 and Q3/Q4. The answer to Q1 yields the identifiers of the reads in which f occurs, while that to Q3 gives also all its positions in the read. This clearly differs since f may occur several times in a read (e.g., if f is a poly-A sequence). Sometimes the positions are needed, sometimes only the reads (see below). 3. Queries Q5-Q8 are versions of Q1-Q4 constrained to a single occurrence of f in the reads. Of course other variants can also be computed, e.g.w h e r et h e number of occurrences is limited by a user defined threshold. Since f is constrained to occur only once in each read, Q6 and Q8 are equivalent, and we will mention only Q6 in the sequel. [START_REF] Burkhardt | q-gram Based Database Searching Using a Suffix Array (QUASAR[END_REF]. The data structure we propose is intended to be kept in memory and used for multiple queries.

Although this paper focuses on the data structure, its efficiency, and on the algorithms to solve these type of queries, it is important to list applications of these queries. In which context of read analysis, can one use such queries? Note that in such context, k is smaller than the read length. Theoretical and empirical investigations show that for instance, with k ≥ 19 or 20, kmers indicate in average a single genomic location in the human genome [START_REF] Philippe | Using reads to annotate the genome: influence of length, background distribution, and sequence errors on prediction capacity[END_REF]. Such values of k can be computed depending on the genome length. Translated to reads or sequences: it is unlikely that two reads sharing a k-mer were not sequenced from the same part of the DNA. In other words, sharing a k-mer is a witness for having a common genomic origin.

Mutation detection

Putative mutations (SNP, somatic mutations, small indels) are indicated by differences between a read and a reference genome. Once the reads have been mapped to the reference genome, one analyzes the sub-collection of reads that covers a genomic position to count how many reads support the variation observed in the read or that observed in the genome. If one considers the two substrings of length k centered on this mutation position, one in the read and one in the genome, answering Q2 for these substrings will give an approximate count of these two haplotypes. If one needs the corresponding reads, then Q1 is the appropriate query. If only a single, or a few reads, share this k-mer, then a sequence error might be suspected [START_REF] Salmela | Correction of sequencing errors in a mixed set of reads[END_REF].

Local coverage

Suppose one is given a target sequence, which can be a read or an external sequence. For each of its k-mer, let us call the local coverage, the number of reads sharing this k-mer (this requires a dichotomic search). The local coverage profile (i.e. a histogram of the local coverage) along the target sequence provides useful information in various contexts. For a known mRNA and an RNA-Seq experiment, the average local coverage on all k-mers is a proxy for the expression level of the target, while the p r o f i l ee n a b l e so n et od i s t i n g u i s ht h et a r g e t 'ss u bregions expressed at different levels [START_REF] Denoeud | Annotating genomes with massive-scale RNA sequencing[END_REF][START_REF] Trapnell | Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF]. In another context, with a genomic library, taking reads as queries and looking at their local coverage profile may help to detect those overlapping the extremity of a repeated or transposable element. This may prove useful to study the distribution and evolution of these elements in the genome.

Clustering and assembly without a reference genome

As for Expressed Sequence Tags, it is suitable to cluster and assemble RNA-Seq reads to compute the various transcripts expressed in the assayed library [START_REF] Denoeud | Annotating genomes with massive-scale RNA sequencing[END_REF][START_REF] Trapnell | Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF]. It is necessary to detect near exact alignment between pair of reads, and this is usually performed efficiently by filtration using seeds. In such case, very efficient and sensible seeds are exact shared k-mers [START_REF] Burkhardt | q-gram Based Database Searching Using a Suffix Array (QUASAR[END_REF]. Here, the subcollection of reads sharing a k-mer with a given read, as well as the k-mer positions, can be obtained using query Q3. The answer to Q4 can help guiding the clustering process.

Similar needs of query occur in the assembly of genomic reads [START_REF] Miller | Assembly algorithms for next-generation sequencing data[END_REF][START_REF] Conway | Succinct Data Structures for Assembling Large Genomes[END_REF]. To know with which reads one can assemble a given read without ambiguity, one may perform query Q7 using k-mers at the 5' or 3' extremities of the read. The obtained occurrences together with their positions will indicate the matching reads and the relative positions of read pairs for assembly.

Our application list provides examples and is by no means exhaustive. We could also mention for instance the estimation of the target genome length in assembly context, which uses k-mer counting [START_REF] Marcais | A fast, lock-free approach for efficient parallel counting of occurrences of k-mers[END_REF]. Clearly, these applications are beyond the scope of this paper. However, these paragraphs underline that the proposed data structure suits the needs of read processing in various application contexts, and will provide a unified framework for building read analysis programs.

Results and Discussion

This section contains the main contribution: a data structure to index large read collections, the Gk arrays.

To describe it, we first introduce the notation, formalize the queries, exhibit the index data structure, give its construction algorithm, and the procedures for answering all queries. This makes the content of the Algorithms section. Then, in the Comparison section we investigate its practical usability compared to two alternatives: one based on a generalized Suffix Array (SA) and another based on a hash table. This includes theoretical and practical comparisons.

Algorithms

Here, we detail the algorithms to build the Gk arrays and to answer the queries. We start by defining more formally the queries we want to answer and introduce the necessary notation.

Notation and definition of the queries

Let Σ be an alphabet of size σ. Σ* denotes the set of words, strings or sequences over Σ and, for any integer n, Σ n denotes the set of words of length n over Σ.F o ra word x,| x|d e n o t e st h elength of x. Given two words x and y,w ed e n o t eb yxy the concatenation of x and y. For every 0 ≤ i ≤ j ≤ |x|-1 ,x[i] denotes the (i +1 ) th element of x,a n dx[i.. j]d e n o t e st h esubstring

x[i]x[i + 1 ]...x[j]
. Let ≤ L denotes the comparison operator for the lexicographic order on words. Lexicographic ranks start from zero and all arrays are indexed from zero. For any finite set A, we denote its cardinality by #A.

The input consists a list R =( r 0 ,... ,r q-1 )o fq short sequences of length m,c a l l e dreads,w h i c ha r en o t necessarily distinct. We know that m, k, q N satisfy m ≥ k >0 .

A k-long substring of a word is called a k-mer. For any u Σ*, we denote by

F k (u) the set of k-mers in u: F k (u) ={ v Σ k | ∃p [0, |u|-k] such that v = u[p..p + k - 1]}. Let f Σ k and let us denote the set of indexes of the reads in which f occurs by Ind k (f)={ j [0, q[| f F k (r j )},
(f) (resp. UPos k (f)). Formally, UPos k (f)={ ( j, ℓ)| r j [ℓ..ℓ + k -1 ]=f and ∀i ≠ ℓ, r j [i..i + k -1 ]≠ f}, and UInd k (f)={ j |( j, ℓ) UPos k (f)}. Let i [0, q[, j″ [0, m -k + 1[,
and let f be the k-mer starting at j″ in read r i .N o t et h a th e r ew er e q u i r et h ek n o w l e d g eo ft h ep a i r (i, j″), which defines the k-mer f.N o w ,t h es e v e nk-mer queries can be formally defined as computing

Q1: the set Ind k (f ) Q2: #Ind k (f ), the cardinality of Ind k (f ) Q3: the set Pos k (f ) Q4: #Pos k (f ), the caradinality of Pos k (f ) Q5: the set UInd k (f ) Q6: #UInd k (f ), the cardinality of UInd k (f ) Q7: the set UPos k (f )
Clearly, it appears (see Additional File 1: Proof and queries' algorithms) that the algorithms to compute UInd k (f), resp. UPos k (f), for answering Q5/Q7, simply filter Ind k (f), resp. Pos k (f), on the fly, and are thus similar to the algorithms for Q1/Q3. For place sake, we will only detail the solutions for Q1-Q4 in the sequel.

The index structure

Our algorithm relies on four arrays that allow to query the k-mers of all reads. Hence, we define a word made of the concatenation of all reads: C R = r 0 r 1 ... r q-1 .O f course, a k-mer that overlaps two reads in C R is not necessarily a k-mer of some read. Hence, we introduce a system to renumber the positions of interest in C R .T h e rationale behind is to save place in the Gk arrays by discarding the positions of overlapping k-mers in C R .L e t us denote by r the number of distinct k-mers of all reads, and for the sake of legibility we set m := m -k + 1 and q := q m = q(m -k +1) (the number of interesting positions in a read and in C R , respectively). We call:

• P-position, a starting position in C R of a k-mer that is not overlapping two reads, i.e.a ne l e m e n to f P pos := 0, qm \ j | j mod m ≥ m . • g, the function that renumbers P-positions in order such that their index are consecutive; g is defined by:

g : P pos → Q pos j → j m × m +(j mod m) .
• Q-position,a ni m a g eo faP-position by g(.), i.e.a n element of Q pos := g(P pos )=[0,q[.N o t et h a tt h es e t Q pos is not a query.

Clearly, P pos and Q pos have the same cardinality q, and as (j ≠ j') implies g(j) ≠ g(j'), g is bijective. Hence, g -1 exists and maps a Q-position back to its corresponding P-position in C R . Proposition 1 explicits the conversion between a positioned occurrence and a P-position.

Proposition 1. Let (j, ℓ) with j [0,q[, ℓ ∈ [0, m[ be a positioned occurrence of a k-mer in a read. The corresponding P-position in C R is jm+ℓ. Conversely, let j' be a P-position, the corresponding positioned occurrence in a read is (⌊j'/m⌋, j' mod m).

This numbering system is important for it allows us to go back and forth between a positioned occurrence in a read, its corresponding P-position in C R , and its Q-position that will be stored in our arrays.

Let j be a Q-position. We denote by s Q (j), resp. f Q (j), the suffix, resp. the k-mer, of C R beginning at the P-

position g -1 (j), i.e. s Q (j)=C R [g -1 (j)..qm -1 ]a n df Q (j)= C R [g -1 (j)..g -1 (j)+k -1].
We call s Q (j)aP-suffix. Note that all suffixes beginning at P-positions have different length and are pairwise distinct; thus, there are q such suffixes and they all have a different lexicographic rank. However, this may, and in real data applications will, not be the case for the k-mers, i.e.t h ef Q (j). We call the set {f Q (j)|j Q pos }t h es e to fP k -factors,w h o s ec a r d i nality is r with our notation. Now, we define the Gk arrays: GkSA (Generalized k Suffix Array) is a modified Suffix Array of C R that lexicographically sorts only the Psuffixes, GkIFA (Generalized k Inverse Factor Array) is a modified Inverse Suffix Array (ISA) that stores for each Qposition, in position order in C R , the lexicographic rank of the P k -factors starting at the corresponding Pposition,

GkCFA (Generalized k Counting Factor Array) is an array that associates to a k-mer (actually, to its rank) its number of occurrences at P-positions in C R , GkCFPS (Generalized k Counting Factor Prefix Sum) stores the prefix sums of GkCFA.S i n c eGkCFA and GkCFPS are equivalent only one of them is necessary at a time.

Formally, the definitions are (see Figure 1 for an example and Figure 2): ). The index is composed of three tables and uses a fourth one during construction (GkSA, GkIFA, GkCFA, and GkCFPS). The first table shows the starting indices of k-mers in the text made by the concatenation of all reads, C R , the SA built on C R , and the function g that renumbers P-positions of C R to make them consecutive. P-positions are {0,1,2,3,4,7,8,9,10,11,14,15,16,17,18}; all other positions, those starting positions where the k-mer overlaps two reads, are displayed with a gray background (lines j and SA[j]). Line SA refers to the usual Suffix Array of C R . The k-mer caa occurs 4 times in C R at positions 2, 7, 12 and 16. Among those, only 2, 7, and 16 are P-positions. The lexicographic rank of the P k -factor starting at position 16 is given by GkIFA[g [START_REF] Denoeud | Annotating genomes with massive-scale RNA sequencing[END_REF]] = GkIFA12 = 7, and the number of occurrences of the P k -factor caa is given by GkCFA7, which equals 3. The positions of these occurrences are thus obtained by the set {g -1 (GkSA[j]) | GkCFPS[7 -1] ≤ j < GkCFPS7} = {2,7,16}. See also Figure 2.

• For i a suffix lexicographic rank and j a Q-position (i.e. i, j Q pos ), GkSA[i]=j iff s Q (j)
GkCFA[i]=# { j Q pos |f Q (j)=f Q (GkSA[i])}, • For i a k-mer lexicographic rank (i.e. i ∈ [0, r[), the definition of the prefix sum is GkCFPS[i]= i t=0 GkCFA[t] and GkCFPS[-1] = 0.
Remark 1. The array GkCFPS is not essential to the algorithm: it is solely there to avoid multiple, time consuming computations of prefix sums over GkCFA (see GkCFPS definition above). Moreover, any value of GkCFA can also be accessed in constant time using

GkCFA[i]=GkCFPS[i]-GkCFPS[i -1]
. Thus, GkCFPS will be kept in memory to replace GkCFA.

We give some useful properties of Gk arrays.

Proposition 2. For i ∈ [0, r[, GkCFPS[i]=# { j Q pos | f Q (j) ≤ L f Q (GkSA[i])} (Proof by induction).
In other words, GkCFPS[i] is the number of P k -factors having lexicographic rank less than or equal to i.S i n c e GkSA is sorted on the lexicographic order of the P-suffixes, it is also sorted on the lexicographic order of the P k -factors. Hence, we get:

Proposition 3. Let f Σ k such that Ind k (f) ≠ ∅.
All occurrences of f have the same rank among the P k -factors, and are stored consecutively in GkSA.

Construction algorithm

First, we detail the algorithm for building GkSA,a n d then the one computing GkIFA and GkCFA. Computation of GkSA We first build the full Suffix Array (SA) of C R using a linear time and space algorithm. Since |C R | = mq this first step can be done in O (mq). Then GkSA is obtained from SA by selecting only the P-positions and by renumbering them to Q-positions using function g.T h i ss e c o n ds t e pi sp e r f o r m e di nO (mq) time and space. Moreover, GkSA is built in place of the Suffix Array: our algorithm allocates only the memory for the SA table. When answering Q1/Q2, each read where a given P k -factor occurs should be counted only once (even if the P k -factor occurs more than once in the read). Similarly, for Q5/Q6, we count only reads where a given P k -factor occurs exactly once. To avoid using masks on the reads, we sort in increasing order the values of GkSA corresponding to P k -factors sharing the same lexicographic rank (see Table in Additional File 1: Proof and queries' algorithms). The values that have to be sorted are Q-positions, i.e. integers, thus the sort can be performed in linear time on values of GkSA using e.g. radix sort [START_REF] Cormen | Introduction to Algorithms[END_REF]. The whole process takes O (mq) time and space. Computation of GkIFA and GkCFA Algorithm 1 shows how to compute jointly GkIFA and GkCFA.I t s correctness proof is given in Additional File 1: Proof and queries' algorithms.

Algorithm 1: Computation of GkIFA and GkCFA. The comparison between two P k -factors (line 8) is naively performed in O(k) time, and is the only instruction of the inner loop that takes more than constant time. Hence, the computation of both GkIFA and GkCFA is performed in O((m-k )qk)t i m e .L e tu s emphasize the simplicity of the algorithm, which explains the fast construction times obtained in practice.

Data: GkSA, C R , k, q Result: GkIFA and GkCFA 1 begin 2 GkIFA[GkSA[0]] 0; 3 GkCFA[0] 1; 4 t 0; 5 foreach i ∈ [1, q[ do 6 j GkSA[i]; 7 j' GkSA[i -1]; 8i f f Q (j) ≠ f Q (j') then 9 t t +1 ; 10 GkCFA[t] 0; 11 GkIFA[j] t ; C R r 0 r 1 r q-1 (a) (b) (c) j g GkIFA g( j) GkCFPS GkIFA[g( j)] GkSA GkCFA[GkIFA[g( j)]]
Remark 2. Once the values of GkCFA have been calculated, one can compute the values of GkCFPS in-place in O((mk)q) time (see Remark 1).

Answering the queries

Assume the Gk arrays have been built in a preprocessing step (see section Construction algorithm); we show how to answer the first four queries, starting with Q4 and Q3. Let i ∈ [0, q[, j ′′ ∈ [0, m[,a n dl e tf be the k-mer starting at j" in read r i . This occurrence of f in C R is found at P-position j':= im + j" and the corresponding Q-position is j := g(j'). Q4: Computing the cardinality of Pos k (f) First, we need to find the lexicographic rank of f among the P kfactors, which we obtain directly by setting t := GkIFA

u f = GkCFPS[t]a n dℓ f = GkCFPS[t -1 ] .
Since GkSA is indexed from 0, the starting Q-positions of occurrences of f are comprised in the range [ℓ f ,u f ] in GkSA. The corresponding P-positions are obtained using g -1 (.) and are then transformed into positioned occurrences with Proposition 1. This proves Theorem 2.

Theorem 2. Let f be a k-mer of a read occurring at Qposition j in C R . Then, its lexicographic rank among the P k -factors is t :

= GkIFA[j]. If we set u f := GkCFPS[t] and ℓ f := GkCFPS[t-1] then 1. the starting P-positions of f's occurrences in C R are {g -1 (GkSA[ℓ]) | ℓ [ℓ f ,u f [}, 2. Pos k (f)={ ( g -1 (GkSA[ℓ])/m , g -1 (GkSA[ℓ]) mod m) | ℓ [ℓ f ,u f [}, 3.# Pos k (f)=u f -ℓ f .
Given Theorem 2, the queries regarding Ind k (f) can be answered as follows:

Q1:

Ind k (f): = { g -1 (GkSA[ℓ])/m | ℓ [ℓ f ,u f [}, Q2:
by counting the elements of Ind k (f) while computing it.

The algorithms for Q1, Q3, and Q4 are given extensively in Algorithms 2, 3, and 4. The algorithms for all other queries are included in Additional File 1: Proof and queries' algorithms.

To answer Q7, one computes Pos k (f)a n ds c a n si to n the fly to remove reads (or the positioned occurrences) having strictly more than one occurrence of f. A similar approach solves Q8, and Q5. Variants of these queries where the number of allowed occurrences is constrained by a parameter can be answered similarly. Complexity Answering Q1-Q3 or Q5-Q8 requires to scan the values in GkSA inside the range corresponding to the k-mer f, which can be performed in O(occ_Reads (f)) time, where occ_Reads(f) denotes the occurrence number of f in the reads. Query Q4 is computed in constant time using GkCFPS.

Algorithm 2:

Q1 (Ind k (f)) Data: f ∑ k , j P pos such that C R [j. .j+ k -1 ]=f Result: The set Ind k (f) 1 begin 2 Ind k empty set; 3 t GkIFA[j]; 4 ℓ f GkCFPS[t-1]; 5 u f GkCFPS[t]; 6 prev -1; 7 foreach i [ℓ f ,u f [ do 8 readIndex g -1 (GkSA[i])/m ; 9i f readIndex ≠ prev then 10
Add readIndex to Ind k ; prev readIndex; 11 return (Ind k );

Algorithm 3: Q3 (Pos k (f)) Data: f ∑ k , j P pos such that C R [j..j+ k-1] = f Result: The set Pos k (f) 1 begin 2 Pos k empty set; 3 t GkIFA[j]; 4 ℓ f GkCFPS[t-1]; 5 u f GkCFPS[t]; 6 foreach i [ℓ f ,u f [ do 7 readIndex g -1 (GkSA[i])/m ; 8 posInRead g -1 (GkSA[i]) mod m; 9
Add the pair (readIndex, posInRead) to Pos k ; 10 return (Pos k ); Algorithm 4: Q4 (The cardinality of Pos k (f))

Data: f ∑ k , j P pos such that C R [j..j+ k -1 ]=f Result: The cardinality of Pos k (f) 1 begin //GkCFA[t]=GkCFPS[t]-GkCFPS[t-1] 2 t GkIFA[j]; 3 return (GkCFA[t]);

Practical considerations: implementation and variable read length

The value of k, which determines the length of k-mers used for querying the collection of reads, is a parameter of our index. However, the Gk arrays remain flexible. If for the simplicity of the presentation we have assumed until now that all reads have the same length, the whole structure can be adapted to a collection of reads having variable length. Indeed, since some sequencing technologies produce variable-length reads (e.g. Roche 454 ® ), this adaptation is an important issue of versatility. Indexing variable-length reads We show how our method can be slightly adapted to tackle this problem. Remind that the Gk arrays consider the string C R ,t h e concatenation of all reads, and save place by discarding positions at which a k-mer overlaps two reads. This was done efficiently by converting any read position, or Pposition, into a Q-position, and conversely, using function g. Up to now, this function relies on the fact that t h er e a dl e n g t hi sf i x e d .T h u s ,w en e e dt om o d i f yi t s definition to accommodate different read lengths. For this, we use a bit vector F,a sl o n ga sC R , to record which positions in C R are P-positions: j is a P-position iff F[j]=1 .W ei m p l e m e n ti ta sav e c t o rh a v i n gr a n k and select capabilities [START_REF] Munro | Tables[END_REF][START_REF] Raman | Succinct indexable dictionaries with applications to encoding k-ary trees and multisets[END_REF]. We define these operations as

• rank 1 (F, i) is the number of ones in F[0..i].
• select 1 (F, i) is the position of the i-th one in F (or | F| if there is less than i ones in F).

These operations can be performed in constant time, and F can be stored in a compressed form needing only |F|H 0 (F)+o(|F|)b i t s ,w h e r eH 0 is the zero-th order empirical entropy of F. Then computing g(j)a n dg -1 (j) can be easily performed with a single rank or select query. Indeed, we have g(j)=r a n k 1 (F, j), and g -1 (j)= select 1 (F, j). Finally, using little extra memory, Gk arrays can also handle variable-length reads. Implementation Gk arrays are available as a reusable C ++ library under a Cecill C licence (GPL compliant). It accepts standard formats for the input read collection (FASTA, FASTQ). Depending on the number of k-mer positions, the user should turn on the 64 bit encoding at compilation. It allows to process data sets of more than 2 31 positions. Default is set to 32 bit encoding. Another compilation option can be activated to handle variable-length reads (typically Roche 454 ® datasets), otherwise by default Gk arrays process fixed length reads.

The data structure construction and queries algorithms are coded in standard C and C++. To reduce memory consumption, the full SA of C R is built using libdivsufsort library https://code.google.com/p/libdivsufsort/, which was chosen for its efficiency and low memory usage (see https://code.google.com/p/libdivsufsort/ wiki/SACA_Benchmarks for a benchmark of up-to-date SA construction algorithms). However, its worst case time complexity is not linear in the length of the input sequence. Also the sort of values in GkSA inside each range corresponding to one P k -factor is performed with the quicksort algorithm. A linear time construction of the array GkIFA is possible by using an LCP array (array storing the length of the Longest Common Prefixes between two consecutive suffixes in the lexicographic order). However, building this array would need at least 9mq bytes with Manzini's algorithm [START_REF] Manzini | Two Space Saving Tricks for Linear Time LCP Array Computation[END_REF].

We implemented two versions of the Gk arrays: one which indexes only fixed-length reads, and another for variable-length reads. When not stated otherwise, Gk arrays refers to the implementation for fixed-length reads. For managing variable-length reads we used Sux http://sux.dsi.unimi.it/, an implementation of bit vectors with rank and select operations.

Theoretical and experimental comparisons

The sequencing capacity of new technologies continues to improve. Managing ever increasing read collections will be a major bottleneck in the bioinformatic analysis of High Throughput Sequencing data. The Gk arrays implement one solution to read indexing. If plain, as well as compressed, indexing data structures have been described in the litterature (cf. Introduction), their ability to handle large read collections have not been investigated. As we seek to optimise in practice the memory consumption, the construction time, and query running time, we will compare Gk arrays to two other uncompressed indexes: a generalized SA (gSA) and hash tables. We choose these two alternatives for they represent different approaches to read indexing. Among the uncompressed text indexes that have been generalized to handle a set of texts, the gSA is reckoned to be one of the most memory efficient and has been preferred to hash tables or the suffix tree in other contexts [START_REF] Gusfield | Algorithms on Strings[END_REF][START_REF] Kurtz | Versatile and open software for comparing large genomes[END_REF]. On the other side, the optimisation of web search engines have triggered recent development of highly efficient hash tables, like Google sparse hash http://code. google.com/p/google-sparsehash or the hash tables from SGI extension of the C++ Standard Library http://www. sgi.com/tech/stl. It is thus instructive to also compare Gk arrays to state of the art hash tables. As explained in Introduction, compressed indexes save memory but induce much longer running times to answer queries compared to plain indexes, and have been excluded from this comparison. Nevertheless, designing efficient compressed read indexes is a challenging future research avenue, which could be addressed by compressing the Gk arrays.

A generalized Suffix Array (gSA) solution

We detail here the solution based on a generalized Suffix Array (gSA) to index a collection of reads, all reads having the same length. We call it the gSA solution.I n fact it indexes the string made of the concatenation of all reads, C R . The preprocessing consists in building the generalized Suffix Array (gSA), the Inverse Suffix Array (ISA), and the Longest Common Prefixes (LCP) array of C R . The gSA is built using the same algorithm than for Gk arrays (libdivsufsort). The ISA is built by scanning the gSA in mq time, while the LCP array is also constructed in linear time using an efficient algorithm [START_REF] Kasai | Linear-Time Longest-Common-Prefix Computation in Suffix Arrays and Its Applications[END_REF]. The tables are built in this order and add up in term of memory footprint.

In Figure 3(a) and 3(b), we compare the time and space complexities of gSA and Gk arrays solutions. Since both start by building gSA(C R ) and this is the dominant term of the time complexity, we obtain O(mq) time complexity: the space occupied during the construction of that table alone is 4.02mq, while it amounts to 4mq once built [START_REF] Puglisi | A taxonomy of suffix array construction algorithms[END_REF]. The last three columns of these tables show how the cumulated memory footprint evolves after each step during construction. We also monitored the memory footprint evolution during the construction of gSA and of Gk arrays and illustrate these graphically in Figures 4(a

Gk arrays generalized Suffix Array Q1 O(occ Reads(f)) O(q + occ C R (f)) Q2 O(occ Reads(f)) O(q + occ C R (f)) Q3 O(occ Reads(f)) O(occ C R (f)) Q4 O(1) O(occ C R (f)) Q5 O(occ Reads(f)) O(q + occ C R (f)) Q6 O(occ Reads(f)) O(q + occ C R (f)) Q7 O(occ Reads(f)) O(q + occ C R (f)) (c)
Figure 3 Comparing the complexities of the Gk arrays and generalized Suffix Array based solutions. Comparing the complexities of Gk arrays and of the generalized Suffix Array solutions. A complexity is an expression that evaluates the running time or memory usage in function of parameters describing the input size. The construction time and space complexities of the index for q reads of length m having r distinct kmers are given for the generalized SA in (a), and for the Gk arrays in (b). We detail the cumulative space complexity during the construction of the gSA, and after the main steps of the construction algorithms. I.e.: once the gSA, the ISA, and the LCP arrays are built in (a), and once GkSA, GkIFA, and GkCFPS are built in (b). In (c) we give the time complexities for answering queries Q1-Q7 with a k-mer denoted by f. The procedures for the gSA depends on occ_C R (f), the occurrence number of f in the text made by the concatenation of all reads (i.e.i nC R ), while those for the Gk arrays depends on occ_Reads(f), the occurrence number of f in all reads, and we know that occ_Reads(f) ≤ occ_C R (f). Indeed, the gain of memory provided by Gk arrays increases with both k and q.I fk is small, each k-mer tends to occur more in average, and thus 4(r +1) ≪ 4mq,m e a n i n gt h a tGkCFPS is much smaller than the LCP array. If k is large then 4(m-k+1)q ≪ 4mq and thus, GkSA plus GkIFA tables occupy much less place than the gSA and ISA tables. This constitutes, in almost all cases, a saving of at least 12(k -1 ) q bytes. Locating a k-mer in the reads can be done with a binary search in O(k +l o gqm) worst case time with gSA using SA and LCP arrays and O(k × log q m) worst case time with Gk arrays using GkSA. (We recall that the binary search is the standard procedure in this context [START_REF] Manber | Suffix Arrays: A New Method for On-Line String Searches[END_REF][START_REF] Gusfield | Algorithms on Strings[END_REF]).

However, Manber and Myers [START_REF] Manber | Suffix Arrays: A New Method for On-Line String Searches[END_REF] mentioned that a simple improvement over the classical binary search (namely remembering the minimum length between the longest common prefix of the left and middle elements and the longest common prefix of the right and middle elements at each step of the binary search) permits to r u ni np r a c t i c ea sf a s ta saO(k +logq m) worst case method (see also [START_REF] Gusfield | Algorithms on Strings[END_REF] Section 7.14.3 page 152).

Thus, starting from a k-mer, rather than from a position, when answering the queries will bring an overhead similar in practice for the gSA and Gk arrays.

11 i i +1 ; 12 until (i ≥ qm) or (LCP[SA[i], SA[i + 1]] <k );
13 return (Ind k ); Nevertheless although we consider the same input, a position j of occurrence of the k-mer in a read, answering queries differ between the Gk arrays and gSA solutions. Indeed, since the gSA stores all positions in C R , we need to filter out positions of k-mers that overlap two reads in C R to keep only P-positions. This adds instructions to the procedure compared to that for the Gk arrays: see line 6 in Algorithm 5, which gives the algorithm for query Q1 with the gSA. For answering queries Q1 and Q2, we must perform another slight m o d i f i c a t i o n :w eu s eab i n a r ym a s kf o rd e a l i n gw i t h duplicate k-mers in a same read. This mask is stored in a binary vector B having q bits, one bit per read. The bit corresponding to a read is set to one whenever the kmer has been found to occur in that read, and subsequent occurrence positions in that read will be filtered out if the corresponding bit is set (lines 8 and 10 in Algorithm 5). Assume we query on a k-mer f from one of its occurrence position j. Let us denote by occ_C R (f) the number of occurrences of f in C R , including those overlapping two reads (i.e., starting at non P-positions), and by occ_-Reads(f) the number of its occurrences that are totally included in a read (i.e., those starting at P-positions). For Q1/Q2, Q5-Q7, we obtain with gSA a complexity of O(q + occ_C R (f)) since one initializes the bit vector B of size q and scan all occ_C R (f) occurrences. While with Gk arrays, the complexity depends linearly on occ_Reads(f) and we know that occ_Reads(f) ≤ occ_C R (f).

For Q3/Q4, there is no need of a bit vector with the gSA method, hence their complexity is O(occ_C R (f)), for one needs to scan positions in the gSA using the ISA and the LCP arrays. However, Gk arrays offer a complexity of O(occ_Reads(f))f o rQ 3a n dO(1) for Q4. We summarize all queries time complexities in Figure 3(c).

Remark 3. To avoid scanning occ_C R (f) entries, an alternative solution consists in delimiting reads inside C R using a separator. This solution would lead to a space overhead of q bytes for lowering the time complexity to occ_Reads(f). However we did not retain this solution since our goal is to diminish the space complexity and this solution would not improve much the time complexity.

A solution based on a hash table

An alternative solution is to index all k-mers in a hash table and to store for each read the list of its occurrence positions in the read collection. This list will contain pairs of integers: the read index in the collection, and the starting position of the k-mer in that read. The read index can be stored on a 32-bit integer, while a 16-bit integer suffices for the starting position. In such a case, storing the text is not necessary. The number of entries is the number of distinct k-mers in the read collection, i.e. our parameter r. Generally, r is small compared to 4 k for values of k in [START_REF] Salmela | Correction of sequencing errors in a mixed set of reads[END_REF]60]. Hence the hash table will be sparsely populated. We tried several implementation of state of the art hash tables: the Google sparse and dense hash arrays, and that from SGI extension of the C ++ Standard Library (called hash map).

Preliminary experiments have shown that Google sparse requires significantly much longer to build than SGI hash map, while having a lower memory footprint. With 20 million 75 bp reads, Google sparse hash occupies one third of the memory needed by the SGI hash map, but it takes thrice more time to build. On the contrary Google dense hash tables takes twice more memory, and offers only similar construction time. Hence, SGI extension hash map exhibited the best compromise in term of memory consumption and construction time compared to Google implementations. Thus, we choose SGI extension implementation for the comparison with Gk arrays.

Experimental settings

We tested index structures on three datasets.

1. We used a collection of 40 million Illumina ® RNA-Seq reads of length 75 from a human K562 library taken from the RGASP data (Accession number GM12878 at http://www.gencodegenes.org/rgasp with permission from B. Wold). We call it the K562 dataset. 2. We compiled several lanes of 454 ® genomic sequencing to obtain a collection of 2.8 million reads ranging from [40,3000] bp with an average read length of 523 bp. These were sequenced on a Roche 454 ® GS FLX platform with Titanium chemistry for the Khoisan genome project [START_REF] Schuster | Complete Khoisan and Bantu genomes from southern Africa[END_REF]. We call it the Khoisan dataset. 3. As much longer fixed length reads are not yet available, we constructed a collection of fixed length reads by slicing the Khoisan reads in non-overlapping pieces of 150 bp. We obtained 25 millions of 1 5 0b pr e a d s ,ar e a dl e n g t ht h a tw i l ls o o nb eg e n e rated on High Throughput Sequencing platforms.

k =15 k =20 k =25 k =30
In the first and third collections, reads have a fixed length, while in the second their length varies. The experiments were performed on an Intel Xeon 2.27 GHz equipped with 48 GB of main memory, and running Linux 2.6.18 with C++ compiled using gcc version 3.4.6 and -02 -funroll-loops options.

Experimental comparison

The use of read indexing raises three questions: how much computing resources does the index demand? Is it scalable? How fast can it answer large number of queries? Clearly the resources will depend on the number of reads (parameter q), their lengths and on the length of k-mers (parameter k). We compare three solutions: a hash table (HT), a generalized Suffix Array (gSA), and Gk arrays. Scalability We measured the construction time and amount of memory taken by all solutions for various numbers of reads and k-mer lengths. Figure 5(a) plots the maximal memory footprint on K562 data. At this scale, the value of k impacts only the hash table size; its influence on the gSA and Gk arrays is not visible on that graph. Second, the solutions can be ordered as follows: Gk arrays take the less memory, followed by the gSA, and then the hash table. This order is irrespective of the read number. For k =2 0e.g., Gk arrays use 10 GB, the gSA uses 20, and the HT 44, and the curves clearly indicate that these differences increase with the number of reads. Whatever the value of q, the hash table requires twice as much memory as the gSA, which itself takes at least 70% more memory than Gk arrays. With 25 million reads the hash table saturates the memory, with 30 million the gSA also does, while the Gk arrays constitute the only solution able to index the whole collection, 40 million reads, on that computer. Note that in both cases, the 64-bit implementation of gSA and Gk arrays have to be used to index that amount of reads. For the whole read collection, Gk arrays needs at most 43 GB (k = 15) and at least 36 GB (k = 30).

For all solutions, construction times increase linearly with the number of reads as expected (Figure 5(b)). It remains very similar between the gSA and Gk arrays, which both takes e.g. <1000 s. for 25 million reads. The influence of k is clearly visible on the hash table for 20 million reads: its construction time decreases with k because the parameter r also does (for a given number of reads). As long as they fit in memory, all compared solutions offer practical construction times.

We examined the behavior of Gk arrays on much longer reads, 150 bp, when variable-length read option is activated and when it is not. Figure 6(a) plots space consumption, while Figure 6(b) records the construction time for both options.

We see that adding a bit vector is not space consuming since there is little difference between the two methods (Figure 6(a)). For 13 million reads, the difference is, at most, of 300 MB between the two methods. In Figure 6(b), we plotted the construction time for both indexes. The variable length read implementation becomes slower when the number of reads grows, compared to the fixed length Gk arrays. This shows that despite a constant-time theoretical complexity for rank and select operations; there is a dependency on the length of the bit vector in practice. However, the construction time remains reasonable in the variable case.

Figures 7(a) and 7(b) plot space and time measured for the hash table and Gk arrays (with variable length reads option) on the Khoisan read collection. The gSA has not been implemented to handle variable length reads; note that the relative cost would have been similar to that observed with Gk arrays between fixed and variable read length options. Here for one million reads, variable length Gk arrays require 470 s. to build vs 428 s. for the hash table, but 8 times less memory (5.5 vs 46 GB). The difference increases strongly with the read number. Above one million reads, the memory footprint of the hash table exceeds the computer memory (which is 48 GB), while Gk arrays index the complete collection of 2.8 million reads on the same hardware with <15.6 GB. Hash tables appear to be more space consuming on the Khoisan dataset than on the K562 dataset. This can be explained by the nature of the data. Roche 454 ® sequencers offer a coverage depth much lower than Illumina's. Hence the number of distinct k-factors in the reads is likely to be greater with the Khoisan dataset. Answering queries We measured the mean time needed to answer 100,000 random queries of Q1-Q4. Since Q5-Q7 are slight variations of Q1-Q3 we do not report on these queries.

Figure 8 shows how the mean time for each solution vary with the number of indexed reads (q)a n dk on the K562 collection. Clearly, the influence of q is similar for all solutions, and small compared to the differences between solutions. Generally, gSA takes always longer than the hash table irrespective of the query type, and it also takes longer than Gk arrays for Q1-Q2 and Q4, and a similar time for Q3. The order between the hash table and Gk arrays depends on the query type. They are equally fast on Q1, the hash table does slightly better on Q2, clearly better on Q3, while Gk arrays is much faster on Q4. Anyway, for both the hash table and Gk arrays, the mean running time is in the order of or less than 10 microseconds for Q1-Q3, and around 0.1 microsecond for Q4 with the Gk arrays, meaning reasonable practical times.

For our comparison of Gk arrays with fixed or variable length read options, we see that the latter is becoming slower than the former (up to 7 times slower) when k is small, i.e. when the number of occurrences of k-mers is large. With larger k, the query time of the latter diminishes and becomes 2 to 3 times slower than with fixed Gk arrays.

With variable length reads (Figure 7(c)) the query times remain practical, but the hash table needs between 1 and 32 fold less time than Gk arrays depending on the query.

In summary, under various conditions Gk arrays are equivalent in construction time to a generalized Suffix Array or to a hash table. Compared to these solutions, they also offer reasonable query times under all circumstances; however, Gk arrays clearly outperform them in terms of memory footprint, the main bottleneck for processing High Throughput Sequencing data.

Conclusions

As High Throughput Sequencing becomes widespread, computational biology will face the challenge of managing astronomical quantities of short sequences. Mining such amount of sequences is feasible if the sequences are indexed in a preprocessing step. An index is a data structure that, like a telephone book, enables one to find easily a piece of information. For some value k,i t records the positions of all k-m e r si nt h er e a d si na n organized fashion to minimize the memory usage. Then finding the reads related to some k-mer takes as long as reading the k-mer and listing the corresponding reads, but not as long as scanning all the reads. In other words, read indexing factorizes the results of searches, which later speeds up the numerous queries made while the index is kept in memory. Our main contribution is to propose such an index: the Gk arrays. They are fast to build, require less space than alternative uncompressed solutions, and can thus handle larger read collections: 40 million vs 20 million reads for the hash tables with a memory limited to 48 GB. It is a key issue in practice.

While being comparable to hash tables in terms of time efficiency, only the Gk arrays can completely index a large read collection (like the K562 dataset) with a memory size available on nowadays computing servers. Moreover, our index remains fast for a wide range of values of parameter k (the length of k-mers). We have also shown that Gk arrays are both faster and smaller than an alternative generalized Suffix Array approach. Similarly, on variable-length reads like a Roche 454 ® dataset, Gk arrays can handle the whole read collection using less than 16 GB while hash tables are limited to a smaller sub-collection (about 1 million reads) on a 48 GB machine.

The Gk arrays answer efficiently different types of queries, but they have been optimised for queries where the searched k-mer is extracted from an indexed read. Sometimes one wishes to know for a given k-mer the reads in which it occurs and its positions inside those (e.g. assembly), while in other contexts one only wants the number of reads sharing this k-mer (e.g. estimation of expression level). Moreover, Gk arrays adapt well to variable length reads. Their scalability and versatility are key advantages, which allows to envisage multiple applications as mentioned in Introduction. However, scaling up to gigantic datasets (terabytes of data), as the ones obtained in large metagenomic projects, will require compressed read indexes. The simplicity of use of our index, and its implementation as a C++ library make it a software brick that can be easily exploited in future programs or further developed by the community.

For mapping reads on a reference sequence, solutions exist that index reads with hash tables [START_REF] Hach | mrsFAST: a cache-oblivious algorithm for short-read mapping[END_REF][START_REF] Li | A survey of sequence alignment algorithms for nextgeneration sequencing[END_REF]. For the error correction problem, other works have indexed reads with classical text indexing solutions: with a generalized suffix trie [START_REF] Salmela | Correction of sequencing errors in a mixed set of reads[END_REF][START_REF] Schröder | SHREC: a short-read error correction method[END_REF], a suffix array [START_REF] Ilie | HiTEC: accurate error correction in highthroughput sequencing data[END_REF], or hash tables [START_REF] Salmela | Querying large read collections in main memory: a versatile data structure[END_REF]. Gk arrays represent a first, attractive read indexing solution; it is specialised for this question and should suit different applications. Nevertheless, one can envisage several research perspectives. Indexing approximate k-mers or spaced seeds will authorize more types of queries, but will certainly increase the construction time and space requirements. Designing a dynamic construction algorithm for Gk arrays would futher enlarge their range of applications. Another challenge is to compress Gk arrays by storing sampled positions and recomputing other positions at run time, as done with the Burrows Wheeler transform [START_REF] Ferragina | Opportunistic data structures with applications[END_REF]. This would enable the user to adapt the index to its computer memory, while sacrificing some of its performance. Queries' running time comparison between the Suffix Array solution, the hash table and the Gk arrays. Queries' running time comparison between the Suffix Array solution, the hash table and the Gk arrays. Answering the queries Q1/Q2/Q3/Q4 on K562 dataset (75 bp reads). The plots represent the average time in μs over the same 100,000 queries of the corresponding type (i.e., Q1, Q2, Q3, or Q4). In all cases, the running time decreases when k increases for there are less occurrences in the reads of a say 30-mer than of a 15-mer. The Gk arrays are always faster than the Suffix Array; they even compute Q4 in constant time.

  and the set of positioned occurrences of f in all reads by Pos k (f)={ ( j,ℓ)|r j [ℓ..ℓ + k -1 ]=f}, where a positioned occurrence is given by the pair made of the read index in R and the beginning position of f in this read. Let us denote the restriction of Ind k (f)( r e s p .Pos k (f)) to subset of read indexes where f occurs only once by UInd k

  has lexicographic rank i among the P-suffixes.• For i a k-mer lexicographic rank and j a Q-position (i.e. i ∈ [0, r[ and j Q pos ),GkIFA[j]=i iff f Q (j)has lexicographic rank i among the P k -factors. • For i a k-mer lexicographic rank (i.e. i ∈ [0, r[), j 0 1234 5 6 7891 01 1 12 13 14 15 16 17 18 19 20 C R [ j] a acaa c t caa t t c a aacaag c SA[ j] 13 0 14 3 17 8 1 15 4 18 9 20 12 2 16 7 5 19 11 61 0 g( j)

Figure 1

 1 Figure1Example of the read index data structure: the Gk-arrays. Example of the read index data structure: the Gk arrays. Example for a collection R =( aacaact, caattca, aacaagc)o fq = 3 reads of length m = 7 and considering 3-mers (k = 3). The index is composed of three tables and uses a fourth one during construction (GkSA, GkIFA, GkCFA, and GkCFPS). The first table shows the starting indices of k-mers in the text made by the concatenation of all reads, C R , the SA built on C R , and the function g that renumbers P-positions of C R to make them consecutive. P-positions are {0,1,2,3,4,7,8,9,10,11,14,15,16,17,18}; all other positions, those starting positions where the k-mer overlaps two reads, are displayed with a gray background (lines j and SA[j]). Line SA refers to the usual Suffix Array of C R . The k-mer caa occurs 4 times in C R at positions 2, 7, 12 and 16. Among those, only 2, 7, and 16 are P-positions. The lexicographic rank of the P k -factor starting at position 16 is given by GkIFA[g[START_REF] Denoeud | Annotating genomes with massive-scale RNA sequencing[END_REF]] = GkIFA12 = 7, and the number of occurrences of the P k -factor caa is given by GkCFA7, which equals 3. The positions of these occurrences are thus obtained by the set {g -1 (GkSA[j]) | GkCFPS[7 -1] ≤ j < GkCFPS7} = {2,7,16}. See alsoFigure 2.

Figure 2

 2 Figure 2 Accessing the occurrences of a k-mer in the index. Accessing the index to get the occurrences of a k-mer starting at position j in the concatenation of the reads (i.e., C R ). Accessing GkIFA, GkCFPS and GkSA: (a) From C R to GkIFA: g(j) is the renumbered position of the Pposition j. (b) From GkIFA to GkCFPS: GkIFA[g(j)] is the lexicographic rank of the P k -factor starting at P-position j in C R , and GkCFPS[GkIFA[g(j)]] is the number of occurrences in C R of the P k -factors of rank less than GkIFA[g(j)]. (c) From GkCFPS to GkSA: The positions of the occurrences of the P kfactor starting at position j are in GkSA in the range [GkCFPS[GkIFA[g(j)] -1], GkCFPS[GkIFA[g(j)]]].

  [j] (by definition of GkIFA). The cardinality of Pos k (f)i s simply the number of occurrences starting at P-positions in C R ,w h i c hi sg i v e nb yGkCFA[t] (by definition of GkCFA). By Remark 1, GkCFA[t]=GkCFPS[t]-GkCFPS[t-1].Q3: Computing Pos k (f) By Proposition 3, all occurrences of f starting at P-positions are stored consecutively in GkSA. It suffices to find the lower and upper indices, denoted by ℓ f and u f respectively. By the ordering of GkSA all occurrences of factors smaller than f in the lexicographic order are stored before its occurrences in GkSA. Hence, by definition of GkCFPS and Proposition 2, we have

  ) and 4(b), respectively. For the gSA the three tables add up in memory and each takes 4mq space. With Gk arrays 1. the GkSA

arrays occupy 8

 8 mq +4(r +1) bytes (with 32-bit integers), and m := m -k + 1 is smaller than m.T h i se x p l a i n sw h y the memory footprint of Gk arrays remains smaller in practice than that of gSA (Figures4(a) and 4(b)), even for varying k values (see Figures 5(a), 4(a) and 4(b)).

Algorithm 5 : 2

 52 Q1 (Ind k (f)) with the generalised Suffix Array solution Data: f ∑ k , j P pos such that C R [j. .j+ k-1] = f Result: The set Ind k (f) 1 begin Ind k empty set; 3 Initialize the whole bit vector, D, to zero; 4 i ISA[j];//starting position of f occurrences in SA 5 repeat 6i f (SA[i]modm) ≤ m then //the occurrence position does not overlap two reads 7 readIndex SA[i]/m ; 8i f D[readIndex] ≠ 1 then //we have not found an occurrence in this read yet 9 Add readIndex to Ind k ; 10 D[readIndex] 1;

Figure 4

 4 Figure 4 Evolution of memory footprint during the construction of the generalized Suffix Array and of Gk arrays. Evolution of memory footprint during the construction of the generalized Suffix Array (a) and of Gk arrays (b) when indexing 15 million 75 bp reads with k = 25.

Figure 5

 5 Figure 5 Memory and construction time comparison between the Suffix Array solution, the hash table and the Gk arrays. Memory and construction time comparison between the generalized Suffix Array solution (gSA), the hash tables (HT) and the Gk arrays. K562 dataset is used for that experiment, with 5 million to 25 million reads. The length of k-mers ranges from 15 to 30. gSA plots have been shifted left and HT plots have been shifted right for easing the reading. (a) Maximal memory usage while constructing the index and querying it. The error bars represent the space consumption depending on the value of k. (b) Construction time for the three indexes on the same data as for the maximal memory consumption. The levels of gray on the plots represent the value of k.

Figure 6

 6 Figure 6 Comparing Gk arrays with fixed and variable length reads. Experiments on the sliced Khoisan dataset (150 bp-reads) with Gk arrays (fixed-length reads and variable-length reads). (a) Maximal space consumption of each index for several read numbers and values of k (the error bars represent the variation in space usage, depending on k). (b) Index construction time (the lighter gray corresponds to the smaller k). (c) Query computation time for 13,000,000 reads depending on different k values.

Figure 7

 7 Figure 7 Comparing hash table with Gk arrays on variable-length reads. Experiments on the Khoisan Roche 454 ® dataset (variable length reads) with a hash table and Gk arrays (with variable-length read option). (a) Maximal space consumption of each index for several read numbers and values of k (the error bars represent the variation in space usage, depending on k). The lower space consumption corresponds to a lower k. (b) Index construction time (the lighter gray corresponds to the smaller k). (c) Query computation time for 600,000 reads depending on the value of k.

1 µ 8

 18 Figure 8 Queries' running time comparison between the Suffix Array solution, the hash table and the Gk arrays. Queries' running time comparison between the Suffix Array solution, the hash table and the Gk arrays. Answering the queries Q1/Q2/Q3/Q4 on K562 dataset (75 bp reads). The plots represent the average time in μs over the same 100,000 queries of the corresponding type (i.e., Q1, Q2, Q3, or Q4). In all cases, the running time decreases when k increases for there are less occurrences in the reads of a say 30-mer than of a 15-mer. The Gk arrays are always faster than the Suffix Array; they even compute Q4 in constant time.

  table replaces gSA(C R )i nm e m o r ya n d takes only 4 mq, 2. GkIFA takes an additional 4 mq while GkCFA occupies 4(r +1) with r denoting the number of distinct P k -factors, and 3. finally the GkCFPS replaces GkCFA and takes exactly the same space.

				In total, gSA takes 12mq bytes of memory, while Gk
	generalized Construction	Memory	Cumulated memory footprint
	Suffix	time	building gSA(C R ) gSA built ISA built LCP built
	Array	O(mq)	4.02mq	4mq	8mq	12mq
			(a)
		Construction	Memory		Cumulated memory footprint
	Gk arrays	time building gSA(C O(mq) 4.02mq	4mq	8mq + 4(r + 1)	8mq + 4(r + 1)
				(b)

R ) GkSA built GkIFA & GkCFA built GkCFPS built

Philippe et al. BMC Bioinformatics 2011, 12:242

Acknowledgements and Funding

This work is supported by a CNRS PEPS grant "Bioinformatique, séquençage haut-débit et transcrits chimères en cancérologie", a CNRS PICS grant, a BioStic grant, the Region Languedoc Roussillon, and the ATGC bioinformatics platform. MS and NP were supported by fellowships from the French Ministry of Research, and NP benefits from a fellowship from the Ligue contre le cancer. ER thanks Dortmund University for the Dortmunder Gambrinus Fellowship. We gratefully thanks A. Mancheron for packaging the Gk arrays library.

Additional material

Additional File 1: Proof and queries' algorithms.

List of abbreviations used

Authors' contributions

All authors have designed the algorithm and contributed to the writing of the manuscript. NP and MS have developed the code. NP, MS, TL, ER have performed the experiments. ER supervised the manuscript redaction and submission. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.


For quite a long time, it has been claimed that the relationship between income and pollution was inverted U-shape, yielding the so-called Environmental Kuznets Curve (EKC hereafter). This claim, based on several early empirical studies (see for example, Grossman and Krueger 1993), has been at the heart of a massive empirical and theoretical literature. Empirical research has mainly consisted in examining a wide variety of pollutants for evidence of the inverted U-shaped pattern, resulting in the conclusion that such a shape is valid for many local and flow pollutants. But it does not seem to be the rule for stock pollutants like CO 2 which do rather generate monotonically increasing relation between wealth and pollution (see [START_REF] Brock | Economic Growth and the Environment: A Review of Theory and Empirics[END_REF], for a survey). Parallely, a great deal of theoretical contributions has been devoted to identify conditions under which the EKC arises. This includes optimal growth models like in [START_REF] Dinda | A theoretical basis for the environmental Kuznets curve[END_REF] and [START_REF] Stokey | Are there limits to growth[END_REF] or equilibrium models in the spirit of [START_REF] John | An overlapping generations model of growth and the environment[END_REF], all with stock pollutant.

In all the contributions mentioned above, the role of the abatement technology is crucial.

For example, in the latter abatement starts only after a large amount of capital, and thus pollution, is accumulated, ultimately generating the decreasing part of the EKC. A similar We shall study how the irreversibility mechanism interacts with the ingredients of the optimal technological switch problem outlined above, with a special attention to the outcomes regarding the capital-pollution relationship.

Our contributions are twofold. First of all, our contribution is technical. It is well known that including irreversibility in the above sense also induces an optimal switching time problem from the reversible to the irreversible ecological regime. Accordingly, our optimal AK growth model involves two optimal timing problems associated to technological and ecological switching times respectively. Original multi-stage optimal control techniques will be developed to solve the model, extending previous works of [START_REF] Tomiyama | Two-stage optimal control problems and optimality conditions[END_REF] on technological switching and Tahvonen and Withagen (1996) on ecological switching. Second, and more importantly, the interaction between the ecological and technological mechanisms generates a large set of potential optimal solutions. These solutions feature different relationships between capital and pollution. For our calibrated model, if a single technological switch is optimal, one recovers the EKC provided initial pollution is high enough.

If exceeding the ecological threshold is optimal, then the latter configuration is far from being the rule.

The problem

We consider an optimal growth AK model with two stocks, physical capital (K) and pollution (P). Two production technologies are available from t = 0. Each can be described by two parameters (A i , q i ), i = 1, 2: A i is marginal productivity of capital since the production function is Y = A i K, while q i measures the degree of dirtiness of the technology i.

Concretely, q i captures the marginal contribution of capital to the flow of pollution. The economy starts with technology i = 1 and has to decide whether it switches to technology i = 2, and when. The economy cares about both the levels of consumption and pollution. So

for the problem to be nontrivial, we shall assume that A 1 > A 2 but q 1 > q 2 : technology 1 is more productive but dirtier. A 1 -A 2 measures the cost of adopting a cleaner technology, we assume it permanent.The problem is at this stage similar to the one tackled in Boucekkine [START_REF] Boucekkine | Environmental quality versus economic performance: A dynamic game approach[END_REF]. There are however two major differences with respect to the latter contribution: the existence of capital accumulation and of pollution decay. We also introduce a key feature: irreversible pollution as in [START_REF] Prieur | Optimal emission-extraction policy in a world of scarcity and irreversibility[END_REF].

More precisely, the two state variables evolve according to the following laws of motion: K = (A iδ )K -C, for i = 1, 2, C being consumption and δ the depreciation rate, and Ṗ = q i Kα j P, with i, j = 1, 2 and α 1 > α 2 ≥ 0. The j subscript indexes the ecological regime, here parameterized by the pollution natural decay rate α j . While the technological switch from regime i = 1 to i = 2 does not require any minimal level of physical capital to take place, an ecological switch from j = 1 to j = 2, where α 1 > α 2 ≥ 0, does naturally entail the idea that nature cannot regenerate in the same way for low and high levels of the stock of pollution. Typically, this is modeled through a threshold value for the pollution stock, say P, above which the decay rate falls. 1

As explained in the introduction section, putting together technological and ecological switches enriches considerably the economic discussion. From the technical point of view, the problem sounds at first glance strongly asymmetric, the technological switch involving explicit timing decisions while ecological switches are essentially based on a threshold pollution level posited. It is however easy to see that this apparent asymmetry can be attenuated: in line with the pioneering contribution of Tahvonen and Withagen (1996), it is quite obvious to reformulate the ecological switching problem also as an optimal timing problem: indeed, if such a switch occurs, it exists a date, say t P , where P t P = P (assuming P 0 < P). As we shall see, this does not mean that the two state variables' laws of motion will imply similar optimality conditions. However, the previous observation legitimates the formulation of the two switchings problems both as optimal timing problems. Hereafter, t K with 0 ≤ t K ≤ ∞ will refer to the timing of technological switching while t P with 0 < t P ≤ ∞ does the job on the ecological side.

Given these provisos, our optimal control problem can be written as: max

{C,t K ,t P } V = ∞ 0 [U(C) -D(P)] e -ρt dt subject to, K =      (A 1 -δ )K -C if t ≤ t K (A 2 -δ )K -C else and, Ṗ =                  q 1 K -α 1 P if t ≤ min {t K ,t P } q 1 K -α 2 P if t P < t ≤ t K q 2 K -α 1 P if t K < t ≤ t P q 2 K -α 2 P if t > max {t K ,t P }
P 0 < P and K 0 are given, ρ > 0 is the rate of pure time preference. The control set includes the two timing variables mentioned above, plus consumption. The social welfare function is the same as in Tahvonen and Withagen (1996). The following standard regularity conditions assure the concavity of the problems we will have to handle along the way:

Hypothesis 1. The utility function, U(C), satisfies:

U(0) = 0, U (C) < 0, 0 < U (0) < ∞ and ∃! C/U ( C) = 0.
The damage function, D(P), satisfies: D(0) = 0, D (P) > 0, D (P) ≥ 0 and D (0) = 0.

The optimal control problem stated above is novel in that the two timing problems are of a different nature. In particular, one involves a threshold level for the state variable and the other no. Problems with multiple timing have already been tackled in the literature (see [START_REF] Saglam | Optimal pattern of technology adoptions under embodiment: A multi-stage optimal control approach[END_REF] but in the latter literature only technological switching is considered.

Here two types of timing problems are mixed in the same framework; needless to say, the interaction of both is very likely to give rise to a richer set of outcomes. Indeed, a quick inspection into the set of optimal regimes allowed in our enlarged problem is enough to get this point. A priori, one can list the following optimal regimes: 1. No switch, 2.

One technological switch, 3. One ecological switch, 4. Two switches: technological then ecological and 5. Two switches: ecological then technological. It is not obvious at all to guess a priori which kind of regime will result optimal given initial conditions, preferences and available technological and ecological menus. Even worse, one can identify within the optimal outcomes with at least one switch (the last five categories), different classes of solutions: interior (that is t K > 0 and/or t P > 0) or corner (t K = 0). Eight regimes are thus possible. Even more, one might be interested in distinguishing the case when there is a simultaneous ecological and technological interior switch, which adds another possibly interesting sub-case.

The next section gives our solution approach to this intricate problem.

The solution approach

We shall proceed as follows. First, for every possible regime k (k = 1 to 10 regimes from the discussion just above), we write the corresponding first-order necessary conditions and compute the resulting welfare function, say V k . Then, we pick the regime which delivers the largest social welfare, that is we identify the global maximum of the problem. The chal-lenge is first analytical because the general optimal control problem involved is nontrivial.

It is also computational because comparing nine possible regimes is highly demanding.

We shall ultimately resort to numerical comparison because our underlying optimal control problem does not admit a closed-form solution in most regimes. In the following, we provide the general control theory foundations to identify the solutions with two interior switches; this covers three cases: 0 < t K < t P < ∞, 0 < t P < t K < ∞, and 0 < t K = t P < ∞.

The other cases can be immediately recovered from the literature with one possible technological switch (see [START_REF] Boucekkine | Technology adoption under embodiment: a two-stage optimal control approach[END_REF] or one possible ecological switch (Tahvonen and Withagen, 1996), including corner regimes. 2

A natural approach is to decompose the problem into several sub-problems for given timing variables, to solve each of them, and then to identify the optimal timings. With one timing variable, two sub-problems are involved corresponding to the resulting two time intervals, before and after the switch. In our case, three would result from the occurrence of two switching times. In the spirit of [START_REF] Tomiyama | Two-stage optimal control problems and optimality conditions[END_REF], we shall use the following recursive scheme, illustrated here below on the case 0 < t K < t P < ∞:

• Third interval sub-problem: the problem in this regime is:

max {C} V 3 = ∞ t P [U(C) -D(P)] e -ρt dt subject to,      K = (A 2 -δ ) K -C Ṗ = q 2 K -α 2 P
where t P and the initial conditions K(t P ) and P(t P ) = P are fixed. The associated hamiltonian is:

H 3 = [U(C) -D(P)] e -ρt + λ 22 K ((A 2 -δ ) K -C) + λ 22 P (q 2 K -α 2 P), where λ i j
v is the co-state variable associated with the state variable v = K, P in the technological menu i and ecological regime j. The resulting value-function is of the form V * 3 (t P , K(t P )).

• Second interval sub-problem: in the next interval, the maximization problem is:

max

{C,t P ,K(t P )} V 2 = t P t K [U(C) -D(P)] e -ρt dt +V * 3 (t P , K(t P ))
subject to the corresponding dynamics, for regime i = 2 and j = 1, where t K , K(t K )

and P(t K ) are given, and t P and K(t P ) are free. Respectively denote by H 2 and

V * 2 (t K , K(t K ), P(t K ))
the corresponding hamiltonian and the resulting value-function.

• First interval sub-problem: This sub-problem considers the interval

[0 t K ]: max {C,t K ,K(t K ),P(t K )} V 1 = V = t K 0 [U(C) -D(P)] e -ρt dt +V * 2 (t K , K(t K ), P(t K ))
subject to the dynamics of regimes i = 1 and j = 1, with K(0) and P(0) given, and with free t K , K(t K ) and P(t K ). Again, denote the hamiltonian by H 1 and it is obvious that

V * 1 = V * . 3
Notice that each optimal control sub-problem is well-behaved, we will not spend space on writing the corresponding standard Pontryagin conditions. Rather, we will focus on uncovering the much trickier optimality conditions with respect to the timing variables and the so-called matching conditions. Matching conditions refer to how hamiltonians and the co-state variables behave at the optimal junction times. This is solved by the following theorem. 4

Theorem 1. Let 0 < t * K < t * P < ∞ be the optimal timing. Then:

H * 2 (t P ) = H * 3 (t P ) and H * 1 (t K ) = H * 2 (t K ), (1)
λ 22 * K (t P ) = λ 21 * K (t P ), λ 21 * K (t K ) = λ 11 * K (t K ) and λ 21 * P (t K ) = λ 11 * P (t K ). (2)
A few comments are in order here. First of all, one can read the five optimality conditions above are continuity or matching conditions at the junction times. In this respect, conditions (1) impose the continuity of the hamiltonian at the optimal junction times while the other conditions ensure the continuity of co-state variables at these times. Interestingly enough, one can observe that while at the technological switching time, both co-state variables are optimally continuous, only the one associated to K is necessarily continuous at the ecological switching time. This points at the major difference between the two switching types: in the latter, pollution is fixed at the switching time, equal to the threshold value, while at the technological switching time, both state variables can be freely chosen. This generally implies discontinuity of the co-state variable associated to P at t P .

Second, one can interpret the matching conditions (1)-( 2) as first-order optimal timing conditions for t P and t K respectively. Generally speaking, the matching condition for timing t i may be therefore written as: the marginal switching cost is nil. Therefore, the matching conditions on hamiltonians do equalize marginal benefits and costs of delaying switching times. Hence they do feature first-order necessary conditions with respect to the latters. 5

H * i (t v ) -H * i+1 (t v ) = 0, for i = 1, 2 and v = K, P.
There are two remaining cases of interest with two interior switchings, which derive quite trivially from the analysis of the benchmark case. If ecological switching precedes technological switching, that is if 0 < t P < t K < ∞, Theorem 1 still applies integrally. It is indeed invariant to the sign of t Kt P as one can infer from the discussion following Theorem 1. The last two (interior) switchings case follows the same logic though the list of corresponding first-order timing and matching conditions is shorter because only two successive regimes are involved, not three: one before t K = t P = t s and one after. Denoting by H i , λ ii K = λ i K and λ ii P = λ i P , i = 1, 2, the hamiltonians and co-state variables corresponding to the sub-problems on the intervals [0 t s ] and [t s ∞) respectively, optimality conditions

reduce to H * 1 (t s ) = H * 2 (t s ) and λ 1 * K (t s ) = λ 2 * K (t s ).
We end this section by assessing briefly the impact of technical and ecological switching on steady state. For any technological menu i = 1, 2 and ecological regime j = 1, 2, the expression of steady state capital and pollution are implicitly given by

U ((A i -δ )K i j ∞ ) = q i D ( q i K i j ∞ α j ) (A i -δ -ρ)(ρ + α j ) and P i j ∞ = q i α j K i j ∞ .
One can directly observe that both K i j ∞ and P i j

∞ are non monotonic in the technological menu (A i , q i ). In addition, if

K i j ∞ is decreasing in α j , P i j
∞ also is not monotonic in α j . Intuition runs as follows. Suppose the planner adopts, at some instant, the new technological menu. The reduction in A i generates the usual income and substitution effects that ambiguously affect investment. Moreover, the associated decrease in the intensity of pollution q i , induced by the adoption of the new technology, lowers the social cost of capital accumulation and stimulates investment. So, there are different forces at play and it is not immediately obvious which effect will prevail in the long run. In the same vein, the impact of technology adoption on pollution is unclear. Of course, the decrease in q i is a means to slow down pollution accumulation. But depending on whether technical change stimulates investment, there is an indirect effect that may go the other way round. In the numerical exercise to follow, we find it reasonable to assume that the adoption of the new cleaner technology allows the economy to reach a less polluted steady state. Regarding the impact of the ecological switch on steady state pollution, two effects are also pushing in opposite direction. Indeed, a decrease in α j makes it more difficult for Nature to regenerate itself.

At the same time however, everything else equal, it tends to reduce the incentive to invest in capital (through the higher social cost of pollution λ i j P ). From now on, we shall consider the more realistic case where the ecological switch translates into higher pollution.

It is now time to apply this theoretical analysis and our solution approach to the optimal growth model described in Section 2.

Numerical investigation

What should be the solution of the optimal growth problem with ecological and technical switches? To what extent does this solution respond to changes in the fundamentals of the economy? What is the relationship between the two state variables, capital and pollution, along the optimal path? These are the questions we shall address in this numerical analysis.

Model calibration

The analysis is conducted with the following utility and damage functions: U(C) = θC( C -C) and D(P) = νP, with θ , ν, C > 0. The set of baseline parameters used is

(3) θ = 14.4, q 1 = 0.04, q 2 = 0.02, z = 120, A 1 = 1/3, A 2 = 0.25, K 0 = 15, P 0 = 107, ρ = 0.05, δ = 0.075, α 1 = 0.005, α 2 = 0.002, c = 10, ν = 2.1
As for the discount rate, the chosen value is close to what most western governments use for most long term investments. The depreciation rate is usually between 5% and 10% depending on the level of economic development. We choose 7.5%. A review of the literature suggests that α 1 = 0.005 (implying a half-life of 139 years, see [START_REF] Hoel | Taxes versus quotas for a stock pollutant[END_REF]. The technology parameter value (A 1 ) is set in order to be consistent with the observed K 11 /Y 11 = 3. Moreover in our model, emissions before technology adoption are E 11 = q 1 K 11 . According to the DICE model (Nordhaus 2008, figure 5-10 p.110), E 11 /Y 11 = 0.12 for 2010, yielding q 1 = 0.04. At the steady state, ν = -λ 11 P∞ (α 1 + ρ).

According to the literature (see for instance Nordhaus 2008), the social damage costs of carbon dioxide emissions cannot be than less than 20$ or 30$ per ton but could probably not exceed 50$ per ton. Choosing 40$ to appraise the shadow price of pollution we obtain ν = 2.1. At the steady state again, θ = γλ 11 P∞ q 1 /(ρ + δ -A 1 ) where γ stands for the relative risk aversion. We assume the latter to be equal to 5 that leads to θ = 14.4. Parameter C in the utility function has been arbitrarily chosen because it is a scale parameter that does not affect conclusions of the simulations. Consistently with the discussion about the impact of ecological and technical switches on steady state, in the benchmark we choose parameters (A 2 , q 2 , α 2 ) so that steady state pollution is lower (resp. larger) after the technical (resp. ecological) switch than before. For the values reported in (3), it appears that steady state capital is higher after the technical switch than before. 6 So, if the adoption of the new tech-nology is costly in terms of investment in the short run, it is ultimately beneficial to capital accumulation.

The last important discussion refers to the set of initial conditions. For the particular issue of climate change, irreversible switches will likely occur if no severe action is undertaken soon. Indeed, there is now growing evidence that oceans (the most important carbon sink) display a buffering capacity that begins saturating. At the same time, the assimilation capacity of terrestrial ecosystems will likely peak by mid-century and then decline to become a net source of carbon by the end of the present century. Finally, the potential collapse of the North Atlantic meridional overturning circulation is drawing much of the attention, since it may happen for a 450 ppm CO 2 concentration while we have already reached 390 ppm [START_REF] Yohe | Reducing the risk of a collapse of the Atlantic thermohaline circulation[END_REF]). These considerations have led us to start the numerical exercise by setting the level of initial pollution close enough to -but below -the ecological threshold.

In addition, it is worth noting that we have chosen a threshold P < P 11 ∞ . It implies that a switch, either technical and/or ecological, will necessarily occur at some point in time. We further assume that a technological switch is a priori a means to avoid the ecological switch because P > P 21 ∞ . Finally, the initial endowment in capital K 0 satisfies

K 0 < K 11 ∞ (< K 21 ∞ ).

Benchmark scenario

Solving our calibrated AK model with pollution and switches, we first observe the existence of multiple optimality candidates. For the set of baseline parameters, we find solutions to the necessary optimality conditions associated with three of the nine possible regimes: 1.

Technological switch alone, 2. Simultaneous technological and ecological switches and 3.

Immediate technological switch alone.

Interestingly enough, as far as the nature of the relationship between K and P is concerned, these three candidates show very distinct features. Along the first regime, there is a sustained capital accumulation during the period before the switch, investment being more efficient than after t K . Consequently, K 21 ∞ is nearly reached at the time of the switch and pol-lution rises as well during this phase (see figure 1, top). At some point in time (t K = 15.61), the accumulated pollution and the level of capital become so high that the economy finds it worthwhile to adopt the new greener technology. From that point on, pollution, starting at a high level, decreases because natural assimilation now prevails on the lesser emissions due to the new technology. In the second regime, capital accumulation is delayed relative to what happens in the first regime in order to reach P as late as possible (see figure 1, second row). Once the double switch has occurred, pollution grows at a high rate because the economy cannot rely anymore on a high regeneration rate. Finally, the economy may choose to adopt immediately the new technology. In this third regime, we obtain a U-shape relationship between K and P (see Figure 1, third row): from the beginning, the economy is able with the new technology to accumulate capital while reducing the pollution stock.

What is the optimal regime? It turns out that it is the regime with a technical switch alone (W 1 = 2464.6 > W 2 = 2432.81 > W 3 = 1948.3). This result is very intuitive. Along regime 1, the planner can avoid the ecological threshold and is therefore better off than under regime 2. Regime 1 also dominates regime 3 because under the former, capital accumulation benefits from a higher investment efficiency. To conclude this analysis, it is worth mentioning that the optimal policy, with technical change alone, exhibits a capitalpollution relationship that has the feature of the EKC. 7

Let us investigate whether these results survive to modifications of critical parameters.

Sensitivity analysis

This section reviews all the possible variations around the baseline scenario. Table 1 summarizes our findings.

Several interesting conclusions hold whatever the scenario considered. Firstly, the regime with immediate technical switch always exists and is always dominated by another solution. Secondly, there is no solution featuring a technical switch followed by an ecological switch. Thirdly, each time a candidate with a technical switch alone exists, this yields the optimum. Last but not least, when the optimum is a technological switch alone, the relationship between capital exhibits a turning point and consequently is inverted U-shaped provided that initial capital and pollution are high enough 8 and that agent are sufficiently sensitive to pollution.

Finally, it is worth mentioning that under some scenarios, our conclusions sensibly differ.

In case where the intensity of pollution remains relatively high after adoption (q 2 = 0.03), the technological switch is not valuable to the planner. It implies that the optimum is the regime with a simultaneous double switch. For a not too damaging pollution (ν = 1), a small ecological threat (α 2 = 0.003), or more impatient agents (ρ = 0.1), less attention is paid to the ecological threshold and solutions to the regime with an ecological switch alone exist. However, it is optimal only in the case where agents are relatively impatient. In such a case, the capital-pollution relationship is not at all an EKC: it encompasses an intermediate stage where pollution increases whereas capital decreases (see figure 1, bottom).

Concluding remarks

This paper investigates the income-pollution relationship within and optimal AK growth model with technological and ecological switches. We show that the EKC, that is usually seen as a description of the relationship between wealth and pollution along the different development stages of one country, can also emerge as a result of the implementation of the optimal policy from the current development stage of the economy.

  scheme can give rise to the EKC under endogenous growth: the shift from insufficient to sufficient investment in abatement in an advanced stage of development curves down the pollution level at that stage. Another stream of the EKC literature puts forward technological progress: if richer countries are supposed to use, say, more energy-saving technologies,

FiguresFigure 1 .

 1 Figures

  Figure 1. Benchmark scenario: optimality candidates and corresponding paths. Last row: optimum with ecological switch (ρ = 0.1)

  This condition is quite common in the literature of multi-stage technological switching (see[START_REF] Saglam | Optimal pattern of technology adoptions under embodiment: A multi-stage optimal control approach[END_REF]. We show here that it applies also to ecological switches. Keeping the discussion non-technical, one may interpret the difference H *

i (t v ) -H * i+1 (t v ) as the marginal gain from extending the regime inherent to the time interval [t i-1 , t i ], with t 0 = 0, at the expense of the regime associated with interval [t i , t i+1 ]. Because there are no direct switching costs,

Table 1 .

 1 Summary of resultsTech. switchEco. switch t K < t P t P < t K Bench.t K = 15.6, V = 2464.6

	Tables

Notes

1 The case α 2 = 0 is the case for strong pollution irreversibility, we allow here for intermediate situations where the decay goes down above the threshold while nature keeps some degree of regeneration. 2 The case of two switchings, one interior and one corner, can be also settled similarly.

For example if 0 < t P < ∞ and t K = 0, one have just to solve for the ecological switch the particular problem with t K = 0 given and the implied capital law of motion.

3 One would use exactly the same scheme to handle a dynamic optimization problem in discrete time over three periods. Here the Bellman principle applies on the three intervals involved by the double timing problem instead of discrete periods of time.

4 A detailed proof of the theorem is available upon request. 5 The same type of arguments could be used to visualize easily the kind of necessary conditions implied by corner switching times: for example immediate technological switching,

6 Steady state values, for the benchmark, are: 366.3, 18.3) and (P 22 ∞ , K 22 ∞ ) = (272.9, 27.3).

7 Would A 2 be closer to A 1 (for instance A 2 = 0.3, see table 1) the investment efficiency gap before and after the switch would be reduced and the capital would be significantly increasing after the switch, therefore leading to a nicer EKC.

8 For a low K 0 , assimilation prevails on emission at the beginning, thus generating a decreasing relationship between K and P. For a low P 0 , the pollution keeps raising even after the adoption of the new technology.