Algorithms #2 — Mapping Fun! Super Short Sequence Mapper

Side Project Notes: CubaTrie! A Super Short Sequence Mapper utilising a Radix Trie Data Structure

Source code can be found on my github page! I also published this on my substack — which I will include ad-hoc results that are not included in the final manuscript.

CubaTrie is a C-based sequence mapper that indexes super short references or sequences, of length 20–30 basepairs, in a 4-bit radix trie and performs seed-and-extend matching on query sequences of similar or longer lengths (trimmed reads, paired-end sequencing data, long read sequencing, reference genome, etc.). CubaTrie quantifies mappings for downstream analysis. In this article, I will briefly explain the Data Structure Algorithms (DSAs) implemented in CubaTrie to map super short sequences blazingly fast (relative to state-of-the-art tools) and a couple of interesting DSA employed to expedite mapping and resolve certain issues such as inexact matching in a Trie Search.

A brief background on how the entire project came about — the basic Trie Data Structure was initially written in Python back in June 2025 to solve a super short sequence mapping issue for CRISPR-Cas 9 pooled screening, mapping guide RNA (gRNAs) of length 20bp to paired end data. I extended it as a personal side project to improve its performance, in terms of speed, and implemented a couple of DSA including a Radix Trie structure, a compressed prefix trie, and a k-mer cacher based on substrings previously checked.

The new implementation was remarkably fast, even in Python, consequently I tried writing it in C. The first iteration in C returned significant improvements in terms of speed, it was >10–20x faster relative to Python’s version, even at 1 CPU Core. The time taken to analyse the same dataset dropped from 40 minutes in the first Pythonic version to under a minute in this C based unoptimised mapper. I left the project for the back burner in Q4 2025 though.

Fast forward 2026 April, I wanted to acquaint myself with agentic AI, which I found to be incredibly powerful, it’s like a super smart and efficient senior engineer. I was also trying to make sense of my role in these always-evolving A.I. developments, especially in the space of Bioinformatics — but more on that in a separate post! I was also leaving my current post and I wanted to wrap up the project since it displayed improved performance.

CubaTrie performed way better than expected, which consequently became a full blown tool with a couple of additional new features such as cut mode, inexact matching — allowing for indels too! This was in parallel as I’m wrapping my projects at my current work. I analyzed a couple of public dataset, compiled the results and wrote a manuscript (which I hope it goes through) all in a span of 4–5 weeks? I guess it’s a pretty good problem to have!

I figured I would write this article to document the process, the DSAs implemented to speed up mapping of super short sequences (I don’t write and read on DSA for work — it’s mostly out of self interest), and current limitations of CubaTrie. This article is not on CubaTrie evaluation against other conventional aligners in terms of its performance, including run time, memory, and alignments. This article also touches briefly on the basic DSA and dives straight to the essential details of the implementation.

An overview of the DSA implemented — CubaTrie constructs a Compressed 4-bit Radix Trie to store its reference, the super short sequences, builds a k-mer bitset, lookup for each substring (per base in the read) in the k-mer bitset, and extends to the trie traversal if a match is found. For inexact pattern matching, CubaTrie implements a Depth First Search (DFS) to store the number of errors (substitutions, and indels when enabled) during Trie Traversing, and prunes any recursion that exceeds the constraints (e.g. max number of errors).

4-bit Radix Trie

Trie or Prefix Tree is a Data Structure whose nodes store a single letter, either bases (Figure 1) or alphabet, where words or prefixes can be retrieved from the structure by traversing down a branch path of the tree, commonly used for autocomplete (think search engines or messaging on your phone), and IP routing by Network Engineers.

Figure 1. Trie or Prefix Tree Data Structure for 4 gRNAs: AGCTGTAG, AGCGGTAT, AGCGGTAA, ACCCGTAA. A trie stores strings base-by-base along a path where each node houses a single character or base. The } indicates which nodes are compressed into a single node as per figure 2.

Radix Trie instead implements path compression where nodes with single children are merged (Figure 2) thereby saving space via reducing the number of internal nodes and pointers. Each internal node has to have 2 or more children — or in the case of CubaTrie, a minimum or two and a maximum of 4. Unlike a regular trie, the edges or pointers of a radix trie can hold a sequence of a string, or just a single character element. Radix Trie is a space-optimized data structure with a space complexity O(NxL), where N is the number of super short references and L is the average length of those sequences. The lookup time complexity for a Radix Trie is O(k), where k is the length of the sequence being searched.

Figure 2. Radix Trie for the same 4 gRNA sequences of 8bp. Radix Trie implements path compression where nodes with single children are merged where each internal node must have 2 or more children. A node could contain substrings of various lengths. We observe that G-C is compressed into a single node housing substring GC, G-G-T-A is compressed to a single node housing substring GGTA. CubaTrie constructs the k-mer bitset based on the Radix Trie. Symbol ∧ marks the direct access of the Trie cursor states

Honestly, there isn’t much to add since there’s no new innovation or speedup technique except that it is a DNA-fixed fanout, specialized to A/C/G/T with an exact maximum of 4 child slots, and, also, that was how the name CubaTrie came about. Cuba (pronounced as Chew-bah and not the country Q-ba) means to try in the Malay Language. Essentially, the tool means Try Trie*. It is surprisingly common to hear the phrase ‘Cuba Try’ spoken amongst the Malays in Singapore and Malaysia. Why the double emphasis? It serves to add emphasis, increase friendliness, and create a colloquial, rhythmic flow, just like the tool to map for super short sequences.

Oh, there’s also the seed-to-cursor continuation (More in the next section). The bitset index stores trie cursor states that can resume inside a compressed edge rather than restarting from root. The tight coupling between compressed radix structure and seed index is what makes cubaTrie’s seed-and-extend relatively efficient.

Bitset

I initially implemented a hash-based cache set to store k-mers, k indicating the length of the substring typically of length 4–12 bases, that was found or searched in a previous scan which avoids expensive Trie Traversal at most positions. A naive approach performs the Trie search from every offset in every read. The hash-based cache set initial implementation led to a remarkable improvement in its speed.

I eventually figured out that we can retrieve the k-mers of interest a priori, instead of checking every offset and adding k-mers to a hash-based cache set, based on the constructed Radix Trie. The cache set evolved into a bitset (via Codex’s suggestion). A bitset is a compact array of bits where each position represents a yes/no state, with O(1) membership testing. If bit i is 1, item i or the specific k-mer of interest is present; if it is 0, it is absent. An example of a construction of a bitset for 4-mer, eg. TCAT.

Firstly, we convert the 4-mer to an integer code. With Start code = 0 and A=0, C=1, G=2, T=3:

1. Read T(0) : code = (0<<2)|3 = 3
2. Read C(1) : code = (3<<2)|1 = 13
3. Read A(2) : code = (13<<2)|0 = 52
4. Read T(3) : code = (52<<2)|3 = 211

Subsequently, we set the bit in words (bitset) for the integer code 211.

1. word_index = 211 >> 6 = 3
2. bit_offset = 211 & 63 = 19
3. words[word_index] |= (uint64_t)1u << bit_offset
4. words[3] |= (uint64_t)1u << 19

The 4-mer TCAT maps deterministically to integer code 211, which is then switched ‘on’ at the 3rd 64-bit word, bit position 19 (0-based).

CubaTrie converts each read substring (the k-mer) into a 2-bit encoded integer code (as shown above), checks for the seed presence in the bitset (0/1). If the seed is present, the bitset points to the state table to pass the trie cursor state, consisting of the node, edge_idx, and edge_off, where Trie Traversal continues or extends from that cursor state instead of the root node, skipping root-to-seed traversal.

A seed-and-extend with a bitset implementation is markedly faster in the case of CubaTrie as it avoids repeatedly starting the search from the radix trie root at every read offset. In a standard trie traversal, the first k bases must be re-traversed for each possible window, which adds substantial overhead across large FASTQ datasets or long genome references. Bitset seeding rapidly encodes the initial bases and rejects most non-matching windows using a simple bit test. For candidate windows that pass this filter, CubaTrie directs the search to a saved cursor state at depth k, bypassing the expensive root-to-seed traversal. This reduces pointer chasing, branching, and cache misses, resulting in fewer full trie traversals and much higher overall throughput.

There are limitations as for all seed-and-extend methods, one of which is that we may miss alignments where the mismatches are within the k-mer region. We’ve observed this in the case of aligning CasRx gRNAs to EV-A71 sequences (refer to the results on the publication). Bowtie maps slightly more than CubaTrie for cases where two mutations are found in the seed region for CubaTrie. It is possible to minimise the number of False misses via utilising shorter k-mers, the lowest that is permitted is 4, with a maximum mismatches of 1 (only substitution error permitted).

CubaTrie supports only one seed mismatch, via generation of Hamming-1 neighbours of the seed code (ie. change in one base at one position within the sequence), and testing those neighbour codes in the same bitset as explained above. If a neighbour is found, extension proceeds from that neighbor’s cursor state with one mismatch already accounted for. This design combines high sensitivity with speed via full trie traversal avoidance at most read offsets.

Why did I limit its seed mismatch to one? Increasing the number of mismatches allowed in the seed adds two major costs. While the bitset checks themselves are pretty cheap, there are more seed checks per window and many more candidate extensions. For example, in the case of an 8-mer, while exact mode test for only one 8-mer in a bitset, a 1 allowance for mismatch in seed adds another 24 (k*3 : 3 bases for every position in a k-mer) additional checks per window. Permitting mismatch seeds will naturally mean more windows pass pre-filter (including false positive or near positive hits), thereby extension runs more frequently and adds to the overall cost. It is especially costly in the case of tracking indels via DFS — more in the next section!

DFS for Inexact Matching

CubaTrie performs inexact trie traversal via a bounded Depth First Search (DFS). Briefly, Depth First Search is a traversal approach in which the traverse begins at the root node (empty circle at the top of Figure 2) and proceeds through the nodes till the terminal node (where gRNA object is stored) or as far as possible, depending on the constraints.

The DFS in CubaTrie starts from a seed-derived cursor state, explores one full path of edge labels in the radix trie as deep as possible, while tracking edit budget, before backtracking to try alternative branches. DFS is a natural way to enumerate all paths in the Radix Trie of the super short sequences, inclusive of those with or without mutations relative to the query sequence, while still controlling the search with strict bounds.

While path compression saves memory, the trie traversal behaves like an uncompressed trie, it checks for one nucleotide at a time. Each DFS state keeps a current trie cursor position, read position, reference-consumed length, edits used, and the operation path. From each state, the DFS state branches into a Match or Mismatch (costing 0 or 1 per base respectively), an insertion in the read (costing 1 per base) and a deletion from the reference sequence (costing 1 per base). The latter two are enabled optionally.

In CubaTrue, recursion paths are aggressively pruned once edit cost violates the length or mismatch constraints (e.g. edit distance > number of errors allowed). A hit is accepted when trie traversal reaches a terminal node within the configured reference-length window (max_len — min_len). Candidates nearby-hits are collapsed based on lowest edit distance, followed by those with fewer indels and a smaller span deviation. Otherwise, if all 3 conditions return the same values, the hit with the earlier start is printed. (On hindsight, perhaps it should print all — allows for further refinement in the downstream analysis)

DFS with indels enabled takes a longer time to execute due to the combinatorial growth. Since the DFS branches into insertion and deletion, it creates a much larger branching factor with many more alternative paths that are explored before pruning, relative to a substitution-only traversal.

Other small implementations

Instead of re-encoding each sliding window (O(nk)), CubaTrie implements a rolling-hash style, closer to O(n), by shifting out the old base and shifting in the new one for each step, where n is the read length. During read scanning, each read window encodes the k-mer and loops over all k bases again which is roughly (n-k+1) k base-to-bit operations. There is a 90% overlap with the previous k-Mer for a sliding window, which is repetitive and wasteful. While the initial approach is simple and robust, it spends extra CPU on repetitive encoding work since

Current Limitations

CubaTrie does have a couple of limitations that I am trying to figure out. Most urgent is the i/o — CubaTrie performs worse in terms of writing to a SAM format file or piping the output to samtools relative to Bowtie2 for the same number of CPU cores. The DSA definitely performs better than BWT for super short sequences as per the manuscript. Though, admittedly I would need to brush up on my multiprocessing knowledge — so that’s one to come!

CubaTrie also does not handle ambiguous bases (N) currently, non-ACGT bases are not represented in seed encoding or Radix Trie construction. nt2bits returns -1 for ambiguous bases (eg. N, R, Y, etc.), so those windows are skipped. Ambiguous bases are commonly found in Viral RNA genomes, which is one of the common query sequences for CubaTrie. These ambiguities arise primarily from high mutation rates and the presence of diverse viral populations (quasi-species) within a single host, making it difficult to assign a single base during consensus sequencing. I am still thinking whether it is worth the implementation since in the case of diagnostics and therapeutics we would rather the gRNA to target a conserved region of the RNA virus.

With regards to the mismatches for CubaTrie’s seed-and-extend strategy, only substitution (hamming-1) is allowed and the number of substitutions permitted is capped at 1 as iterated in a previous section. Current strategy implements a neighbour generation that changes only one base; it does not model for insertions and deletion in seeding. Seed sensitivity is limited for noisier data such as Long Read Sequencing where a higher error rate is prevalent.