SimdMinimizers: Computing Random Minimizers, fast

We aim to solve the following problem as fast as possible: given a bitpacked representation of a sequence of ACGT DNA characters, compute the positions of all (canonical) random minimizers.

This work introduces a carefully optimized algorithm to compute the minimizers of a genomic sequence. Conceptually, it consists of two parts. First, we introduce a scalar algorithm (Section 3) that uses $O(n)$ time and $O(w)$ space. Most of this can be trivially parallelized to $L=8$ independent lanes with AVX2 or NEON SIMD instructions, and in Section 4 we specifically handle the input and output. The parts we discuss are:

1. 1.

   We apply ntHash [25, 16], a pseudo-random rolling hash function for $k$-mers (Section 3.2).

2. 2.

   We compute the sliding-window minima of the hashes using the two-stacks method [12, 32] (Section 3.3).

3. 3.

   We compute canonical minimizers based on refined minimizers [27] to decide the strand of each window (Section 3.4).

4. 4.

   We extend the scalar algorithm to SIMD by using it on $L$ chunks of the (bitpacked) input sequence in parallel (Section 4.1).

5. 5.

   Lastly, we collect and deduplicate the $L$ parallel streams into unique minimizer positions (Section 4.2).

Our method is $3.4\times$ (for large $w=19$) to $6.8\times$ (for small $w=5$) times faster than the fastest non-SIMD algorithm for computing forward minimizers. For canonical minimizers, we only compare against a simple implementation and find over $15\times$ speedup. As a result, we can compute the minimizers of a human genome in 5.2 seconds, and the canonical minimizers in 6.7 seconds. We also adapt our method to support generic plain-text ASCII input ($|\mathrm{\Sigma }|=256$), which is slightly (30%) slower due to the larger input characters.

A Rust implementation of our method is publicly available at https://github.com/rust-seq/simd-minimizers. The packed sequence representation, the splitting into chunks, and the parallel iteration over their bases is extracted to a separate library, packed-seq, available at https://github.com/rust-seq/packed-seq.

In this work, we assume that the input sequence is over the DNA alphabet $\mathrm{\Sigma }=\{\text{A},\text{C},\text{T},\text{G}\}$ and that each letter is encoded using two bits: $\text{A}=\text{00},\text{C}=\text{01},\text{T}=\text{10},$ $\text{G}=\text{11}$. This encoding can easily be obtained from the ASCII representation by applying a mask: (c >> 1) & 3. Additionally, we assume that the whole sequence is bitpacked using this 2-bit encoding, which can be done as a preprocessing step on the input if necessary. Non-ACGT characters have to be handled during this preprocessing as well, and could be skipped, converted to A, or the input sequence could be split at these points. We assume that the hardware is little-endian, and that the integer value of a sequence $x_0{x}_1{x}_2\mathrm{\dots }{x}_{k-1}$ is given by ${\sum }_{i=0}^{k-1}{x}_i\cdot 4^i$.

Given parameters $w$ and $k$, a window $W$ of length $\mathrm{\ell }=w+k-1$ contains $w$ consecutive $k$-mers. The minimizer of the window is the smallest $k$-mer in the window. For random minimizers, $k$-mers are ordered by a pseudo-random order, usually given by comparing hashes of the $k$-mers. In case of ties, the leftmost smallest $k$-mer is chosen.

Our goal is to compute the absolute position of the minimizer of every window $W$ in the input text. Since adjacent windows often have the same $k$-mer as minimizer, we only want each position to be listed once in the output.

Because DNA is double-stranded and most sequencing technologies do not distinguish these two strands, genomic sequences have an additional constraint: a sequence and its reverse-complement (the reversed sequence of complementary bases $\text{A}\leftrightarrow \text{T}$ and $\text{C}\leftrightarrow \text{G}$) should be considered identical. To satisfy this constraint, canonical minimizers should return the same set of $k$-mers regardless of the strandedness of the input. Specifically, if the canonical minimizer of a window $W$ is at position $p$, then the canonical minimizer of the reverse-complement ${\overline{W}}^r$ of $W$ should be at position $|W|-k-p=w-1-p$. In practice, canonical minimizers are often overlooked as an implementation detail and most existing methods simply compare canonical $k$-mers, computed as $x^c=\min (x,{\overline{x}}^r)$, which gives a weaker guarantee [24].

Computing the minimizers of a sequence generally involves a few steps, as shown in Figure 1. First, the $k$-mers are hashed. We do this using ntHash, a rolling hash (Section 3.2). Next, in each window of $w$ $k$-mer hashes, we must find the position of the leftmost $k$-mer with the smallest hash. For this, we use a method based on the two-stacks algorithm (Section 3.3). Lastly, many adjacent windows will have the same minimizer, and thus these positions must be deduplicated.

In this section, we present a scalar algorithm (that will be trivial to parallelize later by using SIMD instructions Section 4), and we also introduce a variant to compute canonical minimizers in Section 3.4.

The rolling hash takes constant time per character, and our sliding window algorithm will also take amortized constant time per character, so that the entire method runs in $O(n)$. It needs $O(w)$ space to store the hashes of the current window. With the SIMD optimizations, this improves to $O(n/L)$ time and $O(Lw)$ space, where $L$ is the number of SIMD lanes.

So far, our algorithms for iterating the input bases (Algorithm 1), hashing $k$-mers (Algorithm 2), and computing sliding window minima (Algorithm 4) are completely scalar. Now, we would like to use SIMD instructions to speed them up. With 256-bit AVX2 instructions, for example, we can process $L=8$ lanes of 32-bit values at a time. A first approach could be to use this to compute the hash of $L$ consecutive minimizer values at the same time, and then to compute the minimizer of the $L$ new windows all at once. Unfortunately, this is tricky due to the sequential nature of rolling hashing and sliding window minima.

Instead of processing consecutive $k$-mers in parallel, we choose to split the input sequence into $L$ equally long chunks that we process in parallel. This way, we compute one hash of each chunk in parallel, and then compute one minimizer position of a window of each chunk in parallel as well. This works, because the computations on the chunks are completely independent of each other. We let adjacent chunks overlap by $\mathrm{\ell }-1$ characters, so that each window is fully contained in exactly one chunk. The total number of windows may not be divisible by $L$. In that case, we round up the chunk length and return the number of elements that was added as padding, so that these can be removed later. For simplicity, we omit that case from the code snippets.

Table 1 shows the time usage for various incremental subsets of our method. To start, iterating the $8$ chunks of the input and summing all bases takes 0.15 ns/bp. Appending all u32x8 SIMD vectors containing the bases to a vector takes 0.30 ns/bp, indicating that writing to memory induces some overhead. Collecting the second $k-1$-delayed stream of characters that leave the $k$-mer (by adding them to the non-delayed stream) has no additional overhead. Computing ntHash only takes 0.02ns/bp extra. The sliding window computation nearly triples the total time. Collecting the minimizer positions to a linear vector (i.e., transposing matrices and writing output for each of the 8 chunks) again incurs 50% overhead, ad deduplicating them actually again adds some time.

Going back a step, using canonical ntHash instead of forward ntHash takes 0.14 ns/bp extra, and determining the canonical strand (via a third $\mathrm{\ell }$-delayed stream and counting GT bases) takes another 0.49 ns/bp. As before, collecting and deduplicating are slow and add around 0.70ns/bp.

In conclusion, we see that iterating the chunks of the input and determining minimizers is quite fast, but that a lot of time must then be spent to “deinterleave” the output into a linear stream. As can be expected, canonical minimizers are slower to compute than forward minimizers, but the overhead is less than $50\%$, which seems quite low given that the ntHash and sliding window minimum computation are duplicated and a canonical-strand computation is added.

We compare against the minimizer-iter crate (v1.2.1) [22], which implements a queue-based sliding window minimum using wyhash [34] and also supports canonical minimizers. For an additional comparison, we optimized an implementation of the remap method with ntHash based on a code snippet by Daniel Liu [20].

Results are in Table 2 and Figure 5. minimizer-iter takes around 26ns/bp for forward and 33ns/bp for canonical minimizers, and its runtime does not depend much on $w$ and $k$, because the popping from the queue is unpredictable regardless of $w$. Rescan starts out at 11.7 ns/bp for $w=5$ and gets significantly faster as $w$ increases, converging to around 5 ns/bp for $w\gg 100$. This is explained by the fact that rescan has a branch miss every time the current minimizer falls out of the window, which happens for roughly half the minimizers at a rate of $1/(w+1)$. Thus, as $w$ increases, the method becomes more predictable and branch misses go down. Our method, simd-minimizers, runs around 1.61 ns/bp for forward and 2.20 ns/bp for canonical minimizers when given packed input, and therefore is $3.4\times$ to $6.8\times$ faster than the rescan method.

In Figure 5, we see that simd-minimizers’ performance is mostly independent of $k$ and $w$ since it is mostly data-independent. Only for small $w\le 5$ it is slightly slower due to the larger number of minimizers and hence larger size of the output.

As we use SIMD with $8$ lanes, we could in theory expect up to $8\times$ speedup. In practice this is hard to reach because of constant overhead and because the overhead to work well with SIMD in the first place. In particular for large $w$, rescan benefits from very predictable and simple code and only outputs unique minimizer positions, making it very efficient. In SIMD, on the other hand, we use a data-independent algorithm, and output the minimizer position for every single window, which then has to be deduplicated. Thus, it is nice to see that even for large $w$, our method is over $3\times$ faster, despite this overhead.

Apart from taking bit-packed input, simd-minimizers also works on ASCII-encoded DNA sequences of ACTG characters directly, which are then packed into values $\{0,1,2,3\}$ for ntHash (in that order) during iteration. This is around 0.25ns/bp slower, mostly because of the larger size of the unpacked input.

The mulHash variant is around 0.20ns/bp slower again, but works for any ASCII input. Performance on 100 MB of the Pizza&Chili corpus [5] English⁸⁸8https://pizzachili.dcc.uchile.cl/texts/nlang/ and Sources⁹⁹9https://pizzachili.dcc.uchile.cl/texts/code/ datasets is nearly identical to performance on the random DNA shown in Table 2.

We also run simd-minimizers on the chromosomes of a human genome (T2T-CHM13v2.0¹⁰¹⁰10Available at https://github.com/marbl/CHM13 [26]), of total size 3.2 Gbp. Here, computing forward minimizers takes 5.19 seconds, and canonical minimizers takes 6.71 seconds for $(w,k)=(11,21)$, which corresponds 1.67 and 2.15 ns/bp, which is within a few percent of Table 2.

For $(w,k)=(11,21)$, we get a density of 0.173 for forward minimizers, and 0.167 for canonical minimizers, which are both close to the expected density of $2/(w+1)=1/6\approx 0.167$. Changing to $(w,k)=(19,19)$, we get density 0.098 for forward minimizers and 0.100 for canonical minimizers, both close to the expected density of $2/(19+1)=1/10=0.1$. Thus, we see that in practice, ntHash is a sufficiently random hash function to approach the expected density of random minimizers with a perfectly uniform random hash function.

We test throughput in a multithreaded setting, by using 6 threads to process the 25 chromosomes (22, X, Y, and mitochondrion) in parallel. This way, processing takes 0.97 s (forward) and 1.27 s (canonical) when the input data is already loaded into memory, showing slightly above $5\times$ speedup. This is just below $6\times$ speedup as the chromosomes don’t perfectly partition the data into 6 equal parts, so that some threads finish before others.

For applications, we recommend to simply call our library in parallel from multiple threads as needed.

We chose to use the data-independent two-stacks method as the core of our algorithm, that returns the minimizer position of every window and then requires deduplication. Given the promising performance of rescan for large $w$ in Table 2, an interesting alternative could be to go the opposite way and speed up the rescan step. This is particularly relevant when minimizers are sparse, as the number of branches may be small enough for linear scans to benefit from SIMD. Accelerating linear scans could also be useful for smaller inputs such as short reads, where splitting into 8 chunks may not be very efficient.

Another approach would be to use 512-bit AVX512 instructions and process 16 lanes in parallel. In theory that could be another $2\times$ faster, but in practice the collecting and deduplicating of values may become an even larger bottleneck. In simple experiments, without further profiling and optimizing the code, it is in fact $20\%$ slower.

Additionally, implementing low-density schemes like the open-closed mod-minimizer [7, 11] would be a valuable extension.