TahcoRoll: fast genomic signature profiling via thinned automaton and rolling hash

Thread-safe shared memory hashing

Jellyfish [22] exploits the compare-and-swap assembly instruction to update a memory location in a multi-threaded environment, and uses the “quotienting technique” and bit-packed data structure to reduce wasted memory. It also provides a function to count only a list of specific k-mers. Squeakr [30] builds an off-the-shelf data structure based on counting quotient filter (CQF). It maintains both global and local CQFs to facilitate updates of each thread. CHTKC [27] constructs a lock-free chaining hash table for multi-thread hash accesses.

Disk-based hashing

Disk-based hashing reduces memory usage with complementary disk space. In general, this method breaks k-mers into bins, and reserves them in files. Each bin is then placed into the memory for calculation. DSK [31] splits k-mers into bins using a specific hash function based on the targeted memory and disk space. MSPKmerCounter (MSPKC) [32] presents a novel technique, Minimum Substring Partitioning, to lower the memory consumption of storing k-mers. Recognizing the fact that consecutive k-mers in a read often share a shorter substring, these consecutive k-mers can be compressed and stored in one bin. It is suggested to index reads with an odd number k less than 64. KMC [33], KMC2 [34], and KMC3 [35] are serial developments of parallel counters. These methods scan reads one block at a time, and utilize several splitter threads to tackle these blocks. KMC2 leverages the concept of minimizer to further reduce disk usage. KMC3 speeds up the running time and optimizes the memory usage by taking a larger part of input data and better balancing the bin sizes.

Probabilistic hashing

To avoid calculating the counts of k-mers with sequencing errors, BFCounter [36] uses Bloom filter to identify all k-mers that are present more frequently than a threshold with a low false-positive rate. The algorithm examines read data in two passes to avoid reporting false-positive counts. khmer [24] uses a streaming-based probabilistic data structure, CountMin Sketch [11]. The algorithm is designed to conduct in-memory counting, and cannot tackle k larger than 32.

Suffix-arrays

Suffix-arrays show the potential of examining arbitrary k-mers without any restriction of k on a single scan. However, constructing a suffix-array on read data can be computationally expensive. Tallymer [37] is tailored to detect de novo repetitive elements ranging from 10 to 500 bp in the genome. The algorithm first constructs an enhanced suffix-array, and indexes k-mers for a fixed value of k. The indexing step needs to be repeated for different k’s.

MSBWT [38] is designed to consolidate raw reads via a multi-string variant of Burrows–Wheeler Transform (BWT). Instead of concatenating all reads and sorting, MSBWT establishes a BWT on each string and merges these multi-string BWTs through a small interleave array. The final structure supports a fast query of k-mers of arbitrary k.

Burst tries

An obvious shortcoming of tree-based and hash-based solutions is that they can suffer from enormous cache misses when working on massive datasets. KCMBT [25] uses a cache efficient burst trie to store compact k-mers. The trie stores k-mers that share the same prefix in the same container. When a container is full, k-mers are sorted and burst. An appropriate balance between the container size and the tree depth is necessary to prevent constant sorting and bursting. As a result, KCMBT requires to employ hundreds of trees.

Unfortunately, it is limited to handle k-mers with k less than 32.

Problem statement

We focus on counting the occurrences of a set of representative k-mers instead of all possible variable-length k-mers because of the following two reasons. First, the number of all variable-length k-mers in a DNA sequence is huge. More specifically, the lower bound of time complexity to investigate occurrences of all possible k-mers with different k is at least O(L ²) for a read length L. Second, some existing works [20, 21] focus on discovering relevant variable-length k-mers for various applications. Given a list of signatures, it is not necessary to profile all possible k-mers.

Suppose that P is the set of representative k-mers as signatures, where the length k is different across the set. Given a set of sequencing reads T , our goal is to profile the occurrences of each signature p ϵ  P . Note that the lengths of occurred patterns are shorter than the read length.

More formally, for each signature p ϵ  P , we aim to develop an efficient algorithm to compute the number of overall occurrences c _p is:

where | p | and | t | indicate the lengths of the signature p and the sequencing read t, and t[i … j] denotes the substring of t from the i-th to the j-th character. For each signature p ϵ  P , occurrences in the read set T are counted as a number c _p , which is the objective of signature profiling. Figure 1 shows an example with five signatures P = { ATT, GA, TTG, AGAT, TC } and two sequencing reads T = { AATTGAGAT, ATTGACATCG } .

Figure 1:

An example with five signatures and two sequencing reads in signature profiling. Each segment represents an occurrence of the corresponding signature in the read.

Framework overview

Figure 2 illustrates the overall framework of our proposed TahcoRoll that consists of the automaton construction phase for pre-processing and the read query phase for profiling. In the automaton construction phase, given a signature set for profiling, we first binarize original signatures into binarized signatures, thereby building the corresponding thinned automaton in the manner of Aho–Corasick automaton. In the read query phase, we traverse the thinned automaton and leverage the rolling hash technique to accelerate the verification process for profiling.

Figure 2:

An example with five signatures and two sequencing reads in signature profiling. Each segment represents an occurrence of the corresponding signature in the read.

Aho–Corasick automaton

To address signature profiling, an intuitive way is to reduce the task into multiple pattern matching [28] by mapping signatures onto patterns and each set of reads onto the input text. Multiple pattern matching algorithms find all occurrences in a read for each signature, so the profiling results can be obtained by aggregating these occurrences. We propose to apply the AC [29], one of the state-of-the-art approaches for multiple pattern matching, to profile signatures.

AC conducts the matching process along a trie that corresponds to patterns. Each node in AC has a failure link that allows fast transitions from one node to the other representing its longest possible suffix without backtracking. Informally, AC constructs a finite state machine (or an automaton) that resembles a trie and failure links. The pattern matching process can be treated as transitions between nodes in the automaton, and failure links provide efficient transitions between failed matches. Figure 3 shows an example of AC with five signatures. For example, the node of signature ACAT has a failure link to the node of AT. Given a sequencing read ATTTC to be profiled, AC will first match the signature ATT in the blue node. Then, it fails to match the third T and transits to the orange node that still has no child of T. After traveling along the failure link again to the yellow node, both the last two characters TC can proceed towards the orange and brown nodes that indicate a match of signature TTC.

Figure 3:

The automaton of AC with five signatures. Black solid links are trie links, and red dashed links are failure links. Colored nodes and thicker links are traversed while profiling the read ATTTG.

The construction of the automaton in AC with signatures p ϵ  P only requires a simple breadth-first search (BFS) with a linear time complexity O ( ∑ p ∈ P | p | ) . To profile signatures in reads O(t ϵ  T ), AC only needs to simulate transitions on the automaton, which also has a linear time complexity O ( ∑ t ∈ T | t | + ∑ p ∈ P c p ) . The space complexity of AC is also linear, O ( ∑ p ∈ P | p | ) , to maintain a node and a constant of links for each character. In theory, AC is a perfect fit for signature profiling.

Thinned automaton with binarized pattern matching

Even though we have shown the theoretical capability of AC for signature profiling, there are still some hurdles for AC in practice. One of the most critical issues is the memory usage when the number of signatures is huge. More specifically, each character in signatures can be referred to as a trie node, which provides plenty information and consumes a considerable amount of memory. For example, the Python implementation of AC requires more than 240 GB of memory to process 24 million signatures whose lengths range from 135 to 151. Especially for signatures with fewer and shorter common prefixes, nodes tend to have more child nodes. The greater width leads to an increase in memory usage.

To reduce both the number of nodes and the width of the automaton, we propose the thinned automaton with binarized pattern matching. More formally, each signature p [ 1 … | p | ] ∈ P is transformed into a binarized pattern p ′ [ 1 … | p | ] before being added into the automaton. The i-th character p′[i] of p′ is defined as follows:

Note that these four characters can be randomly divided into two groups. In practice, we suggest grouping characters so that the constructed trie could be as balanced as possible based on data distribution. Besides, if the data distribution is too complicated to have a suitable grouping method for characters, a limitation of our proposed method is the inevitable effect on the efficiency. In this paper, we use a balanced partition which groups A and G together. Compressing two characters into one bit 0 or 1, binarized patterns improve the representation capability of a depth-d node in the trie from 1 to 2 ^d unbinarized pattern(s), thereby reducing both the width of the automaton and the number of nodes. In this paper, the automaton with binarized patterns is named thinned automaton because of its reduced width. Here, we conduct a theoretical analysis of the improvement of the thinned automaton against the plain AC. For convenience, we assume that each character in a signature is uniformly distributed. To estimate the worst-case scenario, we assume that every signature has the largest length m observed in the set. While inserting a signature into a trie, the number of newly added nodes depends on the presence of its prefixes in the trie. Proposition 1 gives an expectation of finding prefixes for n signatures with c possible characters.

As shown in Proposition 3, the improvement with the thinned automaton is guaranteed under the assumption mentioned above. However, DNA sequences are biased. In this scenario, where the characters of each signature are not uniformly distributed, the improvement can be more pronounced because more duplicated segments lead to fewer trie nodes.

Even though the thinned automaton reduces the number of nodes, compressed representations may lead to collisions. Figure 4 shows an example of binarized results for five patterns and two sequencing reads, and Figure 5 further illustrates the corresponding thinned automaton. Two signatures CA and TG share the same binarized pattern 10 and result in a collision when reaching the yellow node in Figure 5. Substrings with identical binarized representations may also lead to false matches. For instance, ATGC in the second read, which is not a signature, has the same binarized representation 0101 as the signature ACAT.

Figure 4:

The binarized representations for five patterns and two sequencing reads. Two signatures GA and TC share the same binarized pattern (red). A substring in a sequencing read ATCG has the identical binarized form to the signature AGAT (blue).

Figure 5:

An example of the thinned automaton of AC with five signatures. Black solid links are trie links, and red dashed links are failure links. The yellow node represents two signatures GA and TC.

To maintain the correctness of signature profiling, each match to a binarized pattern needs to be verified with the original signatures. In other words, it is very time-consuming if there are serious collisions in certain nodes. A naïve comparison costs O ( ∑ p ∈ { p | p ′ = h , p ∈ P } | p | ) time to verify signatures with the same representation h.

Acceleration by rolling hash

Using hash functions is an intuitive idea to speed up comparisons between strings. As the lengths of signatures vary, arbitrary substrings of the read t ϵ  T is required to compute hash values during verification. However, on-the-fly computation of hash values takes an additional linear time O ( | t | ) for each checkup; pre-computing all possible substrings is also infeasible due to dispensable computations and extensive O ( | t | 2 ) additional memory.

To accelerate verification, we propose to apply rolling hash [39] that alleviates the time complexity for each checkup from linear to constant with a linear time pre-processing and additional linear memory consumption. Rolling hash is a family of hash functions where the input is hashed with a window that moves through the input. A new hash value can be rapidly calculated from the given old hash value in O(1) time. % time complexity. It also allows O(1) query time on the hash value of any substring in the input with content-based slicing. We implement the Rabin-Karp algorithm [40] as the rolling hash function. Formally, the hash value of a length-L input t [ 1 … L ] is defined as follows:

where t[i] is the i-th character of the input; a is a constant multiplier; q is a constant prime modulus. The hash value of a length-i prefix of t can be recursively calculated through the hash value of the length-(i − 1) prefix:

With bottom-up computation, hash values of all prefixes H ( t [ 1 … i ] ) can be preprocessed in both O(L) time and space complexity. Given the hash values of all prefixes, the hash value of a substring t [ i … j ] can be derived in O(1) as follows:

As a theoretical analysis, Proposition 4 gives a theoretical upper-bound of the collision probability. The larger the prime modulus q, the smaller the hash collision probabilities. Note that we employ the Rabin-Karp algorithm instead of cyclic polynomials for hashing because the former method has been demonstrated to be more efficient for general applications [41].

To apply rolling hash for acceleration, each node contains a hash table that maps a hash value onto the original signature. When transitioning to the node, the hash value of the matching substring in the read can be rapidly calculated and verified for its presence in the hash table.

As a result, the average time complexity of each checkup reduces to O(1). The overall time complexity of TahcoRoll is O ( ∑ p ∈ P | p | + ∑ t ∈ T | t | + ∑ p ∈ P c p ) , including the construction of the automaton and the matching process. The only memory overhead is hash tables with exactly | P | values, which is an amortized O ( | P | ) space.

Experimental datatsets

The performance of different algorithms is affected by four factors: signature lengths, number of signatures, read length, and number of reads. The flexibility of synthetic reads allows us to closely examine the effects of these factors. The randomly generated signatures are designed to test the worst scenario as their characters are uniformly distributed and may not share as many common prefixes as in the real sequencing applications. We also generate signatures from both genomic and transcriptomic sequences to analyze real reads from a diverse range of sequencing platforms. Synthetic datasets are available at https://figshare.com/s/6f02feaf89c4ff6ddc9e .

Implementation details

TahcoRoll builds a thinned automaton on a set of signatures represented by their canonical form (i.e., the lexicographical minimum of itself and its reverse complementary sequence).

Each node of the automaton holds an unordered map for average constant time complexity of searches and insertions after querying binarized patterns on the structure. The memory consumption of all automaton nodes is simultaneously pre-allocated for efficient memory operations. This is achieved by estimating the number of nodes through binary searches on sorted patterns. The profiling process scans each read twice, one for its forward sequence and the other for its reverse complementary sequence. The operations of the rolling hash are optimized by pre-computing the powers of the prime modulus.

In addition, the paralleled version applies the multi-threading capability of C++14. % for implementation. Unless otherwise mentioned, most of the experiments in this paper are conducted on a server with 504 GB of memory and two Intel Xeon E5-2680 v2 @ 2.80 GHz CPUs, where each CPU offers 10 cores.

Software adaption

We include all the k-mer counters mentioned in Section “Introduction” for experiments.

We also implement two baseline methods. The first one is a naïve implementation in C++, denoted by “Naïve”. It uses a hash table to store k-mers and scans through the reads multiple times with different window sizes. Theoretically, Naïve is light in memory, but requires an extensive running time. The second baseline is the conventional AC. We test two publicly available implementations written in Python (PlainAC_Py; pyahocorasick 1.1.3) and C++ (PlainAC_C++; cjgdev/aho_corasick).

Automaton construction

The memory of AC is sensitive to the composition of signature patterns such as k-mer lengths, the number of k-mers, and common prefixes shared by different k-mers. Figure 6 compares the computational resources used for automaton construction in TahcoRoll against PlainAC_Py and PlainAC_C++ over 16 sets of signatures. The implementation of PlainAC_C++ uses several additional data structures on each node to facilitate the traversal on an automaton, causing a huge memory overhead. As a result, PlainAC_C++ is fast in automaton construction but requires twice and five times more memory than PlainAC_Py and TahcoRoll, respectively. For the large batch of 24 million k-mers, PlainAC_C++ maxes out the memory capacity (>396 GB) of our server, and thus the recorded run-time is truncated. Our thinned automaton consistently requires less time than PlainAC_Py in construction. As we increase the number of k-mers, the construction time rises. The memory of the thinned automaton is significantly reducing to nearly half of the memory required in PlainAC_Py.

Figure 6:

Run-time and memory for constructing the automaton given 16 sets of k-mer patterns. A lower value represents a more efficient approach. PlainAC$_C++ maxes out the memory capacity while constructing the “large” batch of 24 million k-mers (large_240), and thus its recorded time and memory are truncated. TahcoRoll consistently requires less time and memory than PlainAC_Py.

Pilot study of 13 approaches

We perform a preliminary assessment of the memory footprint and run-time on 11 existing counters, together with Plain_AC and TahcoRoll.

We separate the analyses into two panels as demonstrated in Figure 7. The top panel focuses on different numbers of reads, and the bottom panel focuses on different numbers of k-mers. Methods in the bottom-left corner of each plot indicate being both time and memory efficient. As we predicted, Naïve uses very little memory, but takes a long time to complete. PlainAC is fast but requires a large amount of memory when increasing the number of k-mers. Consistent with the analysis in Figure 6, PlainAC_C++ uses twice as much memory as PlainAC_Py.

Figure 7:

Run-time (x-axis) and memory (y-axis) of counting small batches k-mers on synthetic reads of 75 bp. Each point represents a pair of measurements (run-time and memory). Points on the lower left corner of each plot indicate more efficient approaches. The top panel examines different read sets with 1.2 million k-mers; the bottom panel examines different k-mer sets with 34,497,448 reads.

TahcoRoll is the most efficient approach in five out of these six analyses. KMC3 and Squeakr use less memory when there are 24 million k-mers, but requires more time than TahcoRoll. When we fix the number of k-mers (top panel), the memory and run-time for KMC3 and Squeakr increase with the number of reads, but the memory stays constant for both TahcoRoll and Jellyfish. Jellyfish is memory efficient when counting a given list of k-mers with the same size; however, repeating this process for different k’s makes it more time-consuming than TahcoRoll.

Extensive study on synthetic datasets with both single and multiple threads

We use 1.2 million k-mers ranging from 15–151 bp (wide) to evaluate the scalability on different read lengths and number of reads. We highlight the total run-time and memory consumption of each approach in Table 1. The run-time is further broken down into the automaton construction phase (Prep) and the read querying phase (Query) for TahcoRoll and PlainAC_Py. The two phases of MSBWT and KMC3 include indexing the reads (Prep) and querying the k-mers (Query). Read processing is performed in the querying phase of TahcoRoll and PlainAC_Py, but in the preparation phase of MSBWT and KMC3. Therefore, the run-time of querying is not on the same scale across different approaches. For Jellyfish, we use its function to count the list of k-mers directly, so the run-time cannot be split in details. Its memory usage depends on the size of the list of k-mers, and can be as efficient as TahcoRoll. However, its run-time does not scale well with datasets containing more or longer reads. TahcoRoll consistently outperforms others on different read sets in run-time and memory.

Next, we use 86,976,737 reads of 180 bp to evaluate the scalability on different batches of k-mers, which are designed to test the worst scenario. Table 2 shows that when the k-mers are short (small), PlainAC_Py uses the least amount of time. When k-mers get longer, TahcoRoll is the most efficient approach. This observation is due to less collision in the signature sets. Under a severe condition where there is a large number (12 and 24 million) of k-mers with uniformly distributed characters, TahcoRoll requires more memory than MSBWT in three out of six cases. It is worth mentioning that MSBWT and KMC3 write a huge amount of intermediate files to disk (at least 16 GB for MSBWT and 43 GB for KMC3 for this dataset) to alleviate the memory bottleneck. In contrast, TahcoRoll is an in-memory approach that does not generate any intermediate data.

MSBWT, KMC3, and Jellyfish allow indexing reads in parallel, so we evaluate the parallel settings on the wide batches of k-mers. Table 3 shows the run-time and memory usage of analyzing 86,976,737 synthetic reads of 180 bp across four sets of k-mers. Both Jellyfish and TahcoRoll scale well with the number of threads, but the improvement of MSBWT and KMC3 is marginal. This is mainly due to the limitation of I/O as these two approaches constantly read and write files to disk. The run-time of TahcoRoll remains faster than others across different experiments and threads. The four-thread TahcoRoll demonstrates to be faster than others with 16 threads.

Data from different sequencing platforms

Here, we examine the practical usage by analyzing signatures from real DNA sequences with reads from different sequencing platforms. Table 4 summarizes the nature and analysis of each dataset. TahcoRoll is run with a single-thread and eight-threads; others are run with eight-threads. For the measurement that is less efficient than TahcoRoll, we compute the fold-change to those reported by the eight-thread TahcoRoll. MSBWT is unable to finish indexing for the PacBio data within two days. KMC3 cannot index long reads from Nanopore as the data exceeds the buffer size set by the program; it also consumed all memory available on the machine (32G) for PacBio data. Overall, the run-time of single-thread TahcoRoll is as efficient as Jellyfish with eight-threads and significantly outperforms KMC3 in short reads and MSBWT in long reads. In parallel settings, TahcoRoll runs at least four times faster than MSBWT and Jellyfish and demonstrates a drastic improvement over KMC3.