The number of k-mer matches between two DNA sequences as a function of k and applications to estimate phylogenetic distances

We study the number Nk of length-k word matches between pairs of evolutionarily related DNA sequences, as a function of k. We show that the Jukes-Cantor distance between two genome sequences—i.e. the number of substitutions per site that occurred since they evolved from their last common ancestor—can be estimated from the slope of a function F that depends on Nk and that is affine-linear within a certain range of k. Integers kmin and kmax can be calculated depending on the length of the input sequences, such that the slope of F in the relevant range can be estimated from the values F(kmin) and F(kmax). This approach can be generalized to so-called Spaced-word Matches (SpaM), where mismatches are allowed at positions specified by a user-defined binary pattern. Based on these theoretical results, we implemented a prototype software program for alignment-free sequence comparison called Slope-SpaM. Test runs on real and simulated sequence data show that Slope-SpaM can accurately estimate phylogenetic distances for distances up to around 0.5 substitutions per position. The statistical stability of our results is improved if spaced words are used instead of contiguous words. Unlike previous alignment-free methods that are based on the number of (spaced) word matches, Slope-SpaM produces accurate results, even if sequences share only local homologies.


Introduction
Phylogeny reconstruction is a fundamental task in computational biology [1]. Here, a basic step is to estimate pairwise evolutionary distances between protein or nucleic-acid sequences. Under the Jukes-Cantor model of evolution [2], the distance between two evolutionarily related DNA sequences can be defined as the number of nucleotide substitutions per site that have occurred since the two sequences have evolved from their last common ancestor. PLOS  In the present paper, we study the number N k of word or spaced-word matches between two DNA sequences, where k is the word length or the number of match positions of the underlying pattern for spaced words, respectively. Inspired by Bromberg's Slope Tree approach, we study the decay of N k as a function of k. More precisely, we define a function F(k) that depends on N k and that can be approximated-under a simple probabilistic model of DNA evolution, and for a certain range of k-by an affine-linear function of k. The distance between two sequences-defined as the number of substitutions per site-can be estimated from the slope of F in this range, we therefore call our approach Slope-SpaM, where SpaM stands for Spaced-Word Matches. In an earlier implementation, we calculated N k and F(k) for many values of k, in order to find the relevant affine-linear range of F and its slope therein [43]. In the present paper, we show how two values k min and k max can be calculated within this range such that the evolutionary distance between the compared sequences can be estimated from the values F(k min ) and F(k max ). Since it suffices to calculate the number N k of (spaced) word matches for only two values k = k min and k = k max , this new implementation is much faster than the previous version of Slope-SpaM.
Using simulated DNA sequences, we show that Slope-SpaM can accurately estimate phylogenetic distances. In contrast to other methods that are based on the number of (spaced) word matches, Slope-SpaM produces still accurate results if the compared sequences contain only local homologies. In addition, we applied Slope-SpaM to infer phylogenetic trees based on genome sequences from the benchmarking project AFproject [44].

Design and implementation
The number of k-mer matches as a function of k We are using standard notation from stringology as used, for example, in [45]. For a string or sequence S over an alphabet A, |S| denotes the length of S, and S(i) is the i-th character of S, 1 � i � |S|. S [i‥j] is the (contiguous) substring of S from i to j. We consider a pair of DNA sequences S 1 and S 2 that have evolved under the Jukes-Cantor substitution model [2] from some unknown ancestral sequence. That is, we assume that substitution rates are equal for all nucleotides and sequence positions, and that substitution events at different positions are independent of each other. For simplicity, we first assume that there are no insertions and deletions (indels). We call a pair of positions or k-mers from S 1 and S 2 , respectively, homologous if they go back to the same position or k-mer in the ancestral sequence. In our model, we have a nucleotide match probability p for homologous positions and a background match probability q for non-homologous nucleotides; the probability of two homologous k-mers to match exactly is therefore p k .
In our indel-free model, S 1 and S 2 must have the same length L = |S 1 | = |S 2 |, and positions i 1 and i 2 in S 1 and S 2 , respectively, are homologous if and only if i 1 = i 2 . Note that, under this model, a pair of k-mers from S 1 and S 2 is either homologous-if the k-mers start at the same respective positions in the sequences-or completely non-homologous, in the sense that none of the corresponding pairs of positions is homologous. For a pair of non-homologous k-mers from S 1 and S 2 , the probability of an exact match is q k .
Let the random variable X k be defined as the number of k-mer matches between S 1 and S 2 . More precisely, X k is defined as the number of pairs (i 1 , i 2 ) for which holds. There are (L − k + 1) possible homologous and (L − k + 1) � (L − k) possible background k-mer matches. The expected total number of k-mer matches is therefore In [38], we used this equation to estimate the match probability p for two observed sequences using a moment-based approach, by replacing the expected number E(X k ) by the empirical number N k of word matches or spaced-word matches, respectively. Although in Eq (1), an indel-free model is assumed, we could show in our previous paper, that this approach gives still reasonable estimates of p for sequences with insertions and deletions, as long as the sequences are globally related, i.e. as long as the insertions and deletions are small compared to the length of the sequences. It is clear, however, that this estimate will become inaccurate in the presence of large insertions and deletions, i.e. if sequences are only locally related. Herein, we propose a different approach to estimate evolutionary distances from the number N k of k-mer matches, by considering the decay of N k if k increases. For simplicity, we first consider an indel-free model as above. From Eq (1), we obtain which-for a suitable range of k-is an approximately affine-linear function of k with slope ln p. Substituting the expected value E(X k ) in the right-hand side of (2) with the corresponding empirical number N k of k-mer matches for two observed sequences, we define In principle, we can now estimate p as the exponential of the slope of F, as long as we restrict ourselves to a certain range of k. If k is too small, however, N k will be dominated by background word matches, if k is too large, no word matches will be found at all. The range where we can use the slope of the function F to estimate the match probability p is discussed below. The above considerations can be generalized to a model of DNA evolution with insertions and deletions and to sequences that share only local homologies. Let us consider two DNA sequences S 1 and S 2 of different lengths L 1 and L 2 , respectively, that have evolved from a common ancestor under the Jukes-Cantor model, this time with insertions and deletions, and possibly including additional non-homologous parts of the sequences. Note that, unlike with the indel-free model with global homology, it is now possible that a k-mer match involves homologous as well as background nucleotide matches. Instead of deriving exact equations like (1) and (2), we therefore make some simplifications and approximations.
We can decompose S 1 and S 2 into indel-free pairs of 'homologous' substrings that are separated by non-homologous segments of the sequences. Let L H be the-unknown-total length of the homologous substring pairs in each sequence. As above, we define the random variable X k as the number of k-mer matches between the two sequences, and p and q are, again, the homologous and background nucleotide match probabilities. For realistic data, we also have to take into account homologies between one sequence and the reverse complement of the other sequence. We therefore have twice as many background k-mer matches than in our above simplified model. All in all, the expected number of k-mer matches can be roughly estimated as and we obtain Similar as in the indel-free case, we define where N k is, again, the empirical number of k-mer matches. As above, we can estimate p as the exponential of the slope of F. An important point is that these considerations remain valid if L H is small compared to L 1 and L 2 , since L H appears only in an additive constant on the left-hand side of (5). Thus, while the absolute values F(k) do depend on the extent L H of the homology between S 1 and S 2 , the slope of F only depends on p, q, L 1 and L 2 , but not on L H .

The number of spaced-word matches
In many fields of biological sequence analysis, k-mers and k-mer matches have been replaced by spaced words or spaced-word matches, respectively. Let us consider a fixed word P over {0, 1} representing match positions ('1') and don't-care positions ('0'). We call such a word a pattern; the number of match positions in P is called its weight. In most applications, the first and the last symbol of a pattern P are assumed to be match positions. A spaced word with respect to P is defined as a string W over the alphabet {A, C, G, T, � } of the same length as P with W(i) = � if and only if P(i) = 0, i.e. if and only if i is a don't-care position of P. Here, ' � ' is interpreted as a wildcard symbol. We say that a spaced word W w.r.t P occurs in a sequence S at some position i if W(m) = S(i + m − 1) for all match positions m of P.
We say that there is a spaced-word match between sequences S 1 and S 2 at (i, j) if the same spaced word occurs at position i in S 1 and at position j in S 2 , see Fig 1 for an example. A spaced-word match can therefore be seen as a gap-free alignment of length |P| with matching characters at the match positions and possible mismatches a the don't-care positions. Spacedword matches or spaced seeds have been introduced in database searching as an alternative to exact k-mer matches [46][47][48]. The main advantage of spaced words compared to contiguous k-mers is the fact that results based on spaced words are statistically more stable than results based on k-mers [8,38,[49][50][51][52]. In short, neighbouring spaced-word matches are statistically less dependent from each other than neighbouring word matches; this is also the reason why spaced seeds lead to more sensitive results in database searching [53], compared with seeds that are defined as exact word matches as used, for example, in the original version of BLAST [54].
Quite obviously, approximations (4) and (5) remain valid if we define X k to be the number of spaced-word matches for a given pattern P of weight k, and we can generalize the definition of F(k) accordingly: if we consider a maximum pattern weight K and a given set of patterns {P k , 1 � k � K} where k is the weight of pattern P k , then we can define N k as the empirical number of spaced-word matches with respect to pattern P k between two observed DNA sequences. F(k) can then be defined exactly as in (6), and we can estimate the nucleotide match probability p as the exponential of the slope of F. The relevant affine-linear range of F For simplicity, let us first assume that the two compared sequences have the same length L, and that they are homologous over their entire length. If, for a factor α, we want there to be at least α times more homologous than background word matches, we have to require L � p k � α � 2 � L 2 � q k . We therefore obtain a lower bound for k as If, on the other hand, we want to have at least β expected word matches, we have to require L � p k � β, so k would be upper-bounded by Therefore, in order to ensure that the lower bound given by (7) is smaller than the upper bound in (8), we must require It follows that the match probability p must be large enough for our approach to work. More specifically, a simple transformation of (9) shows that we have to require As an example, with a sequence length of 1Gb, a background match probability q = 1/4 and parameters α = β = 10, we would need p > 0.57, corresponding to a Jukes-Cantor distance of around 0.64 substitutions per position in order to satisfy (9) and (10), respectively.
The inequalities derived in this subsection can be easily generalized to the situation where sequences share only a local region of homology of unknown length L H -as long as we can assume that the ratio between L H and the total sequence length is lower-bounded by some known constant.

Implementation
It is well known that the set of all k-mer matches between two sequences can be calculated efficiently by lexicographically sorting the list of all k-mers from both sequences, such that identical k-mers appear as contiguous blocks in this list. This standard procedure can be directly generalized to spaced-word matches.
To estimate the slope of F, one has to take into account that the values F(k) can be used only within a certain range of k, as explained above. If k is too small, N k and F(k) are dominated by the number of background (spaced) word matches. If k is too large, on the other hand, the number of homologous (spaced) word matches becomes small, and for N k < 2 � L 1 � L 2 � q k , the value F(k) is not even defined. For two input sequences, we therefore need to identify a range k min , . . ., k max in which F is approximately affine linear. To find such a range, we want to use inequalities (7) and (8).
There seems to be a certain difficulty in this approach, since the inequalities that we want to use to estimate the match probability p depend not only on the length of the sequences and the background match probability q, but also on the probability p that we want to estimate. It is easy to see, however, that the left-hand side of (7) is monotonically decreasing as a function of p, while the right-hand side of (8) is monotonically increasing. We can therefore calculate k min and k max for a fixed small probability p 0 and obtain values k min and k max that work for any match probability p > p 0 , since these values are guaranteed to be within the affine-linear region that we want to use to estimate p. For a suitable small match probability p 0 , we therefore define and and estimate the slope of F as The match probability p can then be estimated as the exponential of this slope, provided that p � p 0 holds. In our test runs, we obtained good results with parameter values α = β = 1 and with two different values for p 0 , namely of p 0 = 0.6 in (11) and p 0 = 0.53 in (12) and k max ¼ À ln L ln 0:53 where L is the average length of the compared sequences. With these values, we then estimate p asp with the function F defined as in (3). Finally, we apply the usual Jukes-Cantor correction to the estimated match probabilityp, in order to estimate the Jukes-Cantor distance between the sequences, e.g. the number of substitutions per position that have occurred since they diverged from their last common ancestor. As a numerical example, to illustrate how our approach works, we applied our approach to the genomes of Shigella dysenteriae 1 197 (4.44 Mb) and E. coli strain UTI89 (5.15 Mb). In Fig 2, F(k) is plotted against the word length k for contiguous words. By inserting the average sequence length L into (14) and (15), we obtain k min = 19 and k max = 24. With these values, we can calculate the slope of F in the relevant range as    (3), is plotted against the word length k for contiguous words. From the length of the sequences, the values of k min and k max are calculated with (14) and (15) as k min = 19 and k max = 24.

Distances between simulated DNA sequences
https://doi.org/10.1371/journal.pone.0228070.g002 The number of k-mer matches between two DNA sequences as a function of k length L = 100, 000 each, and we estimated their distances with Slope-SpaM. Here, we ran the program (a) with k-mers and (b) with spaced words based on randomly generated patterns with a probability of 0.5 for match positions. In Fig 3, the average estimated distances are plotted against the real distances; standard deviations are represented as error bars. Our estimates are fairly accurate for distances up to around 0.5 substitutions per site, but become statistically less stable for larger distances. We did the same test runs with simulated sequences of length L = 1, 000, 000; the results of these runs are shown in Fig 4. In both cases, the results were more accurate and statistically more stable when spaced words were used instead of contiguous k-mers. The number of k-mer matches between two DNA sequences as a function of k

Distances between sequences with local homologies
As shown theoretically in the section Design and Implementation, our distance estimator is still applicable if the compared sequences share only local regions of homology, in contrast to other approaches that estimate phylogenetic distances from the number of word or spacedword matches. To verify this empirically, we generated pairs of semi-artificial sequences consisting of local homologies, embedded in non-related random sequences. In each test run, the random sequences that we added to the homologous sequences had the same length in both sequences, but we carried out test runs with random sequences of varying length. We then estimated phylogenetic distances between these sequences with Slope-SpaM, as well as with Mash, Skmer and Spaced, three other alignment-free methods that are also based on the number of word or spaced-word matches. To see how these tools are affected by the presence of nonhomologous sequence regions, we compared the distance values estimated from the semi-artificial sequences to the distances estimated from the original sequences by the same tool.
To generate these sequence sets, we started with a set of 19 homologous genomic sequences from different strains of the bacterium Wolbachia from Gerth et al. [55]. These sequences contain 252 genes; the length of the sequences varies between 165kb and 175kb in the different strains. We then generated 9 additional data sets by adding unrelated i.i.d. random sequences of varying length to the left and to the right of each of the homologous genome sequences. This way, we generated 10 sets of sequences where the proportion of homology was 1.0 in the first set-the original sequences without added random sequences-, 0.9 in the second set, 0.8 in the third set, . . ., and 0.1 in the 10th set. Within each of these 10 sequence sets, we estimated all 19 2 � � ¼ 171 pairwise distances with the four programs that we evaluated. To find out to which extent these programs are affected by adding the non-related random sequences, we calculated for each data set and for each sequence pair the ratio between the estimated distance value and the distance estimated for the respective original, fully homologous sequence pair. The average ratio of these two distance values over all 171 sequence pairs is plotted against the proportion of homology in the semi-artificial sequences in Fig 5. As can be seen, distances calculated with Mash are heavily affected by including random sequences to the original homologous sequences. For the 10th data set, where the original homologous gene sequences account for only 10% of the generated sequences, and 90% of the sequences are unrelated random sequences, the distances calculated by Mash are, on average, more than 50 times as high as for the original gene sequences that are homologous over their entire length. Skmer and Spaced are also affected by the added random sequences, although to a lesser extent. By contrast, Slope-SpaM was hardly affected at all by the presence of the nonhomologous sequences.

Phylogenetic tree reconstruction
In addition to the above artificial and semi-artificial sequence pairs, we used sets of real genome sequences with known reference phylogenies from the as benchmark data for phylogeny reconstruction AFproject [44]. AFproject is a collaborative project to systematically benchmark and evaluate software tools for alignment-free sequence comparison in different application scenarios. The web page of the project provides a variety benchmark sequence sets, and the output of arbitrary alignment-free methods on these sequences can be uploaded to the server. This way, the quality of these methods can be evaluated and compared to a large number of other alignment-free methods, among them the state-of-the-art methods.
To evaluate Slope-SpaM, we downloaded all data sets from the AFproject server that are relevant for our study, namely five sets of full-genome sequences that are available in the categories genome-based phylogeny and horizontal gene transfer: (1) a set of 29 E.coli/Shigella genomes [28], (2) a set of 14 plant genomes [56], (3) a set of 25 mitochondrial genomes from different fish species of the suborder Labroidei [57], (4) another set of 27 E.coli/Shigella [58] and (5) a set of 8 genomes of different strains of Yersinia [59]. These data sets have been used previously by developers of alignment-free tools to benchmark their methods.
We ran Slope-SpaM on the above sets of genomes and uploaded the obtained distance matrices to the AFproject web server for evaluation. For these categories of benchmark data, the AFproject server calculates phylogenetic trees from the produced distance matrices using Neighbor Joining [60]. It then compares the obtained trees to trusted reference trees of the respective data sets under the normalized Robinson-Foulds (nRF) metric. That is, the standard Robinson-Foulds (RF) distance [61] that measures the dissimilarity between two tree topologies is divided by the maximal possible RF distance of two trees with the respective number of leaves; the nRF distances can therefore take values between 0 and 1. On the AFproject server, the benchmark results are ranked in order of increasing nRF distance, i.e. in order of decreasing quality. The results of our test runs with the AFproject benchmark sequences are shown in Table 1.
It should be mentioned that the performance of the compared methods on the set of eight Yersinia genomes is in stark contrast to their performance on the other four data sets. This has also been observed for other methods evaluated in AFproject: methods that performed well on , Slope-Spam (this paper) and Spaced [40] between semi-artificial sequences consisting of 252 homologous genes from 19 strains of Wolbachia, embedded in non-related random sequences. Ten data sets were generated by adding random sequences of different lengths. The x-axis is the proportion of the homologous sequence within the semi-artificial sequences, the y-axis is the ratio between the distance estimated from the semi-artificial sequences and the distances estimated from the original, homologous gene sequences, see the main text for more details. https://doi.org/10.1371/journal.pone.0228070.g005 The number of k-mer matches between two DNA sequences as a function of k other data sets, performed badly on Yersinia and vice versa. It may therefore be advisable to look at this data set and the reference phylogeny used in AFproject in more detail. Table 2 shows the run time of Slope-SpaM on three of the data sets from AFproject that we used in the previous subsection. We used the set of 25 mitochondrial genomes from different fish species (total size 412 kb), the set of 29 E. coli genomes (total size 138 Mb) and the set of 14 plant genomes (total size 4.58 Gb).

Availability and future directions
A number of alignment-free methods have been proposed in recent years that estimate the nucleotide match probability p in an (unknown) alignment of two DNA sequences from the number N k of word or k-mer matches for a fixed word length k. This can be done by comparing N k to the total number of words in the compared sequences [9] or-equivalently-to the length of the sequences [27]. A certain draw-back of these approaches is that they assume that the compared sequences are homologous to each other over their entire length. Unfortunately, these methods cannot be easily applied to sequences that share only local homologies with each other. If two sequences share only a limited region of homology and are unrelated outside Table 2. Runtime in seconds of Slope-SpaM (with spaced words) and six other alignment-free programs on three different sets of genomes that were used as benchmark data in AFproject [44]. On the largest data set, the set of plant genomes, co-phylog and kmacs were unable to produce results. this region, a given number N k of word matches could indicate a high match probability p in an alignment of those homologous regions-while the same number of word matches would correspond to a lower p if the homology would extend over the entire length of the sequences.
In the present paper, we introduced Slope-SpaM, an alternative approach to estimate the match probability p between two DNA sequences-and thereby their Jukes-Cantor distancefrom the number of word matches. The main difference between Slope-SpaM and most other word-based methods is that, instead of using only one single word length k, our program considers a function F such that the values F(k) depend on the number N k of word matches of length k. We showed in that there is a certain range of k where F is affine-linear, and the nucleotide match probability p in an alignment of the input sequences can be estimated from the slope of F. Further, we showed that one can calculate two values k min and k max within the relevant affine-linear range, such that it suffices to calculated F(k min ) and F(k max ) in order to calculate the slope of F in this range. The runtime of our program is therefore essentially determined by the time required to calculate the number of k-mer matches for only two values of k.
From the definition of F, it is easy to see that the slope of F, and therefore the estimated match probability p, are not affected if the compared sequences share only local homologies. This was confirmed in our test runs with semi-artificial sequences, obtained by adding nonrelated random sequences to truly homologous sequences. We generalized this approach to spaced-word matches, i.e. word matches with possible mismatches at pre-defined positions, specified by a binary pattern of match and don't-care positions.
To evaluate phylogenetic distances estimated by our approach, we used simulated DNA sequences with known Jukes-Cantor distances. For these data, we could show that distance estimates obtained with Slope-SpaM are highly accurate for distances up to around 0.5 substitutions per sequence position. By using spaced words instead of contiguous k-mers, we could improve the statistical stability of our results. To evaluate Slope-SpaM and other alignmentfree methods on phylogenetic tree reconstruction, we used all relevant data sets from the benchmark project AFproject. As in the original AFproject study [44], we applied Neighbor Joining to the distance matrices produced by the approaches that we evaluated, and we compared the resulting phylogenetic tree topologies to the reference topologies that are available from the AFproject web page. This is a common way of evaluating alignment-free methods. As shown in Table 1, the performance of Slope-SpaM on these data was in the medium range, compared to the other methods that have been evaluated in AFproject, if our program was used with spaced words. We obtained better tree topologies with spaced words than with contiguous words. In terms of topological accuracy of the produced trees, however, Slope-SpaM was outperformed by some of the top methods evaluated in the AFproject study.
An advantage of Slope-SpaM, compared to most other alignment-free methods, is its high speed. It may be possible to further improve the run time of the program by using sketching techniques [62], to reduce the number of (spaced-)word matches that are to be considered. Conversely, if the number of word matches is too small for short or distantly related sequences, spaced-word matches based on multiple patterns may help to increase the range of sequences where our method can be applied.
In the last few years, a number of alignment-free approaches have been proposed that are able to use unassembled short sequencing reads as input [9,28,35,37,63,64], as a basis for sequence clustering [34], for phylogenetic tree reconstruction [65] or for phylogenetic placement [64]. We are planning to evaluate systematically if Slope-SpaM can be used to find the position of unassembled reads in existing phylogenies, or to estimate phylogenetic distances between species based on sequencing reads instead of assembled genomes.
Since the approach implemented in Slope-SpaM is novel and rather different from existing alignment-free approaches, improvements may be possible by systematically analyzing the influence of the parameters of the program. In particular, there may be more sophisticated ways of finding suitable values of k min and k max to further improve the accuracy of Slope-SpaM. Also, in our study we used the Jukes-Cantor model, the simplest possible model for nucleotide substitutions. Using more sophisticated substitution models may also lead to improvements in future versions of the program.
The source code of our software is freely available through GitHub under the GNU GPLv3 license: https://github.com/burkhard-morgenstern/Slope-SpaM.