Sequence Similarity Estimation by Random Subsequence Sketching

Recently, a sketching method named LexicHash [9] proposes to compare sketches based on the length of their common prefixes, rather than relying on fully matched $k$-mers. This has the effect of sketching with $k$-mers for all lengths $k$ up to a predefined maximum value. However, LexicHash still suffers from the common issue of $k$-mer-based methods, namely, a small number of edits can destroy all long $k$-mer matches between two similar sequences. Furthermore, LexicHash is designed for the task of overlap detection, rather than estimating the similarity between two sequences. In particular, the authors define the LexicHash similarity score between two sequences as the length of the longest matching prefix among their sketches. Hence, a score $k$ only indicates that the two sequences share a common $k$-mer, which may be effective for detecting overlapping reads, but appears to be insufficient for edit similarity estimation (see Figure 4). In their follow-up work LexicMap [26], the authors design a sequence alignment tool that uses the identified common $k$-mers as seeds, effectively employing LexicHash as a seeding method rather than a sketching method. See Section 4 for further discussion of the differences between sketching and seeding.

To the best of our knowledge, the only existing subsequence-based sketching method is Tensor Slide Sketch (TSS) [12]. Instead of picking $k$-mers from the input sequence, TSS aims at producing a sketch by counting the occurrences of all subsequences. Since there is an exponential number of them, TSS has to group subsequences in a smart way to facilitate counting. However, to make it efficient, TSS is restricted to count all short subsequences, which limits its capacity in distinguishing similar and dissimilar sequences.

A closely related line of research focuses on embedding edit distance into metric spaces. Most notably, the CGK embedding [4] maps sequences into Hamming space with bounded distortion and, when combined with locality-sensitive hashing techniques, has been applied to various bioinformatics tasks, including edit similarity joins [29] and overlap detection [27]. However, CGK maps a sequence of length $n$ to one of length $3n$, making it unsuitable as an efficient sketching method, where the goal is to reduce representation size. To take advantage of the structure of the underlying data, machine learning methods have been developed that learn embeddings outperforming CGK in downstream applications [5], despite lacking its theoretical guarantees. These approaches, however, require training data with known edit distances, limiting their applicability to tasks involving long sequences, such as genome phylogeny.

A sequence $x$ of length $kt$ over an alphabet $\mathrm{\Sigma }$ can be viewed as a sequence of $k$ “tokens” each of which is a string of length $t$. We say $x$ is a tokenized subsequence of a length-$n$ sequence $s$ if there is a list of indices such that the length-$t$ substring of $s$ starting at $i_j$ matches the $j$-th token of $x$. Note that when $t=1$, a tokenized subsequence is a regular subsequence; it is not necessarily the case when $t>1$, as the tokens are allowed to overlap, see Figure 1 for an example.

To construct SubseqSketch for input sequences, we first generate a list $L$ of random sequences of length $kt$ each, where $k$ and $t$ are predefined parameters. We call $L$ the list of testing sequences. Two SubseqSketches are comparable only if they were generated with the same list $L$; in this sense, $L$ serves as shared randomness in the sketching process, analogous to the shared random ordering of $k$-mers in MH sketches. Given an input sequence $s$ and a testing sequence $\mathrm{\ell }$ in $L$, SubseqSketch determines the maximum number of prefix tokens in $\mathrm{\ell }$ that form a tokenized subsequence of $s$. The resulting vector consists of $|L|$ integers, one for each testing sequence. This vector is the sketch of $s$, denoted as $\text{SubseqSketch}(s)$. See Figure 2 for an illustration.

A straightforward linear scan computes the $|L|$ sketch entries in $O(|L||s|)$ time. This worst-case time complexity can be improved by preprocessing the input sequence $s$ to build an index that facilitates rapid lookup for the occurrence of the next token of a testing sequence. For example, for token size $t=1$, we can build an automaton on $s$ in $O(|s||\mathrm{\Sigma }|)$ time and space. In the automaton, each character $s_i$ stores $|\mathrm{\Sigma }|$ pointers. The pointer corresponding to $c\in \mathrm{\Sigma }$ points to the next appearance of $c$ after $s_i$ (or null if no $c$ exists after $s_i$). Then for a testing sequence $\mathrm{\ell }$, we can simply follow the pointers according to the characters of $\mathrm{\ell }$, until either a null pointer is encountered or $\mathrm{\ell }$ is exhausted. This takes $O(|\mathrm{\ell }|)$ time for each testing sequence so the total sketching time is $O(|s||\mathrm{\Sigma }|+|\mathrm{\ell }||L|)$.

For larger token size, a similar idea can be applied: we can preprocess $s$ to build a lookup table of size ${|\mathrm{\Sigma }|}^t$ where each entry records the occurring positions of that token on $s$, either in a sorted array or some other data structures that supports quick search. Each testing sequence can then be processed by following this lookup table until all tokens are used or the end of a position array is reached. This allows each integer in the sketch to be computed in $O(k\log |s|)$ time, instead of a $O(|s|)$ linear search. We provide this preprocessing approach as an option in our implementation. However, through experiments we found that the linear search std::string::find provided in the standard C++ library is almost always faster. The overhead of preprocessing may only be justified for a large number of very long testing sequences with a small token size, which is not a recommended setting for our sketching algorithm (see Section 2.4).

The SubseqSketch of a sequence $s$ provides a highly informative representation of $s$. To build intuition, consider two sequences $s$ and $t$. If both sketches show large numbers at the same index, then $s$ and $t$ must share a long tokenized subsequence and are hence likely similar in terms of the edit distance. Conversely, if one sketch has a large value while the other has a small value at the same index, it suggests that the sequences are likely dissimilar.

As with other sketching methods, a similarity measure over the sketches is required to translate the above intuition into a quantitative score that accurately reflects the true similarity between input sequences. Methods that compare sketches for equality at corresponding indices, such as MH and OMH, naturally employ Hamming similarity, which counts the number of matching entries between sketches. SubseqSketch, on the other hand, generates integer-valued vectors, enabling the use of a wide range of well-established distance/similarity metrics. We empirically evaluate a list of metrics using the data from Section 3.1. SubseqSketches are first computed, after which similarity scores are calculated using various metrics. The Pearson correlations between these scores and the ground truth edit similarities are reported in Table 1.

As shown in the table, cosine similarity is among the most effective metrics for producing estimates that are strongly correlated with the true edit similarity. According to its definition,

where $\cdot$ denotes the vector dot product, pairwise cosine similarities between two sketching matrices can be computed using a single matrix multiplication (assuming the rows are normalized), which is highly optimized in modern hardware and numerical libraries. We therefore adopt cosine similarity between SubseqSketches in our implementation for its effectiveness and computational efficiency.

SubseqSketch has three parameters: the token size $t$, the number of tokens $k$ in each testing sequence, and the size $|L|$ of the testing list. The parameter $|L|$ controls the size of the sketches. In particular, a SubseqSketch takes $|L|\log k$ bits space to store. As with other sketching methods, increasing the sketch size improves estimation accuracy but comes at the cost of greater time and storage requirements. In the experimental sections, we compare the sketching methods at the same sketch size.

The parameters $t$ and $k$ are related. In the resulting sketches, each entry is an integer between $0$ and $k$. If $t$ is too large (for example, close to the input length $n$), most entries would be $0$; on the other hand, if both $t$ and $k$ are small, most entries would reach the maximum possible value at $k$, regardless of the input sequence $s$. Neither case is desirable as the sketches cannot provide a strong distinction between similar and dissimilar input sequences. Note that we can always choose a large $k$ to ensure that few, if any, sketch entries reach the maximum value. However, this increases the sketch file size, as each entry requires $\log k$ bits – an inefficient use of space if most entries are significantly smaller than $k$.

We now try to derive an optimal choice of $k$ for $t=1$. In a recent paper [8], the authors motivated their sequence sampling method with an interesting puzzle (paraphrased): is the number of DNA 5-mers containing the substring ACGT the same as that for the substring AAAA? Astute readers will immediately answer “no” because it is impossible for a 5-mer to both start and end with ACGT – taking the union of the two disjoint groups gives the correct number – which is not the case for AAAA whose symmetry would cause the same strategy to double-count the 5-mer AAAAA.

As a curious extension, the same question can be asked, replacing substring with subsequence, namely, we do not require the containment to be consecutive. This seemingly more complicated version turns out to have a counterintuitively nicer answer: the number of $n$-mers containing a given $k$-mer as a subsequence is a function of $n$ and $k$, independent of the choice of the $k$-mer. Consider a length-$k$ sequence $x$, we count the number of length-$n$ sequences $s$ whose subsequence is $x$. To avoid over-counting, we only count $s$ if $(i_1,\mathrm{\dots },i_k)$ is the first occurrence of $x$ in $s$. It means the characters in $s$ before $i_1$ cannot be $x_1$, leaving them $|\mathrm{\Sigma }|-1$ choices each. The same holds for regions in between $i_j$ and $i_{j+1}$, and finally all characters after $i_k$ are free to be anything in $\mathrm{\Sigma }$. This leads to ${(|\mathrm{\Sigma }|-1)}^{i_k-k}{|\mathrm{\Sigma }|}^{n-i_k}$ choices. Note that the expression only depends on $i_k$ (i.e., any combination of $i_1,\mathrm{\dots }{i}_{k-1}$ yields the same number), so we can group the terms and sum over choices of $i_k$ to get the answer

We emphasize that the above calculation is independent of the chosen sequence $x$. An example is shown in Figure 3.

Given the length $n$ of the input sequences, we can use this formula to find a $k$ such that at most a small threshold fraction (e.g., $0.01$) of the sketch entries reach the maximum value $k$. For the example in Figure 3 where $n=100$, $k=36$ would suffice.

For larger $t$, the derivation is not as neat. We can view a regular sequence $s$ over the alphabet $\mathrm{\Sigma }$ as a tokenized sequence over the alphabet ${\mathrm{\Sigma }}^t$ and apply the above formula. But unlike adjacent characters in the original sequence, consecutive tokens with an overlap of length $t-1$ are not independent, causing the formula to significantly overestimate. Since using a small $k$ makes the sketching faster to compute and smaller to store, with an exception in Table 3, we fix $k=15$ in the following experiments (namely, each entry in the sketch fits in 4 bits) and aim to choose $t$ to ensure the sketching entries are neither too small nor maxed out. Table 2 provides empirical recommendations for $t$ across common input sizes $n$.

Using randomly generated testing sequences is the best one can do in a data-oblivious setting, while better performance can usually be achieved if we can afford to adjust the sketches according to the input data. One idea to introduce data dependency is to sample subsequences from the input to form the testing list. This is particularly suitable when the input comprises a small number of sequences – for example, when estimating phylogenetic distances among a group of closely related genomes, as shown in Section 3.3. On the other hand, if the sketches are used to build an index of a large database of sequences to handle queries, it may not be practical to re-sketch the entire database with a new testing list for each query. In this situation, we simply use the data-oblivious version with a fixed list of randomly generated testing sequences and demonstrate in Sections 3.1 and 3.2 that it already achieves good performance.

To directly compare the sketch similarity against the desired but much more expensive to compute edit similarity, we generate $100,000$ random DNA sequences of length $1,000$. Each sequence is randomly mutated (an insertion, deletion, or substitution) for a random number of rounds up to $1,000$ to produce a pairing sequence. For each pair, we compute their exact edit similarity, as well as sketch similarities for SubseqSketch, MinHash (MH), Order Min Hash (OMH), Tensor Slide Sketch (TSS) and LexicHash (LH). For each sketching method the Pearson correlation between the exact edit similarity and the sketch similarity over the 100,000 pairs of sequences is reported. MH, OMH, and TSS use the implementation of [12]. LH uses the implementation of [9].

Figure 4 shows the scatter plots of all the pairs under different sketching methods. Observe that SubseqSketch achieves the best Pearson correlation. Both MH and OMH are good estimators for sequences with high edit similarities but struggle to distinguish dissimilar sequences with edit similarity between $0.5$ and $0.8$. The TSS and LH similarities show a visually more linear relationship with the edit similarity and consequently exhibit higher Pearson correlations than MH and OMH. But they both suffer from extremely large variance, especially for dissimilar sequences, which makes it difficult to interpret their estimation in practical applications. SubseqSketch strikes a balance between the ability to estimate the full range of edit similarity and the estimation variance.

As with other sketching methods, the variance of SubseqSketch can be reduced by using a larger sketch. For all the experiments, we measure the size of a sketch as the number of entries in it (sometimes called its dimension), and all methods are configured to produce the same number of entries (except for TSS, which we follow the suggestion in [12] even though it produces a larger sketch). However, in real applications, the actual space needed to store the sketches is a more relevant measure. Recall that each entry of SubseqSketch can be stored in 4 bits (ref. Section 2.4) which is four times smaller than an entry of MH (16 bits for $k=8$), six times smaller than OMH (24 bits for $k=6$ and $\mathrm{\ell }=2$), and eight times smaller than TSS and LH (32-bit float/int). Thus, given a fixed amount of disk space, SubseqSketch can utilize more testing sequences than the number of $k$-mers MH or OMH can select, thereby achieving a similar or better variance. In the experiments, we do not exploit this practical advantage, opting instead to use the same number of sketch entries across all methods.

The task of nearest neighbor search involves finding the top-$T$ most similar sequences for a query among a large database. It has numerous applications in bioinformatics, including metagenomic classification and gene function prediction. Because computing the exact edit distance between the query and every sequence in the database is computationally prohibitive, a common approach is to embed database sequences into a well-studied metric space where efficient nearest neighbor indexing is readily available (for example, the hierarchical navigable small world index [18]). The query is embedded in the same space, and the nearest neighbors identified by the index serve as approximations to the true nearest neighbors under edit distance. Given the high correlation between SubseqSketch and edit distance, we demonstrate in this section that it can serve as an embedding method and yields improved nearest neighbor search performance.

In this experiment, we restrict our comparison to CNN-ED [5], a deep convolutional neural network model for sequence embedding that has been shown to outperform non-machine learning approaches such as CGK embedding [4]. For evaluation purpose, the CNN-ED pipeline splits an input dataset into three disjoint sets: a training set with $1,000$ sequences, a query set with $1,000$ sequences, and a base set containing the remaining sequences. It computes pairwise edit distances among the training set to generate training labels. After the model is trained, it generates embeddings for all sequences in the query and base set. All-vs-all embedding distances as well as the ground-truth edit distances between the query set and the base set are then computed so that the top-$T$ neighbors identified by the embedding can be compared with the top-$T$ neighbors according to the edit distance. Following this pipeline, we apply SubseqSketch to all sequences in the query and base set, and report top-$T$ neighbors using the all-vs-all sketch distances. It is worth noting that, unlike in practical applications, indexing is omitted in this experiment to enable a direct comparison of the embedding/sketching quality, independent of variations in index accuracy.

We show results on two widely used datasets GEN50kS and GEN20kL from [29] which are also benchmarked in the CNN-ED paper. The GEN50kS dataset contains $50,000$ sequences with an average length $5,000$. The GEN20kL dataset contains $20,000$ sequences with an average length $20,000$.

In Figure 5 and Figure 6, we plot the commonly used recall-item curves for both SubseqSketch and CNN-ED. For a figure labeled top-$T$, the $T$ nearest neighbors of a query in the base set according to the edit distances are considered true neighbors. The horizontal axis represents the number of neighbors (items) each method is allowed to report (according to their respective sketch/embedding distances) and the vertical axis marks the fraction of true neighbors being reported (recall). The CNN-ED pipeline presents full-range results – from reporting a single item to reporting all items – which, while not practical for typical use cases (where only the top-$T$ neighbors are retrieved), allows for plotting complete performance curves.

The CNN-ED results are obtained by the implementation of [5], which we trained for 50 epochs following the reported hyperparameters in the original paper. For a fair comparison, SubseqSketch is configured to produce vectors of the same length as the reported embedding dimensions of CNN-ED. Observe that SubseqSketch consistently outperforms CNN-ED by a large margin. This is a surprising result. It is commonly believed (which is often, though not always, justified) that machine learning models can outperform traditional algorithmic methods because the models can learn data-dependent features that the data-oblivious algorithms cannot take advantage of. In [5], the CGK embedding [4] was shown to produce a worse result than CNN-ED on this task, even though it is an edit distance embedding with theoretical guarantees. Our result here demonstrates that there is a gap between theoretical bounds and practical performance which warrants further investigation. In particular, we conjecture that SubseqSketch can also provide some guarantees on the distortion as a randomized embedding function for the edit distance, though a theoretical proof seems difficult.

Phylogeny reconstruction is another task commonly used to evaluate the performance of sketching methods. Given a set of biologically related genomes, the goal is to build a phylogeny on them based on pairwise similarities/distances estimated by the sketches. The result can then be compared with a ground truth tree constructed from some biological model or multiple sequence alignment. We test on two datasets for this task: one is a simulation of a simple mutation model similar to that used in [20]; the other is a set of 29 assembled E. coli genome sequences collected in [28].

The simulation aims to model both point mutations and mobile genomic elements, commonly found in bacterial genome rearrangements, known as insertion sequences (IS). The simulated sequences form a perfect binary tree. The root of the tree is a random sequence of length $10,000$; it is considered as the $0$-th generation genome. To obtain the $i$-th generation, each sequence in the $(i-1)$-th generation produces two children genomes by independent and random point mutations with mutation rate $0.01\%$. Then a random IS of length $500$ is inserted at a random position for each newly generated $i$-th generation genomes. Note that the IS is shared among all sequences in the same generation, but the inserting positions can be different. See Figure 7 for an illustration. Although simple and somewhat unrealistic, this model produces a solid ground truth phylogeny and allows us to investigate the effectiveness of different sketching methods to recover the mixed history of point mutations and large insertion events.

Each method computes a pairwise distance matrix for all sequences in each generation. The matrices are used to build the phylogenies with the neighbor-joining algorithm implemented in the biotite package [13]. The normalized Robinson-Foulds (nRF) distances between the constructed trees and the ground truth tree are then calculated with the ETE toolkit [10]. The nRF distance measures the dissimilarity of branching patterns between two trees and ignores branch lengths. A value of $0$ means the two phylogenies have the identical tree topology, whereas a value of $1$ indicates the two trees are maximally dissimilar.

Figure 8 shows the nRF distances achieved by each method on progressively larger inputs from the simulated dataset. The horizontal label $i$ means all the $2^i$ sequences from the $i$-th generation are used as input sequences. Not surprisingly, pairwise edit distance (ED) most accurately captures the mutation history, at the cost of significantly longer computation time (see Figure 9). Among the sketching methods, SubseqSketch constructs the best phylogeny for generation $6$ and larger inputs. Furthermore, the nRF distances obtained by SubseqSketch exhibits a strong correlation with those achieved by the exact edit distances, indicating it can be used as a faithful approximation of the expensive edit calculation. In contrast, although MH and OMH produce trees with smaller nRF distances for the smaller input sets, they both show some inverse relation with the nRF using edit distances (e.g., from generation $3$ to $4$, the nRF distances of trees constructed by edit distance increased, but the nRF distances for MH decreased; similarly from generation $4$ to $5$ for OMH). LH is omitted from this experiment because its implementation choice for boundary handling tends to assign the maximum similarity score to pairs sharing a short matching suffix (see the line at normalized similarity score $1$ in Figure 4). While this may be appropriate for the overlap detection task that LH is designed for, it hinders accurate phylogeny reconstruction on our datasets.

We also plot the running time of each sketching method in Figure 9 to demonstrate the efficiency of SubseqSketch. As expected, all the sketching methods are much faster than computing the all-vs-all exact edit distances. Among them, SubseqSketch is consistently the fastest, regardless of the number of input sequences. More specifically, SubseqSketch achieves a $6\times$ speedup compared to the second fastest method (MH).

For testing on a real dataset, 29 assembled E. coli genome sequences are obtained from the AFproject [30], a benchmark project for alignment-free sequence analysis tools. It provides a web interface where the phylogenies can be computed from the uploaded distance matrices. The nRF distances are then reported by comparing the resulting trees against a ground truth tree built from multiple sequence alignment. It also provides the normalized Quartet Distance (nQD) as an additional measure for topological disagreement.

Results on this E. coli dataset are summarized in Table 3. On the AFproject website, nearly 100 tools (include different configurations for the same tool) are ranked based on the nRF distance. SubseqSketch is ranked 7th and there are 12 tools that achieve smaller nRF distances due to ties. It is worth pointing out that the higher ranked ones are tools designed specifically for the task of phylogeny reconstruction, which are often based on some sketching method but also apply biological and algorithmic heuristics to adjust the sketch distance matrix. Since SubseqSketch is a sketching method rather than a complete tool for phylogeny, here we aim to evaluate the sketch quality without those adjustments. With the unadjusted distance matrices, SubseqSketch constructs the best phylogeny (closest to the ground truth) among MH, OMH, and TSS.

In this task, since there are only $29$ genomes, we can afford to sample the testing sequences from the input to further improve the quality of SubseqSketch. Because the inputs are all closely related, this sampling strategy also enables us to use a much larger token size $t=40$ to achieve an even better result than the recommended $t=15$. From Table 3, it is evident that setting $t=40$ significantly improves accuracy.