Mapping sequences by parts

Background: We present the N-map method, a pairwise and asymmetrical approach which allows us to compare sequences by taking into account evolutionary events that produce shuffled, reversed or repeated elements. Basically, the optimal N-map of a sequence s over a sequence t is the best way of partitioning the first sequence into N parts and placing them, possibly complementary reversed, over the second sequence in order to maximize the sum of their gapless alignment scores. Results: We introduce an algorithm computing an optimal N-map with time complexity O (|s| × |t| × N) using O (|s| × |t| × N) memory space. Among all the numbers of parts taken in a reasonable range, we select the value N for which the optimal N-map has the most significant score. To evaluate this significance, we study the empirical distributions of the scores of optimal N-maps and show that they can be approximated by normal distributions with a reasonable accuracy. We test the functionality of the approach over random sequences on which we apply artificial evolutionary events. Practical Application: The method is illustrated with four case studies of pairs of sequences involving non-standard evolutionary events.

which allows one to compare sequences with shuffled or inverted elements. The main idea of their work is to combine local and global alignments. During the first stage, the method selects conserved segments (local) and during the second stage, it chains an optimal subset of the pairs of segments previously selected (global). Many specialized approaches have been developed to model the specific evolutions "by blocks" of certain elements of sequences: minisatellites [5] or swaps in proteins sequences [6]. In the latter, the method is based on selection of common segments by local alignments scores and can be applied in a more general framework. The approach proposed in [3], mostly applied to more than two sequences, proceeds in similar manner in its first stage, then performs post-treatments and a graph representation of common elements of sequences.
For simplicity, we present the N-maps without taking into account inversions. Some hints will be given about how to extend the definitions and the algorithms in order to handle this type of evolution. Under this restriction, the (optimal) N-map of a sequence s over a sequence t is basically the way of cutting s into N parts that maximizes the sums of the scores of the gapless alignments of all the N parts against t. The gapless alignments can be local or global and so can be the N-map. This approach can be seen as a generalization of the "alignment with a fixed number of gaps" method initially introduced in [7] and recently studied in [8,9]. As this method, our approach is an attempt to avoid the introduction of some arbitrary costs on the transformations between sequences (like gap penalties in the case of alignment). For this purpose, we need a concrete way to determine the "best" number N of parts for mapping a sequence s over a sequence t. As in [9], we define this problem from a probabilistic point of view. Practically, we choose the number of parts leading to the most significant optimal score. The significance is empirically evaluated among pairs of independent identically distributed (iid) random sequences of same lengths and symbols distributions as s and t.
The rest of this paper is organized as follows. Section 1 is devoted to formal definitions and basic properties of the N-maps. We present the algorithms computing the optimal scores and correponding N-maps of a sequence s over a sequence t in Section 2. The algorithmic complexities of these computations are O (|s| × |t| × N) in time and O (|s| + |t| × N) in memory space. These complexities have an extra factor N with regard to the classical pairwise alignment algorithms. However typical values of interest of N are small compared to the lengths of the sequences: choosing a number of parts of the same order as the lengths of the sequences does not make any sense. The choice of the number of parts is discussed in Section 3, in which we investigate the distributions of the scores of the optimal N-maps of random sequences. In particular, empirical evidences lead us to approximate these distributions by normal ones and to measure the significance of optimal scores in terms of Z-values. The approach is evaluated in Section 4 by applying artificial evolutionary events over random sequences and by measuring the ability of the approach to retrieve the corresponding homologous segments. Section 5 shows four case studies of sequences (two pairs of proteins, a pair of DNA sequences of transposon elements and a pair of sequences of genes of microbial genomes) in which the homologies cannot be reported by a classic alignment. Finally in Section 6, we discuss the approach and present some research directions we plan to explore.
The sources of the software computing N-maps are available at [10]. We also provide additional utilities to estimate Z-values, represent N-maps as pictures (see Section 5), filter, merge and extract common segments.

Notations and Definitions
We consider sequences (or strings) over some finite alphabet of elements called letters or symbols. In practical applications, symbols can represent nucleotides, amino acids or genes. The elements of a sequence s are indexed from 1 to |s|, where |s| denotes the length of s, i.e. s = s 1 s 2 ... s |s| . For 1 ≤ i ≤ j ≤ |s|, the notation s [i, j] designates the substring s i s i+1 ... s j . We note the reverse sequence of s, i.e. = s |s| s |s| -1 ...s 1 . The set of all sequences of length l over is noted . Let s and t be two sequences. A pair of intervals of positions ([a, b], [c, d]) is a diagonal of (s, t) if 1 ≤ a ≤ b ≤ |s|, 1 ≤ c ≤ d ≤ |t| and b -a = d -c. The first (resp. the second) interval of a diagonal D of (s, t) will be designated as the s-interval (resp. the t-interval) of D. In order to avoid to deal specifically with some "pathological cases", we allow diagonals to be empty (of length 0).

Definition 1 Let s and t be two sequences. A N-map of s over t is a N-tuple of diagonals of
Without loss of generality, we assume in the following that the diagonals of a N-map of s over t are indexed according to the positions of their s-intervals. In particular, the first diagonal (resp. the last diagonal) is the one with the smallest (resp. the greatest) start position of s-interval. http://www.almob.org/content/2/1/11 A N-map of s over t is nothing but a peculiar type of map from a subset of positions of s to the set of positions of t.
In other words, it associates at most one position of t to a position of s; and none, one or several positions of s to a position of t (See Figure 1 or Figure 2 for dotplot representation).
A given classical alignment (which is also a map between positions) can be seen, for a certain positive integer N, as a N-map, both of s over t and of t over s. More precisely, an alignment with a fixed number K of gaps, like studied in [7][8][9], is a (K + 1)-map which, with notations of Definition 1, verifies the additional conditions: For 0 <N ≤ K ≤ |s| and a given N-map of s over t, there is at least one K-map defining the same map from positions of s to positions of t.
Let be a scoring scheme, i.e. a map from to ‫.ޒ‬ The score associated to a N-map following is: As in classical alignment methods, we will consider in the following only additive scoring schemes, i.e. defined from a substitution matrix π, for all lengthes l and all pairs of sequences (u, v) ∈ × , by: The score of an empty diagonal (l = 0) is 0.
The maximum of the scores over all the N-maps of s over t is noted (s, t): An optimal N-map of s over t is a N-map with score (s, t). By convention, a 0-map is the empty set and (s, t) = 0.
Depending on whether the substitution matrix contains negative values or not, the optimal N-map is said to be local or global. These concepts are used by analogy with the case of alignment. When the matrix contains only non-negative values (global case), a corresponding optimal N-map of s over t will attempt to associate each position of s with a position of t, as in a global alignment. When the matrix contains some negative values (local case), the optimal N-map will be reached by considering only subparts of s which lead to a positive contribution of the total score when associated with a segment of t, once more as in local alignment. Basically, a global N-map of s over t spans the entire length of s (except possibly some boundary positions) while a local N-map identifies N non-overlapping segments of s with maximum scores against t.
Some pathological situations could arise in the local case.
In particular there could be some positions i of s such that π [s i , t j ] is negative for all positions j of t. Without considering empty diagonals, (s, t) would be not always growing with N.  . . |t| r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r

Algorithms
Given two sequences s, t and a positive integer N, we address two problems: • Problem 1: computing the optimal scores (s, t) with K running from 1 to N.
• Problem 2: outputting the diagonals of an optimal Nmap.

Computing the optimal scores
Let Best [i, j, K] be the maximal score obtained by a K-map of s [1, i] [i, j, 0] to 0 for all integers i, j and K, we have the following recurrence relation: The correctness of this relation is straightforwardly proved by induction. Let us consider the maximum involved in the right part of this equation. It is equal to: • Best [i, j, K] , if the greatest score is obtained by incrementing the length of the last diagonal of an optimal K-map of s [1, i] over t ending at (i, j) -then the last diagonal will end at (i + 1, j + 1).
To compute the entries of Best referred to the index (i + 1), we only need to know the entries referred to the index i. Thus computing the optimal scores of all the K-maps of s over t, with K from 1 to N, can be done in O (|s| × |t| × N) time using O (|t| × N) memory space to store the dynamic programming variables. We can introduce now the formal algorithm Alg_1 which solves the problem of computing the optimal scores for global or local N-maps without inversion.
Algorithm Alg_1 takes as input two sequences s, t and a number of parts N and returns: The correctness of Algorithm Alg_1 is proved by induction over the positions of s. The time and memory space complexities are straightforwardly analyzed.
The variables B l and M d are not involved in the computation of the maximal scores of K-maps for 1 ≤ K ≤ N (they will be used by the algorithm in charge of outputting the diagonals). If we are only interested in solving Problem 1, these variables as well as the lines 7, 10, 14 and 19 can be deleted. Algorithm Alg_1 will still return the optimal scores of K-maps with 0 ≤ K ≤ N in the array M s .

Outputting the diagonals of an optimal N-map
Before presenting the formal algorithm, we need to introduce some additional notations and results about "dividing maps".
We say that a position m of s is inside a diagonal of a N-map Figure 2a). This notion excludes two cases: 1. when m is not contained by any diagonal (this is usual with local N-maps), 2. when a diagonal is exactly ending at m in its first interval.
We denote as the maximal score obtained by a K-map of s [i,|s|] over t [j,|t|] starting at (i, j), i.e. such that its

Lemma 1 Let s and t be two sequences, m a position of s and
N a positive integer. The optimal score of a N-map (s, t) is equal to the maximum of the two following quantities: There are two possibilities: either the position m is inside a diagonal K of Γ, or not. In the first case, there are a K-map Γ' ending at [m, c K + m -a K ] and a (N -K + 1)-map Γ" starting at posi- which implies . In the second case, let K be such that m = b K or b K <m <a K+1 . There is a K-map Γ' of s [1, m] over t and a

For all
is an optimal N-map of s over t.
We are now able to introduce the formal algorithm Alg_2 which solves the problem of outputting the diagonals of an optimal global N-map without inversion.
Algorithm Alg_2 takes as inputs two sequences s and t, two positions i and j bounding a substring of s, and a number of parts N. It outputs the diagonals of an optimal N-map of s [i, j] over t ordered according to their first intervals.   ] s tN

Correctness of Algorithm Alg_2
Let us consider Best and defined for s [i, j] as follows.
For all r such that i ≤ r ≤ j and v such that [1, v]  to compute recursively an optimal (K -1)-map of over t, and a (N -K)-map of over t.
If the maximum is reached in "Loop ", the variables D max and contain the two diagonals on both sides of position m of a N-map with score . The diagonal D max (resp. ) is possibly NULL -and not outputted -if m is smaller than the first position (resp. greater than the last position) of the N-map. Applying again Remark 1 leads to the correctness of the algorithm.

Time and space analysis of Algorithm Alg_2
Let us consider the recursion tree of an execution of Alg_2 which outputs an optimal N-map of s over t. The root of this tree is the initial call to Alg_2 with the parameters (s, 1, |s|, t, N), its two children are the two recursive calls in lines 36 and 38, and so on. The depth level of recursion of the initial call/root is 0. The depth level of another call is recursively defined as the incremented depth level of its direct ancestor. Before the two recursive calls at lines 36 and 38, the execution time of a call to Alg_2 with the parameters (s, i, j, t, N) is bounded by c × (j -i + 1) × |t| × N, for a constant c. Time is spent essentially in the two calls to Alg_1 at lines 5 and 6. The two recursive calls are done with the parameters (s, i, j L , t, N L ) and (s, i R , j, t, N R ) where: • (1) Let us remark that because of the possibly unbalanced repartition of N into N L and N R between the subcalls, the Master Theorem [11], generally used to evaluate complexity of divide and conquer algorithms, cannot be applied to prove the desired time complexity.
Since the initial call is done with the parameters (s, 1, |s|, t, N), the following assertions can be proved by induction over the depth level of recursion.
• From Inequalities (1), the length of the substring of s bounded by the two parameters "positions" in a call of depth level d is smaller than .
• From Inequality (2), the sum of the parameters "number of parts" of all the calls of depth level d is smaller than (N -d).
Thus, the total time spent at a level of recursion d is smaller than c × × |t| × N. By summing over all the possible levels (at most N levels), it comes that the total (including all the recursive subcalls) execution time of a call to Alg_2 with the parameters (s, 1, |s|, t, N) is smaller than 2 × c × |s| × |t| × N. This ends the time analysis.

Theorem 2 Algorithm Alg_2 outputs the diagonals of an optimal N-map of a sequence s over a sequence t in time O (|s| × |t| × N) using O (|s| + |t| × N) memory space.
The algorithm taking into account inversions follows the same general outline with additional and symmetrical dynamic programming variables for "reverse diagonals".
A similar idea can be used to compute an optimal alignment with a fixed number N of gaps in O (|s| × |t| × N) time complexity using O (|s| + |t| × N) memory space. It improves the "SANK_AL" algorithm described in [9], which needs O (|s| × |t| × N) memory space.

Choice of the number of parts
Given two sequences s and t, the score of an optimal Nmap of s over t increases with N. The maximum of the optimal scores is reached at most with N = |s| and the corresponding maps generally do not make sense. Some a priori knowledge could help us to decide whether the increase of the score between the K-and the (K + 1)-map deserves to consider an extra diagonal, for instance by introducing a penalty growing linearly with the number of parts.
Without such a priori knowledge, a natural choice is to consider the most significant optimal N-map: here the one which minimizes the probability of observing an optimal score greater than (s, t) between a pair of iid random sequences with the same lengths as s and t, and with the probabilities of symbols set to the frequencies observed over s and t. This choice needs to have informations about the probability distributions of the optimal scores of N-maps. Even if the problem could sound more homogeneous than the alignment case, we failed to derive an analytical approximation of this distribution. However, two cases are quite simple to check: N = 1 A 1-map is nothing but a gapless alignment of s and t. The distributions of the maximal scores were well studied in the local case and are known to converge to extreme value (EV) distributions [12]. where p x is the probability of the symbol x in s. The opti-mal score turns out to be a sum of |s| random variables of this type. Thanks to the Central Limit Theorem, its distribution converges with |s| to the normal distribution . If t is not fixed but random, this distribution becomes a mixture of with weights depending on the probabilities of sequences t.
With reasonable assumptions about the length and the probability distribution of t, we can neglect all the components of the mixture except the one which has distribution (µ|s|, σ 2 |s|) where and . Figure 3 shows the evolution of the empirical density functions of optimal scores of N-maps with N in the range of 1 to 15. In Figure 4, we can see that the empirical density function corresponding to N = 1 in the local case is well approximated by an extreme value distribution. As N increases, even for small values, the empirical distributions differ more and more from extreme value distributions and approach quickly normal distributions both in local and global cases. For a given N, the empirical distributions of global optimal scores are closer to the normal approximations than the ones of local optimal scores.
[ , ]   Empirical and approximated density functions of the optimal scores of local (a) and global (b) N-maps (N = 1, ..., 15 from left to right)  Even if the distribution of optimal scores is of unknown form for intermediate values of N, the empirical observations show that normal approximations fit well except for very small values. This fact leads us to measure significance of the score of an optimal N-map in terms of Z-values as in [9]. The estimated Z-value of an optimal score (s, t) is the number of standard deviations separating this score from the mean: where and denote respectively the mean and the standard deviation of the optimal scores estimated from a given number of trials of pairs of random sequences with the same lengths and the same frequencies of symbols as s and t. The higher the Z-value of (s, t), the lower the probability of observing a greater score in the normal approximation. So an optimal N-map of s over t with a higher Z-value will be consider more significant.
The Z-values must be taken with caution for small N -let us say less than 5 -because the corresponding probabilities are underestimated (the tails of the empirical distributions are heavier than the tails of the normal ones when N is smaller). This point is not crucial because we use Z-values to select a relevant number of part N rather than to assess an absolute significance of N-maps, but it could cause an underestimation of the "real" optimal number of parts. When the estimated most significant number of parts is small, it may be useful to check one or two next values.
Because it is time consuming and its accuracy is not rigorously evaluated, the way of estimating the significance of an optimal score is not fully satisfying. Analytical approximations of the distributions of the optimal scores should be pretty much better but they are beyond the scope of this article.

Evaluation
To evaluate the ability of the approach to retrieve segments of sequences related by evolution, we apply a given number of evolutionary events (mutations and shuffles) to random sequences and we measure the intersection between the homologies known from the artificial evolution and the ones reported by the most significant N-map.
More precisely, given a length L, a number of parts K and an identity proportion α, the protocol follows the steps below for a fixed number of trials.
1. Generate a reference random sequence s a iid with uniform probabilities over symbols of length L over an alphabet of 4 or 20 symbols ("random DNA" or "random protein").
2. Split s a into K equal parts and let s b be the sequence obtained by shuffling these parts with respect to the reverse permutation: (1, 2, ..., K) → (K, K -1, ..., 1). This step defines a reference K-map of s a over s b .

Let
Empirical and approximated density functions of the optimal scores of: a) local 1-maps, 5-maps and 15-maps, b) global 1-maps, 5-maps and 15-maps We do not apply insertion/deletion events over sequences in the protocol because the approach deals with this type of evolution exactly in the same way as the "split and shuffle" of Step 2. Figure 5 shows the evolution of the means of the intersection ratios for K = 2, 5, 10 and 15 parts, as functions of the identity proportion conserved in Step 3, over: • 500 random DNA sequences of length 500 (Figure 5a), • 500 random DNA sequences of length 250 (Figure 5b), • 500 random protein sequences of length 250 ( Figure  5c).
The error bars displayed in Figure 5 report the corresponding standard deviations.
The agreement of the results is perfect or almost perfect when the identity proportion is high. The identity proportion, the number of parts, the length of the sequences and the number of symbols in the alphabet affect the intersection ratio. This can be explained by the fact that the ability of the approach to associate a given segment of s a with its artificially evolved counterparts in depends on the probability of observing another segment in with a better score (here identity proportion). For instance in the case of DNA sequences the expected identity proportion of two segments is 0.25 under the random model used in the protocol. So it is not surprising to observe that the intersection ratio is 0 when the identity proportion artificially required in Step 3 is smaller than this value ( Figure  5a and Figure 5b). In the case of protein sequences, the expected identity proportion is 0.05. Eve tunes smaller than in the DNA case and we observe better results for small values of the identity proportion in Figure 5c. Clearly the identity proportion and the length of the sequence affect the probability of associating with the artificial counterpart. The role played by the number of parts is twofold. First, since it determines the length of the segments split in Step 2 of our protocol, it has a direct effect on the preceding probability. Second, it increases the number of boundaries and the possibility of an error when associating positions which are located at the beginning or at the end of the segments.
Evolution of the intersection ratio with the identity propor-tion for: a) random DNA sequences of length 500, b) random DNA sequences of length 250, c) random protein sequences of length 250 Figure 5 Evolution of the intersection ratio with the identity proportion for: a) random DNA sequences of length 500, b) random DNA sequences of length 250, c) random protein sequences of length 250.

Case studies
In the three first case studies, N-maps are represented as pictures where horizontal bold lines represent the sequences compared. The names of the sequences are specified over and under the lines. Each diagonal is represented as two boxes connected by an edge, where each box corresponds to a segment of one of the sequences. The height of the two boxes depends upon the score of the diagonal divided by its length (see Figures 6, 7 and 8). This type of graphical representation is also used in [1,6]. In the last case study we display N-maps as dotplots in order to make the results easily comparable with the ones of [13]. For convenience reasons, the scores represented in all the figures are normalized by being divided by the greatest entry of the substitution matrix.

Proteins 1
We begin with a case study from [3]. It compares SHK1 protein present in Dictyostelium (SwissProt ID Q9BI25) with ABL1 protein present in human (SwissProt ID ABL1_HUMAN). These proteins share two common domains which occur in a different order in each protein.
When comparing these sequences, the most significant optimal scores are obtained for: • N = 5 for both local maps of SHK1 over ABL1 and of ABL1 over SHK1 (Z-values respectively 36.67 and 36.79 - Figure 6a With the global approach, each part of the sequence to map is associated to the segment which maximizes its gapless alignment score in the second sequence, even if this a) An optimal local 5-map of CRK over NCK Figure 7 a) An optimal local 5-map of CRK over NCK. b) An optimal local 6-map of NCK over CRK. The N-maps are computed using BLOSUM62 substitution matrix. b a) An optimal local 5-map both of SHK1 over ABL1 and of ABL1 over SHK1 Figure 6 a) An optimal local 5-map both of SHK1 over ABL1 and of ABL1 over SHK1. b) An optimal global 9-map of SHK1 over ABL1. c) An optimal global 15-map of ABL1 over SHK1.
Local and global maps are computed using BLOSUM62 substitution matrix (made positive by adding a constant term in the global case). c score is small and does not correspond to a "real" homology. As a result, a most significant global N-maps contain a greater number of diagonals and look more confusing than with the local case. However, despite the fact that not all the homologies reported are relevant, the global Nmaps are interesting because they provide a more complete representation of the common elements. They allow us to consider diagonals formed by segments which are not homologous or not long enough to be selected in an optimal local N-map but which can be meaningful and can suggest an evolutionary history when taken in the whole context. Let us illustrate this point with the 15-map of ABL1 over SHK1 (Figure 6c). The diagonal with the ABL1-interval [1,27] (the first segment of ABL1 in Figure  6c) is too short to be selected in the most significant optimal local N-map but it can make sense when taking into account the larger diagonal with the ABL1-interval [113,198] (the third segment of ABL1 in Figure 6c) that follows it -not consecutively -in the two sequences. This could suggest the deletion or the insertion of the ABL1interval [28 -112] along the evolutionary history of this protein.
A simple solution to make a global N-map clearer is to select only diagonals with scores greater than a given threshold and/or long enough (see Figure 8 -case study 5.3). In particular, by considering only diagonals with average scores over a threshold in the two global N-maps of this case study, we would obtain pictures very similar to the local 5-map.
Finally, the common domains reported in [3]

Proteins 2
We compare two proteins sequences from [6]. A CRK like protein (SwissProt ID P46109) and a NCK adaptor protein (SwissProt ID P16333). This example is given to illustrate the way of pointing out repeated elements and we consider only local N-maps. The most significant optimal scores are obtained for: • N = 5 for local N-maps of CRK over NCK (Z-value 17.29 - Figure 7a), • N = 6 for local N-maps of NCK over CRK (Z-value 18.91 - Figure 7b).
In Figure 7 we can see that the most significant optimal local N-map of NCK over CRK has an extra diagonal with regard to the N-map of CRK over NCK. Apart from this extra diagonal, these two maps share almost the same diagonals set. There are only some small changes on their boundaries, essentially because the non-overlapping constraint of Definition 1 applies either to one or the other sequence. The extra diagonal is composed of a segment of NCK (positions 226-247), which was not included in the diagonals of the reciprocal map, and a segment of CRK which is also part of another diagonal formed with the a) An optimal local 5-map of P-element over P-repressor and reciprocally (BLAST substitution matrix) Figure 8 a) An optimal local 5-map of P-element over P-repressor and reciprocally (BLAST substitution matrix). b) An optimal global 15-map of P-element over P-repressor (Identity substitution matrix). c) An optimal global 24-map of P-repressor over P-element (Identity substitution matrix). Diagonals with average score (here identity proportion) smaller than 0.6 were removed from the global N-maps. c positions 35-58 of NCK. Since they are both homologous to a same segment of CRK, we have a clue that these two segments of NCK are repeated elements. Note that retrieving all the repeated common elements of two sequences needs generally to map one sequence over another and reciprocally to make sure that all the associations of segments are reported.

Transposons
We consider here DNA sequences of two transposons elements occurring in two species of Drosophila and studied in [14]: P-element (GenBank ID AY116625.1) and Prepressor (GenBank ID AF169142.2). We use the BLAST substitution matrix for nucleotides [15] for local N-maps and Identity matrix for global ones.
The most significant optimal scores of N-maps are obtained for: • N = 5 for local maps both of P-element over P-repressor (Z-value 425.83) and of P-repressor over P-element (Zvalue 420.24) corresponding to a same optimal 5-map (Figure 8a), • N = 19 for global maps of P-element over P-repressor (Zvalue 132.75 - Figure 8b), • N = 24 for global maps of P-repressor over P-element (Zvalue 119.96 - Figure 8c).
The corresponding maps are represented in Figure 8 in which we keep only the diagonals of the global N-maps with more than 60% of identity. As expected, filtering the diagonals according to their scores makes the pictures clearer and closer to the local one.
Once more, many diagonals are shared between these three N-maps with small variations in their-boundaries. The two global maps show an extra homologous region formed by several diagonals probably too short to be taken into account in the most significant local N-map.
In Figures 6, 7, and 8 we can remark series of diagonals composed of intervals of positions which seem contiguous and occur in the same order in the two sequences. They cannot be replaced by a unique diagonal because they are separated by small gaps (too small to appear at the scale of figures). In other words, N-maps computing acts over these positions like a classical alignment.

Microbial genomes
This case study illustrates how the approach can be applied to comparative genomics. We compare two microbial genomes: Chlamydia trachomatis (GenBank ID AE001273) and Chlamydophila pneumoniae (GenBank ID AE001363) studied in [13].
Each genome is represented by the sequence of its coding genes in the order they occur. Genomes of Chlamydia trachomatis and Chlamydophila pneumoniae contain respectively 895 and 1052 genes. A gene is identified with the sequence of amino acids of the corresponding protein.
Thus, there are as much different symbols as the total length of the two genomes (except the unlikely case where several genes share exactly the same sequence of amino acids).
We compare two sequences/genomes s and t of symbols/ genes which are themselves sequences of amino acids and we need to define a substitution score π between genes (actually this is only required between the genes of the first genome and the genes of the second one). For two sequences of amino acids p a and p b , we set π[p a , p b ] to the (highest) identity proportion of an alignment of p a and p b . As this substitution score is non-negative, we will consider global N-maps.
Because of the particular type of sequences studied here, the estimations of the empirical means and of the standard deviations of the Z-values are computed in a slightly different way from the one described in Section 3. To estimate the significance of a N-map score of a genome s over a genome t, we compute over a given number of trials, the empirical mean and the standard deviation of the optimal scores obtained by mapping a random shuffle of s, over t.
The empirical distributions of the optimal scores observed by shuffling the first genome depend a lot on the nature of the substitution scores between the genes of the first genome and the genes of the second one. But in non degenerated cases (when the substitution levels between genes are not all the same) we observe a behaviour close to the one described in Section 3. In this case study the most significant optimal scores of global N-maps are obtained for: • N = 94 for the global map of Chlamydia trachomatis over Chlamydophila pneumoniae (Z-value 289.61 - Figure 9a), • N = 97 for the global map of Chlamydophila pneumoniae over Chlamydia trachomatis (Z-value 265.23 - Figure 9b).
Because of the number of diagonals involved in the rearrangement, which is relatively complex and includes several inversions, we represent N-maps as dotplots (see Figure 9). The authors of [13] use this type of representation and show similar figures.
The N-map approach allows us to perform genomes comparison without the initial step of identification of clusters of orthologous genes which is generally a necessary (and sometimes a critical) stage before comparing genomes [16,17]. However, the N-map approach is different to methods such as sorting by reversals because it does not construct an evolutionary history (in the sense that it does not provide a sequence of evolutionary events transforming the genomes). It is rather a way to connect conserved segments and can be seen as an alternative to identify orthologous genes. The fact that two genes are associated in a N-map does not depend only on the level of homology between these genes, but also benefits from the levels of homology between their respective neighbourhoods.

Discussion and future work
Mapping sequences by parts is a simple and effective way to find out similarities between two sequences in the presence of evolutionary events that do not preserve their linear order. This first version was written in order to introduce the idea of "computational mapping of sequences" and needs some technical improvements and extensions such as dealing differently with the bounds of the parts which are mapped or distinguishing different costs for mutational events, to become more realistic from a biological point of view.
In the local case, the optimal N-map of s over t is close to the selection of the N gapless alignments with higher scores. So the results obtained with local N-maps should be generally close to the ones obtained by methods based on local (gapped or not) alignments [2,3,6]. The main difference stands in the non-overlapping constraint of Definition 1. From our point of view, the originality of the method actually makes sense with global N-maps. Strengths and weaknesses of global versus local N-maps are analogous to the alignment case. The local approach allows us to report only significant homologies. But a drawback is that the level of significance needs to be fixed a priori, generally by shifting the entries of the substitution matrix more or less negatively. On the other hand, with the global approach (a positive matrix), adding a same positive constant to all the entries of the substitution matrix leaves the resulting optimal N-maps unchanged. A first drawback is that even weak homologies are reported, but this is not a real problem since they can be easily filtered. A more serious concern is that a strong homology can be possibly diluted in a longer (but weaker) one.
The method can be extended in several directions. A first natural way is to allow gaps while mapping each part of the first sequence. Basically it can be done by extending the definition of diagonal to not constrain the lengths of the two segments to be equal and by defining the score of an "extended diagonal" as the alignment score (penalizing gaps) of its two segments. The algorithms computing the maximal scores and optimal N-maps with extended diagonal scores (computed with linear or affine gap pen-a) Dotplot representation of an optimal global 94-map of Chlamydia trachomatis over Chlamydophila pneumoniae Figure 9 a) Dotplot representation of an optimal global 94-map of Chlamydia trachomatis over Chlamydophila pneumoniae. b) Dotplot representation of an optimal global 97-map of Chlamydophila pneumoniae over Chlamydia trachomatis. Darker is a diagonal, higher is its length-normalized score. alties) are essentially the same as Alg_1 and Alg_2. In particular, their orders of time and memory space complexities do not change. In fact, the current implementation of the method provides an option to align parts with a linear gap penalty. Nevertheless, we presented here the method with the gapless case because it appears conceptually clearer and does not need any parameter such as a gap penalty (this parameter is critical for the distributions of the optimal scores and they appear more confusing in the gapped case).
Further in the same direction, an interesting possibility of extension is to associate different kinds of penalties for insertions/deletions, inversions and shuffling, and to compute the greatest score of a map of s over t according to a substitution matrix and these penalties. From an algorithmic point of view and with reasonable kinds of penalties, this can be done by Dynamic Programming equations analogous to the ones used in Alg_1. These equations could be directly applied to compute the best score and an optimal set of diagonals of a "penalized map" of s over t with complexity O (|s| × |t|). Setting the different values of penalties is a natural way to introduce biological knowledges in the approach but this needs a strong expertise in sequence analysis. We are interested in collaborations in this direction.