Syncmers are more sensitive than minimizers for selecting conserved k‑mers in biological sequences

K-mers, submers and minimizers

Next-generation sequencing (NGS) has enabled dramatic advances in fields ranging from human functional genomics (Morozova & Marra, 2008) to microbial metagenomics (Gilbert & Dupont, 2011), generating vast amounts of data currently at petabase scale and increasing exponentially (Schmidt & Hildebrandt, 2017). Computational methods for analysis of NGS data are often based on k-mers, short subsequences of fixed length k. The number of k-mers in an index is often very large, motivating methods for selecting subsets as a time and space optimization. I use the term submers for a method designed to select a common subset of k-mers from similar sequences. The canonical example of submers is minimizers (Roberts et al., 2004). Minimizers are used in applications including read mapping (Li, 2018), taxonomy prediction (Wood, Lu & Langmead, 2019) and assembly (Sommer et al., 2007; Ye et al., 2012). Minimizers are defined by a choice of k, a window length w, and a coding function which assigns an integer value (code) to a sequence of length k. A coding function implies a total or partial order on the k-mers. The minimizer in a given window is a k-mer with smallest code value (ties may be handled in different ways), and the set of minimizers for a sequence is constructed by taking the union of minimizers over all windows. A simple and popular coding function is to consider letters in the nucleotide or amino acid alphabet to be base-4 or base-20 digits respectively (lexicographic coding). Often the coding function hashes the lexicographic code to improve the pseudo-randomness of the codes and/or to reduce the range of codes, for example, for hash table indexes, especially with larger values of k where a table of size 4^k or 20^k would exceed available memory. The read mapper Winnowmap (Jain et al., 2020) introduced the concept of weighted minimizers which assigns a weight to each k-mer biasing its probability of being selected in the event of ties. The density of a submer rule is defined as the fraction of k-mers selected in a long random string. For given k and w, the coding function can be optimized with the goal of reducing density (Zheng, Kingsford & Marçais, 2020; Marçais, Deblasio & Kingsford, 2018). Submers can also be defined by a universal hitting set (UHS) (Ekim, Berger & Orenstein, 2020; Orenstein et al., 2016, 2017), an explicitly enumerated set of k-mers such that any sequence of length w must contain at least one member. Minimizers and UHS submers provide a window guarantee, that is, there is necessarily at least one submer in every string of fixed length w. Here, I introduce syncmers, a family of alternative methods which select k-mers by inspecting the position of the smallest-valued substring of length s < k within the k-mer. For example, a closed syncmer selected if its smallest s-mer is at the start of end of the k-mer. At least one closed syncmer must be found in every window of length (k − s) k-mers (see Fig. 1 and proof in “Methods”). Technical terms are summarized in Table 1 for the reader’s convenience.

Submer conservation

Methods based on k-mers are generally used for comparison of similar sequences where differences may have been introduced by biological mutations and/or sequencing error. Detection of similarity by a k-mer requires that at least k letters and perhaps more are exactly conserved. Any submer is deleted if one of its letters is mutated. A minimizer is also deleted by a mutation in a position outside its k-mer but within its window if that mutation causes some other k-mer to have a smaller code. It is therefore useful to classify submers as context-free, that is, selected by a rule which considers only the k-mer sequence, and context-dependent, where the rule also considers flanking sequence. Minimizers are context-dependent, while UHS submers and syncmers are context-free. Context-free submers are generally better conserved under mutation providing there are not too many overlaps between neighboring submers.

Submer distance distribution

The length distribution of distances between consecutive submers gives insights into the advantages and disadvantages of a submer rule. An upper bound on the distance provides a window guarantee. Also desirable is that short distances have low frequencies because these correspond to submers with long overlaps which increase the rate of submer deletions under mutation. Ideally, the modal distance would have frequency 1.0 and all other distances would have zero frequency. This would uniformly tile the sequence and thereby provide a window guarantee, maximally suppress short distances, and maximize conservation. Of course, this is not achievable in practice, but provides a conceptual standard for comparison.

Coding function

Nucleotide k-mers are converted to integers by lexicographic coding. Extension to amino acid sequences is straightforward but not considered further in this report. Unless otherwise stated, lexicographic codes were hashed by murmur64 (https://en.wikipedia.org/wiki/MurmurHash).

Minimizers

Mincode submers

Modulo submers

Closed syncmers

Open syncmers

Circular syncmers

Circular syncmers consider a k-mer sequence to wrap around and thereby to contain k distinct s-mers rather than k − s + 1. Non-circular syncmers are described as linear if needed to distinguish these two types. For example, the 5-mer ACGTA contains 2-mers AC, CG, GT and TA if considered to be linear, plus AA if circular. With lexicographic coding, ACGTA is a (k = 5, s = 2) linear open syncmer because the first 2-mer AC is its smallest, but not a circular open syncmer because the smallest 2-mer (AA) is not the first. By increasing the number of s-mers, circularity enables higher compression for given k and s. For example, the compression factor is ~(k − s + 1) for linear open syncmers which increases to ~k for circular open syncmers (by similar reasoning to Eq. (3)). Circular open syncmers do not have a window guarantee, as shown by the following counter-example with k = 4 and s = 2. The repeating string ACCACCACCA… contains three distinct 4-mers ACCA, CCAC and CACC where the position of the smallest 2-mer is 4, 3 and 2 respectively. For circular closed syncmers, I do not have a proof of a window guarantee or a counter-example.

Offset parameter

Syncmer variants can be defined by introducing an offset parameter t specifying the one-based position of the s-mer to be considered. For example, an open syncmer with offset t = 2 has the smallest s-mer in its second position rather than the first. This has the effect of eliminating immediately consecutive (distance = 1) linear syncmers, as shown by the following argument. Suppose κ₁ is an open syncmer with t = 2, and append one more letter to give a new k-mer κ₂. The smallest s-mer in κ₁ (σ⁺) must be in its second position, and this is now the first s-mer in κ₂. Either σ⁺ is the smallest s-mer in κ₂, or the new s-mer at the end of κ₂ is the smallest; in both cases the smallest s-mer is not in its second position and κ₂ is therefore not a syncmer. If the syncmers are circular, consecutive syncmers are strongly suppressed but not completely eliminated. Similarly, for any value of t, distances < t are suppressed. If t > 1, the window guarantee is lost, as shown by the counter-example of a long string of the same letter (e.g., AAAAA….) which cannot contain a syncmer. However, this could be considered a desirable property in practice by avoiding seeds in low-complexity sequence.

Down-sampled syncmers

Prefix submers

Submer properties on random sequences

Submer properties including the compression factor (c) were measured on a random sequence R of length 10⁶ letters {ACGT}. The compression factor was measured as the number of k-mers divided by the number of submers.

Submer conservation with random mutations

Conservation was measured by introducing random substitutions into R at a given frequency, giving a mutated sequence R′. The mutation rate is expressed as percentage identity of R–R’, for example, frequency 0.2 is identity 80%. A submer is conserved if the k-mer at a given position is unchanged (no substitutions) and this k-mer is selected as a submer in both R and R’. The fraction of letters which are both indexed and conserved (denoted ConsPctid, for example, Cons80) is defined as the number of positions appearing in conserved submers divided by the length of R. If a given position appears in multiple conserved submers, it is counted only once. This definition of conservation is designed to quantify the fraction of the mutated sequence covered by alignment seeds.

Evaluation on whole-genome alignment

I implemented a pair-wise whole-genome aligner as follows. Submers and their positions in two genomes are identified and indexed in a hash table. For each pair of submers with the same sequence, one in each genome, an ungapped seed-and-extend alignment (Altschul et al., 1990) is constructed. The seed alignment implied by the matching submers is extended until (maximum score) − (current score) > X, where X is a fixed parameter, and the alignment with the maximum score is reported. After selecting submers, this algorithm has only three free parameters relative to the match score (1): mismatch score (−3), X (16) and minimum score to report an alignment (100). The minimum alignment score is set to a high value to suppress short and/or low-identity alignments where homology is ambiguous. Submer sensitivity was compared by measuring the aligned fraction (AF) of the input sequences obtained by different submer types. AF is calculated as (n₁/L₁ + n₂/L₂)/2 where n₁ is the number of bases in the first genome which appear in an alignment, L₁ is the length of the first genome, and similarly for the second genome. The value of AF ranges from 0 (no alignments) to 1 (all bases in both genomes appear in at least one alignment). Gapped alignments were not considered because alignment length is sensitive to gap penalty values, and the difference in AF between two submer types may therefore vary with gap penalty choice. With gapless alignments, if one submer type has substantially higher AF, it is reasonable to conclude that this submer is more sensitive to homologous seeds. If two submer types have similar AF, then the result is inconclusive because it cannot be ruled out that one submer type would find more low-scoring alignments that are in fact homologous.

Genome pairs

Whole-genome alignments were assessed using the bacterial assemblies shown in Table 2. Pair-wise average nucleotide identity (ANI) and AF were measured from MEGABLAST (Morgulis et al., 2008) alignments having E-values ≤ 10⁻⁹ and length ≥100. Gap-open and -extension penalties were set to high values to ensure that alignments were gapless. I manually selected pairs of genomes with a representative range of ANIs, preferring pairs with higher AF. These pairs were used to measure submer sensitivity in whole-genome alignment.

Representative parameters

To select parameters representative of those used in practice, I chose minimizers (k = 15, w = 10) and (k = 31, w = 15) used by minimap2 (Li, 2018) and Kraken v1 (Wood & Salzberg, 2014), respectively. Comparable syncmer parameters were identified as those giving compression factors equal or better than the minimizers with equal or better seed conservation at 90% identity (Cons90).

Comparison of syncmers with minimap2-like minimizers

Table 3 shows properties for minimap2-like minimizers with selected k = 15 syncmers that achieve better compression and/or better conservation. Notice for example that open syncmers with s = 9, t = 3 achieve compression of 7.0 vs. 5.5 for minimizers (27% lower density) with conservation of 0.312 vs. 0.301 (4% better).

Comparison of syncmers with Kraken v1-like minimizers

Table 4 shows properties for Kraken v1-like minimizers with selected k = 31 syncmers that achieve better compression and/or better conservation. For example, open syncmers with s = 31, t = 5 achieve compression of 11.0 vs. 8.5 for minimizers (29% lower density) with conservation of 0.081 vs. 0.077 (5% better).

Comparison of distance distributions

Figure 2 shows distance distributions for selected k = 8 submers to illustrate how the distribution changes with submer type and parameters. As noted in the “Introduction”, an ideal distribution would have modal frequency 1.0, but this is not possible in practice. A desirable feature is an upper bound w so that all distances >w have frequency zero; this is equivalent to the window guarantee. Also desirable is that short distances have low frequencies because these correspond to submers with long overlaps which are more likely to be deleted under mutations. With minimizers, all distances have approximately equal frequencies (see “Methods”), and short distances are therefore not suppressed. Open syncmers with offset t > 1 strongly suppress long distances and eliminate short distances, as expected (see “Methods”). Compare the second panel in row (A) with the second panel in row (G), which show that minimizers with w = 8 have compression 4.5 and conservation 0.47, while open syncmers with s = 3 and offset t = 2 have the same conservation with substantially better compression (5.9), which can be understood as a consequence of the more desirable distance distribution of the syncmers in addition to their context-independence.

Whole-genome alignment

Table 5 shows results for whole-genome alignment on 19 pairs of bacterial assemblies using (k = 15, w = 10) minimizers and (k = 15, s = 5) syncmers. On 15/19 of these pairs (79%), syncmers had higher AF (more aligned bases, indicating higher sensitivity), there was one tie, and on two pairs minimizers had 0.1% higher AF. As expected from the earlier results, syncmers show a strong tendency to have better sensitivity at lower identities. The AFs obtained with submers at low identities are much lower than MEGABLAST, indicating that many alignable segments of length ≥100 do not contain a conserved submer of either type.

Maximum distance under mutation

Table 6 shows maximum distances for minimap2-like and Kraken-like minimizers compared to selected syncmers tested on random sequences of length 1 M (comparable to a short bacterial genome) and 100 M (short mammalian chromosome) at identities of 95% and 90%. For k = 15, these are linear open syncmers with s = 9 and offset t = 4; for k = 31 they are linear open syncmers with s = 21, t = 5. These were chosen because they have both better compression and better conservation than minimizers (Tables 3 and 4). Neither of these syncmer types has a window guarantee. In identical sequences, the minimizers have a maximum distance (window size w) of 10 (k = 15) and 16 (k = 31), but under mutation the maximum distance between conserved minimizers is much greater, ranging from ~20w (k = 15, 1 Mb, 95%id) to ~300w (k = 31, 100 Mb, 90%id). Thus, the window guarantee is degraded by more than an order of magnitude even at a relatively modest 5% mutation rate in a relatively short sequence (1 Mb), getting worse as the identity is reduced and/or the sequence length is increased. As shown in Table 6, syncmers often have lower maximum distance under mutation despite the lack of a window guarantee, which can be understood as a consequence of their better conservation. This observation calls into question the practical relevance of the window guarantee (see Discussion).

Parameter sweep

Table S1 gives quality metrics for a wide range of submers, enabling programmers and end-users to choose parameters giving a preferred trade-off between compression factor and conservation.

Submer rules

The modulo submer method seems obvious and has probably been reinvented many times. It is a natural standard of comparison for more elaborate methods such as minimizers, and it is puzzling that no previous comparison has been published to the best of my knowledge. Possibly the fact that arbitrarily long gaps can occur, albeit with exponentially low probability, was considered a self-evident disadvantage, and naively this appears difficult to overcome without considering flanking sequence in a window. The syncmer method demonstrates that correlations between neighboring k-mers can be exploited to achieve a window guarantee with a context-free rule. Many variations on this technique are possible, allowing fine-tuning of characteristics such as the spacing distribution at a desired density. Perhaps surprisingly, it is possible to obtain a window guarantee even for submers that are sparse in the sense that some letters are not covered; this can be understood as additionally exploiting the finite (and typically small) number of distinct s-mers (see the “strictly decreasing series” proof for open syncmers). Introducing an offset gives up the window guarantee in favor of suppressing short distances (that is, long overlaps) between neighboring syncmers, which improves conservation. The improved conservation achieved by this technique tends to compensate for its lack of a strict window guarantee so that it may exhibit a shorter maximum distance under mutation than similar rules which do have a guarantee for identical sequences.

Context-free submers do not have edge bias

A k-mer at the start or end (edge) of a sequence appears in only one window, and thus has a probability 1/w of being a minimizer which is roughly half the probability 2/(w + 1) for a k-mer far enough from the edge that it appears in w windows (see Eq. (1)). This issue can be significant in practice for short sequences such as next-generation reads, for example, for identifying overlaps in assembly. This was noted in (Roberts et al., 2004), where it was proposed that the problem can be addressed by introducing a special case at sequence terminals they called end-minimizers. This issue does not arise with context-free submers such as syncmers.

Choosing submer parameters

An application must choose submer parameters such as k and w or s. Typically, these choices are driven by maximizing sensitivity and reducing computational resources, which are conflicting goals. The primary parameter is k, which gives a trade-off between sensitivity and specificity. Small k is more sensitive because mutations are less likely to delete a k-mer, but gives more seeds which increases processing time and may result in more false positive predictions. Larger k improves specificity (seeds are more likely to be homologous) and saves time (fewer seeds will be found in diverged sequences), but reduces sensitivity because more k-mers are deleted by mutations. As an example, bacterial genomes are ~10 Mb when both strands are included. If k ≤ 11, then the probability of finding any given k-mer in a bacterial genome is ~1 because 4¹¹ ≈ 4 Mb. Therefore, an application such as Kraken which requires submers to be unique within a clade must set k ≫ 11 to achieve acceptable specificity. By contrast, a bacterial whole-genome aligner might reasonably set k ~ 11 because the number of false positive seeds between one pair of genomes should be tolerably low. In the following, I will assume that k has been fixed based on such considerations and focus on the choice of submer rule and the rule’s parameters.

The window guarantee is a weak heuristic

The window guarantee is designed to ensure that similar subsequences have seeds, that is, conserved submers. This is a reasonable heuristic, but when mutations are introduced, the guarantee is necessarily lost for any type of submer. As illustrated by the examples in Fig. 3, submers, especially context-dependent submers such as minimizers, are deleted much more rapidly than letters as sequences mutate, with a large majority deleted at identities around 90% or less. Therefore, even if a submer rule has a window guarantee for identical sequences, its seeds will often have gaps ≫w when sequences have ≲95% identity (Table 6). If the primary goal of the submer method is to identify seeds in homologous sequences, then the window guarantee does not provide an upper bound on spacing; it is at best a weak heuristic and should not be regarded as an essential property. If a rule without a window guarantee conserves more seeds at a given density, then it is likely to give better sensitivity in many practical applications.

Density is not the appropriate optimization metric

Several recent papers have focused on minimizing the density of minimizers for given k and w; see (Zheng, Kingsford & Marçais, 2020) and references therein. This would be an appropriate optimization strategy if submers were used to find identical longer substrings in different sequences, but this is rarely the primary goal of an application and other methods are better suited to this task (e.g., Burrows–Wheeler indexes). In practice, most algorithms use minimizers to find homologous sequences which may have diverged due to mutations and/or sequencing error; for example, minimap2 and Kraken are algorithms of this type. Such applications generally aim to maximize sensitivity to diverged sequences under time and space constraints. If two minimizer rules with the same k and w have different density, then the rule with higher density will conserve more submers and is therefore more sensitive. Tuning minimizer density is thus similar to tuning the window length in that it gives a trade-off between sensitivity and index size. If lower minimizer density is needed to save memory or increase speed, then the window length could be increased, or the rule could be changed to one giving lower density. Either strategy could give better conservation, depending on the rules. For most applications, the appropriate optimization metric for submers is conservation, which should be maximized while keeping density fixed and allowing other parameters (e.g., w) to vary. Density should be held fixed because this parameter determines time and space resources for a given k. Maximizing conservation at a fixed density will tend to maximize the sensitivity achieved by an application with resources determined by the density. Conservation should generally be assessed by the number of conserved letters covered by a seed rather than the number of conserved submers to avoid double-counting when submers overlap.

Well-conserved submers

If it is accepted that the optimization goal is maximizing conservation for fixed k and density, then context-free submers have a natural advantage over minimizers because a mutation anywhere in the window can delete its minimizer, while context-free submers are deleted only when a mutation occurs within the k-mer. The advantage is greater at lower identities, because mutations in minimizer flanking sequences are more common, and also at lower densities (for a given minimizer rule) because lower densities require longer windows, which result in more deleted minimizers due to flanking mutations. The apparent advantage of the window guarantee is therefore a liability in practice for minimizers. If the desired density is high enough that it can be achieved by closed syncmers, then they offer an unambiguously superior alternative to minimizers for a wide range of densities and k values because they simultaneously achieve lower density and higher conservation (Tables 3 and 4; Table S1) while giving a similar window guarantee. Strictly, these results are established only for my own implementation of minimizers, but they almost certainly apply to minimizer variants with different coding functions or tie-breaking rules if submers with the same density are compared rather than the same window length. At lower densities, variants such as open syncmers with an offset can again simultaneously achieve lower densities and higher conservation compared to minimizers (Tables 3 and 4; Table S1). Some variants lack a window guarantee, but these syncmers are nevertheless clearly superior to minimizers unless a practically relevant advantage of the guarantee can be shown.

Connections between syncmers, UHSs and minimizers

To the best of my knowledge, all published practical algorithms based on minimizers have used either pseudo-random order (hashing) or lexicographic order; in particular minimap2 and Kraken use hashing. Several theoretical papers have described minimizer variants based on non-random orders; see (Zheng, Kingsford & Marçais, 2020) and references therein. In particular, a UHS can be used to construct minimizers by defining an order such that all k-mers in the UHS have lower code values than all other k-mers. By definition of a UHS, every window of length w will contain at least one k-mer from the UHS, and the minimizer rule will therefore select one member of the UHS from every window. A syncmer rule with a window guarantee defines a UHS, and an order defined by this UHS can be used to select minimizers. Using my terminology, the miniception algorithm (Zheng, Kingsford & Marçais, 2020) can be described as selecting minimizers via a UHS constructed from closed linear syncmers with offset 1. From the perspective of (Zheng, Kingsford & Marçais, 2020), lower density is preferred for fixed k and w, and miniception minimizers are therefore better than using the UHS directly as a submer rule because the minimizer density is lower while preserving the window guarantee (there may be more than one member of the UHS in a window, and the minimizer rule therefore selects fewer submers than using the UHS directly as the rule). From my perspective in this paper, the window guarantee is not important, minimizing density is the wrong objective, and syncmers are preferred because they give better conservation than minimizers. For future work, I would propose optimizing the k-mer order for mincode submers rather than minimizers, with the goal of maximizing conservation with fixed k and density. Free parameters are in the coding function (order) and code value threshold. This provides a fully general framework because any context-free submer rule can be re-formulated as mincode submers by assigning lower values to submers than other k-mers.