Algorithms for efficiently collapsing reads with Unique Molecular Identifiers

Connected components

One method for identifying groups of unique UMIs is based on connected components (known as clusters in UMI-tools (Smith, Heger & Sudbery, 2017) and Calib (Orabi et al., 2018)). A connected component in a graph G = (V, E) is a set of vertices S⊆V where there is a valid path between each pair of vertices u, v ∈ S, and there is no edge (u, v) ∈ E where u ∈ S, but v⁄ ∈ S. Each connected component is a group of UMIs, and the UMI that has the highest frequency is chosen to represent the group. This algorithm allows UMIs that are similar under the Hamming distance metric to be grouped together, but it can lead to an underestimation of the total number UMIs. There may exist an UMI w for two UMIs u and v, where d(u, w) ≤ k and d(w, v) ≤ k, causing u and v to be grouped together even though d(u, v) > k (essentially “bridging” the Hamming distance “gap” between the two different UMIs). We show the entire connected components algorithm in Algorithm 1.

The time complexity of finding all connected components is O(N) with respect to each query to remove_near, since only one remove_near query is allowed per vertex (UMI) in the graph, and each vertex is only visited once.

 
__________________________________________________________________________________________ 
Algorithm 2 The adjacency algorithm on a set of UMIs V . 
__________________________________________________________________________________________ 
  procedure get_groups_adjacency(V , k) 
     Initialize a data structure that implements remove_near and contains with V 
  and k. 
     adj ←{}                             ⊳ Resulting set of groups of adjacent UMIs. 
     Sort V by decreasing UMI frequency values. 
     for u ∈ V do 
         if contains(u) then 
             adj ← adj ∪ remove_near(u, k, ∞) 
         end if 
     end for 
     return adj 
  end procedure 
______________________________________________________________________________________________    

Adjacency

To avoid the underestimation issue of the connected components algorithm, we can sort the UMIs by their frequencies, and examine UMIs from higher to lower frequency. We assume that UMIs that appear at a higher frequency are more likely to be correct, and allow adjacent UMIs in G that have lower frequencies to be grouped with a higher frequency UMIs. Two UMIs u and v are adjacent if there exists an edge (u, v) ∈ E. The lower frequency UMIs are assumed to be due to substitution errors. The full algorithm is presented in Algorithm 2.

Since we must first sort the N UMIs in O(NlogN) time and then query remove_near at most N times, the overall time complexity is O(NlogN + NQ), where Q is the time complexity of remove_near. If Q ≥ logN, then the time complexity with respect to each query is O(N).

Directional

The directional algorithm allows UMIs that are not adjacent to be grouped together, and it also minimizes underestimation from examining connected components. Like the adjacency method, it involves iterating through the sorted UMIs from highest to lowest frequency. Instead of grouping together all adjacent vertices like in the adjacency method, we add a vertex v that is adjacent to vertex u to the group iff f(v) ≤ ϵ[f(u) + 1], where 0 ≤ ϵ ≤ 1 and ϵ is a parameter that represents the threshold percentage at which adjacent UMIs are grouped, which is usually ϵ = 0.5 (this is a rearrangement of UMI-tools’ inequality 2f(v) − 1 ≤ f(u) (Smith, Heger & Sudbery, 2017)). Then, the process is recursively repeated starting at v by examining the adjacent vertices of v. This allows UMIs that are not directly adjacent to a high frequency UMI to be visited and grouped with the high frequency UMI if they have a low frequency. Algrorithm 3 represents the full algorithm.

 
__________________________________________________________________________________________ 
Algorithm 3 The directional algorithm on a set of UMIs V . 
__________________________________________________________________________________________ 
  procedure get_groups_directional(V , k, ϵ) 
     Initialize a data structure that implements remove_near and contains with V 
  and k. 
     dir ←{}                                    ⊳ Resulting set of groups of UMIs. 
     Sort V by decreasing UMI frequency values. 
     for u ∈ V do 
         if contains(u) then 
              dir ← dir ∪ dfs_dir(u, k, ϵ) 
         end if 
     end for 
     return dir 
  end procedure 
  procedure dfs_dir(u, k, ϵ) 
     g ←{}                         ⊳ Group of UMIs found by the directional algorithm. 
     for v ∈ remove_near(u, k, ϵ[f(u) − 1]) do 
         if v ⁄= u then 
              Add all elements of dfs_dir(v, k, ϵ) to g. 
         end if 
     end for 
     return g 
  end procedure 
_____________________________________________________________________________________________    

The overall time complexity of directional is O(N), which is similar to that of adjacency, since they both require an initial sorting of the N UMIs, and each UMI requires a remove_near query.

Constant-time Hamming distance calculations

Let ⊕ denote the bitwise XOR operation, and let the pop_count function return the number of set (1) bits in a binary string. For two binary strings a and b where a, b ∈ {0, 1}^M, the Hamming distance (total number of mismatches) between them can be easily calculated by pop_count(a⊕b). This operation can be done in O(1) time if both a and b are shorter than the word length in a computer, which is typically 64-bits. We show that this idea can be extended to all four nucleotides ({A, T, C, G}).

The main property we want for an encoding of the four nucleotides is that the encodings are pairwise equidistant. One such set of encodings with the shortest possible length for each encoding is ${\displaystyle \begin{array}{c}\hfill e_{\text{A}}=110e_{\text{T}}=011\\ \hfill e_{\text{C}}=101e_{\text{G}}=000\end{array}}$ These encodings can be viewed as the vertices in a regular 3-simplex (tetrahedron) in Hamming space, and they are discussed in Liu (2019). The Hamming distance between each pair of different encodings is exactly two, which means that for two UMI sequences u and v, we can translate them into bit strings a and b using the encoding method above, and the overall Hamming distance between the two bit strings is exactly 2d(a, b). This can be used to easily infer d(a, b). For example, let “AAT” and “AAA” map to 110110011 (concatenate e_Ae_Ae_T) and 110110110 (concatenate e_Ae_Ae_A), respectively. Then, pop_count(110110011⊕110110110) = 2d(AAT, AAA) = 2. Thus, the edit distance d(AAT, AAA) = 1.

Encoding each UMI sequence saves time and memory, since multiple encoded nucleotides can be stored in a computer word, and the Hamming distance between two UMIs can be calculated in constant time per word by using the XOR and pop_count operations. $\lfloor \frac{64}{3}\rfloor =21$ nucleotides can be packed into one computer word of length 64, and an array of words can be used to handle longer UMI lengths. Note that we separately encode the positions of undetermined N nucleotides in sequences that contain it. We treat this undetermined nucleotide as different from the other four nucleotides.

For all other data structures and algorithms that we discuss, we assume that the UMI lengths are short, and Hamming distances are calculated with this method, which takes O(1) time per calculation. This also speeds up hash operations and equality comparisons between UMI sequences to O(1) time, since they can process the entire encoded UMI sequence at once.

Combinations

Instead of a brute-force search through all N UMIs during remove_near queries, it may be faster to directly generate UMIs sequences that are within a certain edit distance and checking whether each of them exists using a hash table. The run time complexity of this “combinations” algorithm, for each remove_near query with an alphabet Σ = {A, T, C, G} of length |Σ|, is $O\left(|\Sigma {|}^k\left(\genfrac{}{}{0.0pt}{}{M}{k}\right)\right)=O\left(|\Sigma {|}^k{M}^k\right)$. This is because any k nucleotides in the UMI can be edited, and each nucleotide can be edited to any nucleotide in Σ. Each combination can be checked in O(1) time on average with the hash table that contains all N UMIs. This method is efficient if k is small, but it scales exponentially with k, and it does not make use of our constant-time Hamming distance calculations. This method is used in umis (2019) for UMI deduplication.

Trie

The main problem with the combinations algorithm is that it may examine multiple UMI combinations that do not actually exist in the input UMIs sequences. To address this problem, a string trie (De La Briandais, 1959) can be built to store each of the N UMIs. Each node in a trie consists of up to |Σ| children—one for each character in Σ. Also, each unique path from the root of the trie to a leaf node at depth M represents a unique UMI of length M, with each node on the path representing a nucleotide in the UMI sequence. A trie allows prefixes that are shared among multiple sequences to be collapsed into one prefix path, which saves space. However, the main benefit of a trie is that while generating UMI combinations, the current prefix of the UMI combination can be checked against nodes that exist in the trie. This allows prefixes of UMI combination that cannot possibly construct an existing UMI to be pruned, which speeds up queries in practice. An example of a trie is shown in Fig. 1.

To build the trie, around MN operations are necessary to insert every UMI sequence into the trie, for a total of O(MN) time. Each remove_near query will take at least O(MR) time, if R total UMI sequences are returned by the query, and the exact time taken depends on how many UMIs share prefixes and k, the number of edits allowed. The factor of M is due to the trie’s height being exactly M nodes.

Some subtrees in the trie may become completely removed after some remove_near queries, when all UMIs in that subtree are removed. These subtrees can be completely skipped during future remove_near queries to save time while processing future queries, by keeping track of a flag that indicates whether the subtree of each node still exists. This flag can be recursively propagated up towards the root of the trie during remove_near queries (i.e., each node’s subtree is completely removed iff all children subtrees of that node are completely removed). Then, an entire subtree can be pruned if the “completely removed” flag is set, during remove_near queries.

Similarly, we can keep track of the minimum frequency of any UMI within each subtree of the trie. After each query and removal, the minimum frequency of each node’s subtree can be recursively updated for non-leaf nodes, to exclude the UMI frequency of any removed UMIs (i.e., each node’s subtree’s minimum frequency is the minimum across the minimum frequencies of that node’s children subtrees that are not completely removed). This allows an entire subtree of the trie to be discounted if its minimum frequency is higher than F, since that implies that no UMI within the subtree has a lower UMI frequency than the maximum threshold F.

n-grams

An alternative method for speeding up the naive O(N) time remove_near query is by filtering UMIs based on partial exact matches of UMI sequences. The benefit of exact matches is that they can be efficiently indexed in a hash table, and therefore search operations would take O(1) time on average. The general idea of the n-gram algorithm is to first decompose each UMI sequence into contiguous, non-overlapping sequences (n-grams) of length n during the initialization of the n-gram hash table. Each unique n-gram, which is represented by both its location in an UMI sequence and its actual nucleotide sequence, maps to a bin (list) of UMI sequences that contains the n-gram at its specific location in the UMI sequence. Then, to search for all UMIs within k Hamming distance away from the queried UMI sequence, we can first decompose the query sequence into n-grams with the same method used during initialization, and only search the bins that corresponded to the n-grams created from the queried UMI. This is used by later versions of UMI-tools (Smith, Heger & Sudbery, 2017) to speed up the construction of its networks for UMI deduplication.

The UMI sequences must be split into k + 1 contiguous and non-overlapping n-grams, where the n-grams are all approximately equal in length. For two UMI sequences that are within k edits apart, and are split into k + 1 n-grams, there must always be at least one n-gram that exactly matches in both sequences. This must be true, since in the “worst-case scenario” when there are exactly k edits and k of the n-grams have an edit within them somewhere, there is still exactly one n-gram that perfectly matches since it must have no edits due to the edit threshold of k. We also want to maximize the length of each individual n-gram, in order to maximize the number of distinct n-grams and prune more of the search space. Therefore, the length of each n-gram must be $\lfloor \frac{M}{k+1}\rfloor$, except for the last n-gram, which may be slightly longer if Mmod(k + 1) ≠ 0. Implementation-wise, each n-gram is represented as a view on a portion of a UMI sequence, to avoid creating unnecessary copies of the UMI sequence data.

When initializing the hash table of UMI n-grams, the time complexity for N UMIs is simply O(MN), since each n-gram of each UMI sequence must be extracted and hashed. To estimate the expected run time of each remove_near query, we assume that the UMIs are uniformly sampled across the sequence space Σ^M. Then, for each possible n-gram segment location, there are around ${\displaystyle N\frac{{\left|\Sigma \right|}^{M-\frac{M}{k+1}}}{{\left|\Sigma \right|}^M}=\frac{N}{{\left|\Sigma \right|}^{\frac{M}{k+1}}}}$ different UMI sequences at each n-gram location, assuming M is divisible by k + 1. Therefore, the time complexity of each query is O(M + kN|Σ|^−M∕k) on average, since each bin of the k + 1 n-grams of the queried UMI must be examined. This algorithm scales very well with longer UMI lengths, since the number of possible n-gram combinations increases, which decreases the likelihood of two UMI sequences having the same n-gram and allows each n-gram bin to be smaller.

Subsequences

Another method for decomposing a UMI sequence into a more general representation during both initialization and queries is by extract subsequences of length M − k (SymSpell, 2019). Alternatively, this can be thought of as picking any k nucleotides in the UMI, and replacing all of them with a new placeholder character (i.e., a subsequence of “ATCG” is “A⋆CG” if k = 1 and the placeholder is ⋆). It is easy to see that if two subsequences with exactly k placeholder characters exactly match, then their corresponding UMIs must match within k edits. After a hash table that is indexed by these subsequences is built in $O\left(N\left(\genfrac{}{}{0.0pt}{}{M}{k}\right)\right)=O\left(M^{k}N\right)$ time, we only need to examine the $\left(\genfrac{}{}{0.0pt}{}{M}{k}\right)$ bins of UMIs that correspond to each of the subsequences of the queried UMI during a query, in O(M^k + R) time, to get R UMI results. Note that some UMIs may have to be examined multiple times during a single query, since each UMI is placed in $\left(\genfrac{}{}{0.0pt}{}{M}{k}\right)$ different bins when building the hash table. This is faster than the combinations algorithm because it spreads out the work of generating sequence combinations within k edits over the initialization and the queries phases of the algorithm. This method is implemented in Calib (Orabi et al., 2018) and referred to as “locality-sensitive hashing”.

 
__________________________________________________________________________________________________________ 
Algorithm 4 The remove_near operation for the BK-tree data structure. 
___________________________________________________________________________________________________________ 
  procedure remove_near(u, k, F) 
     return dfs_bk_tree(root, u, k, F)        ⊳ Start at the root node of the BK-tree. 
  end procedure 
  procedure dfs_bk_tree(c, u, k, F) 
     r ←{}                                       ⊳ Resulting list of removed UMIs. 
     Δ ← d(c,u) 
     if c is not removed and Δ ≤ k and f(c) ≤ F then 
         Remove c. 
         r ← r ∪ c 
     end if 
     for each child v of node c where Δ − k ≤ d(c,v) ≤ Δ + k do      ⊳ Pruning by 
  Hamming distance. 
      if ∃w in the subtree of v, where w is not removed then    ⊳ Pruning by removed 
  nodes 
          if ∃w in the subtree of v, where f(w) ≤ F then        ⊳ Pruning by frequency 
  threshold F 
                  Add each element in dfs_bk_tree(v, u, k, F) to r. 
            end if 
          end if 
     end for 
     return r 
  end procedure 
________________________________________________________________________________________________________________    

BK-Tree

The BK-tree (Burkhard & Keller, 1973) is a type of metric tree that can make use of constant-time Hamming distance computations for handling remove_near queries. Each node in a BK-tree represents an UMI sequence, and each child node in a BK-tree is indexed by the Hamming distance between that child node’s UMI and its parent node’s UMI in an array in the parent node. When inserting an UMI, if a parent node v already has a child node at a specific Hamming distance d(u, v) for some inserted UMI u, then the insertion operation is recursively continued with the child node as the new parent node. Otherwise, a new node with the u is created and attached to the parent node at the index d(u, v) in parent node’s array. A fully built BK-tree is shown in Fig. 2.

When querying, we recursively visit a subset of nodes in the BK-tree. For each node we visit starting from the root, we find the edit distance d(u, v) between the queried UMI u and the UMI at current node v. Then, we only visit the 2k + 1 children of the node of v that are indexed by some edit distance i, where d(u, v) − k ≤ i ≤ d(u, v) + k. This is due to the triangle inequality that mandates that ${\displaystyle \begin{array}{c}\hfill d\left(u,w\right)+d\left(v,w\right)\ge d\left(u,v\right),\\ \hfill d\left(u,v\right)+d\left(u,w\right)\ge d\left(v,w\right)\end{array}}$ for some node w in the subtree of v which is considered as a candidate for the query’s result. Note that i = d(v, w). Rearranging, we get ${\displaystyle \begin{array}{c}\hfill d\left(u,w\right)\ge d\left(u,v\right)-d\left(v,w\right),\\ \hfill d\left(u,w\right)\ge d\left(v,w\right)-d\left(u,v\right)\end{array}}$ Since we only care about a candidate node w if k ≥ d(u, w), we can substitute and arrive at ${\displaystyle \begin{array}{c}\hfill k\ge d\left(u,v\right)-d\left(v,w\right),\\ \hfill k\ge d\left(v,w\right)-d\left(u,v\right)\end{array}}$ which can be rearranged to obtain ${\displaystyle \begin{array}{c}\hfill i=d\left(v,w\right)\ge d\left(u,v\right)-k,\\ \hfill i=d\left(v,w\right)\le d\left(u,v\right)+k\end{array}}$ This allows us to restrict the number of nodes visited during a query operation and speed up the query (Burkhard & Keller, 1973).

The time complexity of initializing and querying a BK-tree depends heavily on the height of the tree. On average, inserting all UMIs during initialization requires O(NlogN) time, since the height of the tree is around O(logN). For each query that results in R UMIs in total, the average time complexity is O(kR + klogN). The factor of k is due to extra children nodes that are examined while traversing the tree. BK-trees are efficient since the height of a BK-tree built on random UMIs is expected to be very low compared to the total number of unique UMI sequences.

If all nodes in a subtree are removed, then that entire subtree can be skipped during each remove_near query, which saves some time. To keep track of these subtrees, we can propagate a flag that indicates whether an entire subtree is completely removed, from child nodes to parent nodes (i.e., a node’s subtree is completely removed iff that node is removed and all of its children subtrees are completely removed).

Similarly, if all nodes in a subtree have a UMI frequency greater than the threshold F, then that subtree does not need to be examined in the query. The minimum frequency of each subtree can be stored and updated after each remove_near operation to prune subtrees from the search space (i.e., a node’s subtree’s minimum frequency is the minimum of that node’s UMI frequency if it is not removed, and the minimum frequencies of each of the node’s children subtrees that are not completely removed). Subtrees that have a minimum frequency greater than F are pruned since they have cannot contain any node that has a UMI frequency less than or equal to the threshold F.

Typically, nodes are inserted in arbitrary order during initialization. However, to ensure that more subtrees are pruned by their minimum frequencies, the UMIs can be inserted in increasing order, based on the frequency of the UMIs. By inserting higher frequency UMIs later, they end up closer to the leaf nodes of the tree, and allow lower frequency UMIs to end up closer to the root of the tree. This allows more subtrees to be pruned by the frequency threshold F.

A sketch of the full algorithm for handling remove_near queries with a BK-tree is presented in Algorithm 4.

Fenwick BK-trees

In remove_near queries, if the number of high frequency UMIs significantly outnumber the number of low frequency UMIs, it may be beneficial to first search for all UMIs that have a frequency lower than F, and then examine those UMIs to find ones that are within k edits of the queried UMI. This is an alternative approach for solving remove_near queries compared to other methods like directly using BK-trees (Burkhard & Keller, 1973), since we build an overarching data structure for pruning UMIs based on UMI frequencies instead of edit distance. Finding UMIs with frequencies lower than a threshold can be done with a Fenwick tree (Fenwick, 1994) on an array of BK-trees that is indexed by UMI frequencies. A Fenwick tree allows for fast queries on a prefix of an underlying array described by the Fenwick tree, as each node in the Fenwick tree describes a range of elements in the array, and only logN nodes are visited per prefix query. Figure 3 visualizes a Fenwick tree during a query operation. If the underlying array is indexed by the frequencies of the UMIs, then a prefix query up to index F on the Fenwick tree will return some property for all frequencies ≤F. In our case, each element in the underlying array is a BK-tree, and thus each node in the Fenwick tree also contains a BK-tree. The BK-tree of a parent node in the Fenwick tree will include UMI sequences from the BK-trees of each of the child nodes. When querying, BK-trees in the nodes of the Fenwick tree that represent the prefix up to F in the frequency array are independently traversed to find UMIs that are within k edits of the queried UMI.

To insert an UMI sequence, we must update logN BK-trees, each of height logN on average. Thus, initializing the Fenwick BK-trees data structure requires O(Nlog²N) time. Each remove_near query will take O(kR + klog²N) time, since logN BK-trees must be examined. However, in practice, we expect the heights of the BK-trees to be very short, and we expect this method to prune a significant portion of the UMIs based on their frequencies.

n-grams BK-trees

We propose a new method for faster remove_near queries based on a combination of BK-trees (Burkhard & Keller, 1973) and n-grams. The original n-gram algorithm requires a linear scan through each UMI that shares an n-gram with the queried UMI. We replace the bin of corresponding UMIs for each n-gram with separate BK-tree to speed up queries. This allows the height of each BK-tree to be very short, and also allows us to make use of methods for pruning subtrees, like ignoring subtrees that are completely removed and subtrees with minimum frequencies larger than F. Overall, this method is faster than using just the n-grams method on average, taking only ${\displaystyle \begin{array}{c}\hfill O\left(\right.M+k{R}_1+klog\left(\frac{N}{|\Sigma {|}^{\frac{M}{k}}}\right)\\ \hfill +k{R}_2+klog\left(\frac{N}{|\Sigma {|}^{\frac{M}{k}}}\right)+\dots +k{R}_{k+1}+klog\left(\frac{N}{|\Sigma {|}^{\frac{M}{k}}}\right)\left)\right.& \hfill \\ \hfill =O\left(M+kR+k^{2}log\left(\frac{N}{|\Sigma {|}^{\frac{M}{k}}}\right)\right)\end{array}}$ time, where the sum of the number of query results across each unique n-gram’s corresponding BK-tree is R (i.e., R₁ + R₂ + … + R_k+1 = R), and each O(log(N|Σ|^−M∕k)) is due to the height of the BK-tree corresponding to one of the k + 1 unique n-gram of length $\lfloor \frac{M}{k+1}\rfloor$. If a significant portion of the UMI sequences share the same n-gram, then the algorithms will run as fast as just using one large BK-tree for all UMI sequences, which takes O(kR + klogN) time, not O(N) time. This method is visualized in Fig. 4.