Design of Worst-Case-Optimal Spaced Seeds

We introduce the different optimization problems that we consider in this work. For a given mask $\kappa$ of shape $(k,w)$ and given sequence length $n$, it asks how many changes $C^{\ast }=C^{\ast }(\kappa ,n)$ can be made at most in order to guarantee at least one unchanged $\kappa$-mer (for all possible placements of $c$ changes among $n$ positions). Equivalently, we may find the smallest $C^{†}$ such that $C^{†}$ changes, placed strategically, change all $\kappa$-mers; then $C^{\ast }=C^{†}-1$. Masks of the same shape $(k,w)$ may have different $C^{\ast }$-values (Figure 1), but even masks with the same (highest) $C^{\ast }$-value can be further distinguished.

This leads to another optimization problem for a given mask $\kappa$, sequence length $n$, and number of allowed changes $c$ (where $c\le C^{\ast }(\kappa ,n)$): Find a worst-case positioning of the changes among the $n$ sequence positions. There are two variants “MinHits” and “MinCov” of the problem that respectively ask to place the changes to minimize

1. 1.

   (MinHits) the number of (starting positions of) unaffected $\kappa$-mers (“hits”);

2. 2.

   (MinCov) the number of covered positions.

The corresponding decision problems, “does there exist a placement of $c$ changes such that there are zero hits (equivalent to zero covered positions)?”, has been proven to be NP-complete [29], but exact solutions can be obtained in practice by ILP solvers [14], as we demonstrate below.

Given a shape $(k,w)$, i.e., weight $k$ and window length $w$, we seek a (symmetric) mask $\kappa$ of the given shape that maximizes, among all such (symmetric) masks, the minimized worst-case number of hits (MaxMinHits) or covered positions (MaxMinCov). These are maximin problems, and the approach that we take to solve them is to solve MinHits and MinCov exhaustively for all masks of the given shape, and then pick the best mask.

We formulate the optimization problems MinHits and MinCov as integer linear programs (ILPs). The used integer variables are binary, i.e., take only values in $\{0,1\}$. We introduce the complete set of parameters and variables here, even though some of the problems use only a subset. The given (constant) parameters are:

• $\mathrm{■}$

  $n$: length of the sequence (typically $n=100$, corresponding to the length of a short read);

• $\mathrm{■}$

  $c$: number of allowed changes;

• $\mathrm{■}$

  $k$: weight of the mask, number of significant positions, typically $k\ge 21$ on a mammalian genome;

• $\mathrm{■}$

  $w$: window length, $k\le w\le n$

• $\mathrm{■}$

  $\kappa$: $k$-tuple of offsets with $0\in \kappa$, $w-1\in \kappa$, satisfying the symmetry condition $i\in \kappa \iff w-1-i\in \kappa$;

We use the following binary variables:

• $\mathrm{■}$

  $x={(x_i)}_{i\in [n]}$, indicators of changed positions: $x_i=1$ if and only if sequence position $i$ is changed;

• $\mathrm{■}$

  $y={(y_p)}_{p\in [n-w+1]}$, indicators of (starting positions of) unchanged $\kappa$-mers: $y_p=1$ if and only if the $\kappa$-mer starting at position $p$ is unchanged (i.e., none of its significant positions is changed),

• $\mathrm{■}$

  $z={(z_i)}_{i\in [n]}$, indicators of covered positions: $z_i=1$ if and only if sequence position $i$ is covered by an unchanged $\kappa$-mer.

As seen here, we use $i\in [n]$ for indexing sequence positions, $p\in [n-w+1]$ for indexing starting positions of $\kappa$-mers, and $j$ for an element of $\kappa$.

We point out a few examples of properties of contiguous vs. designed spaced masks with the same value of $k$.

For $n=100$ and (standard) $21$-mers, placing 4 changes at (0-based) positions $(20,41,62,83)$ affects all 21-mers, so the standard $k$-mer mask with $k=w=21$ (mask D in Table 2) only tolerates $C^{\ast }=3$ changes. In contrast, the $(k,w)=(21,25)$-shaped masks E and F in Table 2, #####_####_###_####_##### and ####_#####_###_#####_####, both tolerate 4 changes and then guarantee at least 8 hits and a minimum of 44 covered positions out of 100, for every distribution of 4 changes. It is noteworthy that the change distribution(s) that achieve(s) the worst case of 8 hits may be very different from the distribution(s) that achieve(s) 44 covered positions. In fact, one can show (with modified ILPs; not explained here) that both of these masks, when we have only 8 hits, we have at least 60 covered positions, and when we have only 44 covered positions, we have at least 20 hits. In any case, every distribution of 4 changes gives at least 8 hits and 44 covered positions.

Better still, both masks E and F tolerate even 5 changes, with at least 3 hits and 33 covered positions. There is no mask with $k=21$ and $21\le w\le 41$ that guarantees more than 8 hits for 4 changes, but the $(21,33)$-shaped mask G, ###_##_#__#_###_#_###_#__#_##_###, guarantees at least 55 covered positions but only 4 hits.

These examples illustrate the benefit of spaced seeds. Several other interesting masks are collected in Table 2, which are also used for our computational experiments. Full tables of optimal masks for different shapes $(k,w)$ for $c\in \{3,4,5,6,7\}$ changes on $n=100$ positions are included in the Supplement.

We used Gurobi 10.0.3 [14] to solve the ILPs from Table 1 on a AMD Ryzen 9 5950X 16-core Processor with hyperthreading (32 logical threads) and 128 GB RAM (which is not needed). Figure 2 shows an overview of solving times for all 1287 symmetric masks of shape $(19,29)$ for $c=5$ changes and sequence length $n=100$, and all 4368 symmetric masks of shape $(25,35)$ for $c=4$ changes (none of these masks tolerates 5 changes) and $n=100$. There are a few observations to be made.

First, solving times are very fast (under half a second) for the MinHits problem and reasonable (under 10 seconds) for the MinCov problem for the chosen parameters. Indeed, the MinHits problem is considerably smaller (Table 1), easier and faster to solve than the MinCov problem.

Second, the solving time depends on the optimal objective value. An objective value of zero (there exists a change distribution for which all $\kappa$-mers are changed) is found quickly. As a tendency, the better the mask (i.e., the higher the MinHits or MinCov objective value), the harder it is to find the corresponding worst-case change distribution or prove its optimality.

Third, unfortunately, increasing $c$ or $n$ can have a drastic effect on solving times. This should not be surprising: We are looking for the worst-case placement of $c$ changes among $n$ positions, for which there are $\left(\genfrac{}{}{0pt}{}{n}{c}\right)$ possibilities, where $\left(\genfrac{}{}{0pt}{}{100}{4}\right)\approx 4\cdot {10}^6$ and $\left(\genfrac{}{}{0pt}{}{100}{5}\right)\approx 75\cdot {10}^6$. While the time increase between 4 and 5 changes seems moderate in Fig. 2, solving for much larger values of $n$ and $c$ than provided in the Supplement is currently not feasible in a comprehensive manner. For example, solving the ILPs for the best $(19,29)$-mask ###_###_#__#_###_#__#_###_### for $n=100$ and $c=5$ (guaranteeing 5 hits or 48 covered positions) takes 0.28 seconds for MinHits, and 9.03 seconds for MinCov. Doubling these to $n=200$ and $c=10$ increases the time to 41 seconds for MinHits (a factor of 146) and to 5200 seconds for MinCov (a factor of 576).

However, we may argue that using $n=100$ (corresponding to a typical short read length) is sufficient because then the guarantees hold for every substring of length 100 of longer reads, and by the pigeonhole principle, if 50 changes are distributed over 1000 positions, there must exist at least one length-100 substring with at most 5 changes. In general, ILP solving times may depend on the decisions taken by the solver and may vary even for the same problem (if randomization is used) and across problems of comparable size and difficulty, and also change from version to version of the solver.

Clearly, established read mappers, such as bwa-mem2 [39], minimap2 [22, 23], or more recent developments like strobealign [37] work well in most settings; otherwise, they would not be used. We acknowledge that the optimal masks presented in Table 2 and their thresholds are mostly a theoretical result and perhaps of limited practical value for production pipelines in bioinformatics core facilities. Nonetheless, to investigate the case of highly diverse sequence regions, we decided to create datasets that pose a challenge to seed-based mappers (including our designed masks) because no long exact seeds exist.

We created two datasets, each of 5 million short reads of length $n=100$, from the autosomes of the t2t human reference genome [32]. The intervals of origin were chosen such that the reads, when unchanged, are easy to map back to the genome, in the sense that repetitive regions were excluded. Precisely, we ensured that all canonical $k$-mers of each selected read for $k=27$, are unique in the reference. In fact, we picked only regions where each 27-mer is not only unique, but also does not have a Hamming distance 1 neighbor elsewhere in the reference. Such $k$-mers have been called strongly unique [42]. Unique $k$-mers that are not strongly unique are weakly unique. As a rule of thumb, for $k\ge 23$, a large fraction of the $k$-mers in the t2t reference genome is strongly unique, and for $k\ge 27$, this fraction saturates [42], which motivated our choice of $k=27$.

Figure 3 contains statistics about how many positions in the t2t reference genome we find a strongly unique, weakly unique or non-unique $\kappa$-mer for selected masks from Table 2. Clearly, $k=19$ (A, B, C) does not yield many strongly unique $\kappa$-mers, but still an overall high fraction (over 60%) of unique $\kappa$-mers. For $k=23$ (I, J, K), we already have strongly unique $\kappa$-mers at over 60% of the genome positions, and a few more positions with weakly unique $\kappa$-mers. Increasing $k$ to $25$ (L, M, N, O) or to $27$ (Q) only slightly increases the fraction of strongly unique $\kappa$-mers further; the fraction of non-unique (“multi”) $\kappa$-mers stays more or less constant. We see that selecting reads from regions with only strongly unique 27-mers is not overly restrictive.

In the unchanged reads, any single contiguous 27-mer would be sufficient to reliably place the read at the correct location in the reference genome. While the unmodified reads would thus be trivial to map, we introduce a large number of changes ($c=5$ or $c=6$) to obtain the actual datasets.

In the first dataset, we introduced $c=5$ changes (A$\leftrightarrow$T, C$\leftrightarrow$G) into each read at approximately equidistant positions. More precisely, for each read, we modified each of the five 0-based positions $(16,33,50,67,84)$ by a random offset uniformly chosen from $\{-3,\mathrm{\dots },0,\mathrm{\dots },3\}$, and then complemented the base at the 5 obtained positions. In the second dataset, we proceeded similarly, but with $c=6$ changes at randomly modified positions obtained from $(14,28,43,57,72,86)$ with uniformly chosen random offsets between $-3$ and $+3$.

Clearly, this mapping task is challenging for all of the approaches: Using 5 or 6 changes makes it hard for the existing mappers to find a good seed, and none of our masks offers any guarantees for 6 changes either: All MinHits and MinCov objective values are zero for $c=6$. Note that some of our masks do offer guarantees for $c=5$.

To attempt to map the modified reads back to their interval of origin in the reference genome, we used bwa-mem2 v2.2.1 [39], minimap2 v2.26 [22, 23] and strobealign v0.11.0 [37] with their default parameters. For these tools, a read was considered uniquely mapped if the resulting BAM file contained a unique alignment; it was considered correctly mapped if that unique alignment started at the correct position (with a tolerance of half the read length).

For our selection of masks from Table 2, we did not create a read mapper for spaced $k$-mers within the scope of this work, but used a simple alignment-free approach. We created a positional index, implemented as a multi-way bucketed Cuckoo hash table [40, 41] for each mask, that for each unique canonical $\kappa$-mer stores its unique location (chromosome and position) and whether it is strongly (vs. weakly) unique. For non-unique $\kappa$-mers, it stores the special value non-unique.

To attempt to “map” a read with a mask, we query the genome position of each $\kappa$-mer in the read and, if unique, subtracted the relative position within the read, to obtain a virtual start position of the read in the genome. We count (for at most 20 different candidate starting positions) how many strongly and weakly unique $\kappa$-mers indicate each candidate position, and aggregate counts to the dominant position (highest count) inside an interval of half the read length. The read is considered uniquely mapped if we find at least two strongly unique or four unique $\kappa$-mers at the dominant position, and that position is the only remaining one with strong $\kappa$-mers after aggregation. The read is correctly mapped if the unique position is the correct one in the genome (with a tolerance of half the read length).

Figure 4 shows the fraction of correctly mapped reads for selected masks and each tool for both datasets. The alignment tools bwa-mem2 and minimap2 struggle particularly with the dataset using 6 changes, because they often do not find a sufficiently long seed to initiate alignment. For the dataset with 5 changes, bwa-mem2 still performs relatively well. Both minimap and strobealign suffer more because they do not examine all $k$-mers seeds but use a minimizer-based sampling approach.

Masks representing contiguous $k$-mers (A: 19-mers, L: 25-mers) also perform extremely badly because almost no unchanged 19-mer or 25-mer exists in these reads. While some 19-mers remain unchanged when only 5 changes are introduced (depending on the chosen random offsets), 19-mers are not specific enough to uniquely place a read within the genome using an alignment-free approach. Note that bwa2-mem also requires an exact match with a minimum length of 19 to initiate further alignment. Since bwa-mem2 also considers non-unique 19-mers and investigates up to 1000 positions, it is more successful than our simplistic approach that only uses unique 19-mers. The change-tolerant spaced 19-mers (B, C) perform much better than the contiguous 19-mer (A), but still not very well. The limiting factor here is the unavailability of sufficiently many (strongly) unique spaced 19-mers (Fig. 3).

Our masks with $k=23$ (I, J, K) perform almost equally well for both 5 and 6 changes, even though they guarantee only 2 hits for 5 changes and provide no hit guarantee for 6 changes. Because the change distribution is semi-random, we typically do not see the masks’ worst-case distributions and therefore usually obtain the required strongly unique $k$-mer hits even for 6 changes.

The spaced 25-mer masks (M, N, O) and the spaced 27-mer mask (Q) do not give any guarantees for 5 or 6 changes. Nonetheless, a large number of reads is placed correctly in both cases, again because the semi-random placement of changes does with high probability not correspond to the worst-case distribution.

In summary, we find that selected spaced masks with weight $k=23$ or $k=25$ are well suited to map reads with high substitution rates (5%, 6%) back to the correct location in genome, even in an alignment-free manner, whereas this task does not work well with contiguous $k$-mers (even for smaller $k$), which are used as seeds for most established tools. The low mapping rates of strobealign are surprising, but we note that it was designed mostly with indel tolerance in mind. Of course, the practical significance of this evaluation is limited: The datasets were designed to point out the difficulties of classical tools with certain distributions of large numbers of substitutions. Nonetheless, it is noteworthy that by using certain spaced masks as seeds instead of contiguous $k$-mers, one can achieve a high built-in substitution tolerance in alignment-free methods.

With this article, we provide the following tools and data, further explained in the repository’s README file¹¹1https://gitlab.com/rahmannlab/seed-optimization:

1. 1.

   a Python program that calls the Gurobi optimizer [14], for which a (free academic) license is required for computing the worst-case change distributions of a mask or a set of masks of the same shape,

2. 2.

   a just-in-time compiled Python program to index the $\kappa$-mers of a genome in FASTA format (gapmap index), which is a modification of our spaced $k$-mer counter hackgap [41],

3. 3.

   a just-in-time compiled Python program to map a FASTQ file of reads using a pre-computed $\kappa$-mer index (gapmap map).

   We furthermore provide the following data:

4. 4.

   a list of worst-case optimal masks for different shapes $(k,w)$ and number of changes $c$ for $n=100$ in the Supplement,

5. 5.

   our modified t2t reference genome (“analysis set”) consisting of the t2t genome [32], an Epstein-Barr virus (EBV) sequence and a PhiX sequence,²²2https://kingsx.cs.uni-saarland.de/index.php/s/CcggyrNZWpFoZWj,

6. 6.

   the challenge datasets of Sec. 3.3 of 100-bp reads with 5 changes³³3https://kingsx.cs.uni-saarland.de/index.php/s/nHJbaAxi3ZGApEG and 6 changes⁴⁴4https://kingsx.cs.uni-saarland.de/index.php/s/mA9cXMYcgHsZXFb. The FASTQ headers contain the true origins and the positions of the introduced changes.