Fuzzysplit: demultiplexing and trimming sequenced DNA with a declarative language

Availability

Fuzzysplit is implemented in Java without any external libraries. Its compiled binaries and source codes are available here: https://github.com/Daniel-Liu-c0deb0t/Java-Fuzzy-Search. It is available under the MIT license.

Fuzzy patterns

As each fuzzy pattern may have multiple subpatterns, multiple fuzzy string searches are required to find the matching locations of all subpatterns in a fuzzy pattern.

Fuzzysplit uses the following well-known recurrence to search for a string subpattern P of length |P| in some text T of length |T| with the Levenshtein distance metric (Levenshtein, 1966): (1) $$d(i,0)=0$$ (2) $$d(0,j)=j$$ (3) $$d(i,j)=\{\begin{array}{l}d(i-1,j-1),\text{\hspace{0.17em}if\hspace{0.17em}}{T}_i=P_j\hfill \\ 1+\min \left\{d(i-1,j-1),d(i-1,j),d(i,j-1)\right\},\text{\hspace{0.17em}otherwise}\hfill \end{array}$$ where d(i, j) indicates the Levenshtein edit distance between P_{1 … j} and T_{1 … i} for a match in T that ends at index i (Sellers, 1980). This property allows us to find match locations within a certain number of edits, where substitutions, insertions, and deletions all cost one edit. With dynamic programming, finding all matches runs in O(|P| · |T|) time. Fuzzysplit supports allowing transpositions to cost one edit (Wagner & Lowrance, 1975) and allowing custom wildcard characters by using variations on the algorithm above. Furthermore, Fuzzysplit allows partial overlaps between each subpattern and the ends of the text. These options allow the matching algorithm to be fine-tuned for specific use cases.

There are many ways to improve the efficiency of the algorithm above, especially with bit-vectors (Myers, 1999; Hyyrö, Fredriksson & Navarro, 2005; Wu & Manber, 1992b). We use Ukkonen’s cutoff heuristic (Ukkonen, 1985) to improve the average time complexity to O(k · |T|), where k is the maximum number of edits allowed. We choose this instead of the faster bit-vector techniques in order to retain length information about the matches. This technique is also used by Cutadapt (Martin, 2011).

For fuzzy searching with Hamming distance (Hamming, 1950), Fuzzysplit uses the Bitap algorithm (Baeza-Yates & Gonnet, 1989; Wu & Manber, 1992b) to allow matching with the Hamming distance metric. It has a worst-case time complexity of $O\left(k\cdot \left|T\right|\cdot \lceil \left|P\right|/w\rceil \right)$, where w is the length of a computer word. We use 63-bit words in Fuzzysplit.

We adopt the use of n-grams to quickly eliminate subpatterns before attempting one of the more time-consuming algorithms above. Fuzzysplit splits each subpattern into overlapping n-grams. If the text does not contain any n-grams of a specific subpattern, then we do not need to match that subpattern. However, for shorter subpatterns, edits within the text may cause none of the n-grams to match. Therefore, we only use the n-gram method to eliminate a subpattern P if its length satisfies the following inequality: (4) $$\left|P\right|>(k+1)(n-1)+k$$ where n is the n-gram size. The inequality represents the worse-case scenario where the edits are evenly spaced in the text and they split the text into segments that are just one short of the n-gram size. It is easy to see that if the inequality is satisfied for a subpattern, then it must have at least one n-gram that matches with a n-gram from the text if they are supposed to match within an edit distance of k by the Hamming (1950) or Levenshtein (1966) metric. This equation only covers Hamming and Levenshtein distance, but extending it to consider transpositions is trivial.

By default, Fuzzysplit automatically chooses the maximum n-gram size possible for all subpatterns in order to maximize the number of subpatterns eliminated. This is because the frequency of each unique n-gram is non-increasing as the n-gram size gets longer. For a subpattern P, the maximum n-gram size possible can be obtained by solving the inequality above for n: (5) $$n<\frac{\left|P\right|-k}{k+1}+1$$

We use the minimum of all such maximum n-gram sizes in order to satisfy the n-gram size constraint for all subpatterns. In Fuzzysplit, n-grams should not be used if there are wildcard characters within the subpatterns or the text, as they are not considered when generating n-grams.

Fixed-length wildcard patterns

Searching for a fixed-length wildcard pattern in some text T is trivial in O(|T|) time by simulating a sliding window the length of the pattern and checking if the characters in the window are all part of the pattern.

Matching multiple interval-length wildcard patterns

Unlike matching fixed-length wildcard patterns and fuzzy patterns, we match multiple consecutive interval-length wildcard patterns at the same time. For a list of interval-length wildcard patterns $\mathcal{P}$ of length $\left|\mathcal{P}\right|$, we wish to match each pattern ${\mathcal{P}}_i\in \mathcal{P}$ to contiguous, non-overlapping segments of the text T, in order. Each interval-length wildcard pattern ${\mathcal{P}}_j$ only matches some segment of text T_{a … b} if min_j ≤ b − a + 1 ≤ max_j, and every character in T_{a … b} is present in the set of allowed characters for ${\mathcal{P}}_j$ (i.e., ∀i ∈ {a…b}, the character T_i ∈ A_j, where A_j is the set of allowed characters of ${\mathcal{P}}_j$). For example, the j-th interval-length wildcard pattern represented by "A-Z", where min_j = 0 and max_j = 3, matches "ABC" and ""(empty sequence), but not "123" and "ABCD", as it matches any sequence of uppercase letters of some length between zero and three.

We formulate the solution to the problem of matching multiple interval-length wildcard patterns $\mathcal{P}$ with a segment of text T as a recursive function. We define the recursive function as f(i, j), which returns a nonzero value iff ${\mathcal{P}}_{1\dots j}$ matches some T_{1 … i}: (6) $$f(0,0)=1$$ (7) $$f(i,j)={\displaystyle \sum _{k=i-{\text{max}}_j}^{i-{\text{min}}_j}}f(k,j-1)\cdot [T_c\in A_j,\forall c\in \{k+1\dots i\}]$$ where [F] evaluates to 1 for some expression F if F is true. Otherwise, it evaluates to 0. The answer to whether the entire ${\mathcal{P}}_{1\dots \left|\mathcal{P}\right|}$ matches the entire T_{1 … |T|} is given by whether $f(|T|,|\mathcal{P}|)\ne 0$.

The intuition behind the algorithm is that it attempts to continue a match for a previous pattern by matching the current pattern right after the previous pattern’s match location in the text. A continuation is only possible if there is a valid number of characters in the text after the previous pattern and all of those characters are allowed by the current pattern. More formally, for each i ∈ {1 … |T|} and $j\in \{1\dots |\mathcal{P}\left|\right\}$, we assume that the current pattern ${\mathcal{P}}_j$ must end at i in the text. That implies that ${\mathcal{P}}_j$ must start at some index k + 1 where i − max_j ≤ k ≤ i − min_j in order to satisfy the length restriction for interval-length patterns. Furthermore, each character in T_{k + 1 … i} must be allowed by the current pattern ${\mathcal{P}}_j$. This whole transition relies on calculating the answers for the previous pattern ${\mathcal{P}}_{j-1}$ that ends at different indexes in the text, which is essentially solving the exact same problem. This eventually leads to the base case where an empty list of interval-length patterns matches an empty piece of text.

A naive implementation using dynamic programming to cache intermediate function calls results in a loose upper bound of $O(|T{|}^3\cdot \left|\mathcal{P}\right|)$ for the time complexity, which is very inefficient. We improve the time complexity by using two methods:

• We can speed up the check for whether a pattern allows all of the characters in a segment of text by precomputing the longest length of contiguous characters that ends at a certain index and only has characters that are allowed by a certain pattern. We show this value for each interval-length pattern ${\mathcal{P}}_j$ and ending index i as l(i, j), where (8) $$l(0,0\dots |\mathcal{P}|=0$$ (9) $$l(i,j)=\{\begin{array}{l}l(i-1,j)+1,\text{\hspace{0.17em}if\hspace{0.17em}}{\mathcal{P}}_j\text{\hspace{0.17em}allows\hspace{0.17em}}{T}_i\hfill \\ 0,\text{\hspace{0.17em}otherwise}\hfill \end{array}$$

• We notice that the overlap between the segment in the text where the current pattern ${\mathcal{P}}_j$ has to start its match ([i − max_j + 1, i − min_j + 1]) and the segment where the text characters are all allowed by ${\mathcal{P}}_j$ ([i − l(i, j) + 1, i]) is the only region in the text where ${\mathcal{P}}_j$ can begin, if it must end at index i. That means that we only need to query a consecutive region for whether the previous pattern that ends anywhere in that region matches the text. For this task, we keep separate prefix sums of f(i, j) for each interval-length pattern: (10) $$p(0,0)=1$$ (11) $$p(i,j)=p(i-1,j)+f(i,j)$$

Finding the sum of some segment f(a … b, j) for some pattern ${\mathcal{P}}_j$ becomes p(b, j) − p(a − 1, j), an O(1) time operation.

The optimizations above directly lead to the following recursive solution: (12) $$f(i,j)=p(i-{\text{min}}_j,j-1)-p(\text{max}\{i-{\text{max}}_j,i-l(i,j)\}-1,j-1)$$

We implement this in Fuzzysplit for a total run time complexity of $O(|\mathcal{P}|\cdot |T|)$ by using $\left|\mathcal{P}\right|\times \left|T\right|$ dynamic programming matrices for l, p, and f.

Overarching algorithm

Fuzzysplit uses a greedy overarching algorithm that matches a list of arbitrary patterns ${\mathcal{P}}_{1\dots \left|\mathcal{P}\right|}$ from one line of the template file with its corresponding line of input text T_{1 … |T|}. First, it partitions the list of patterns into continuous chunks of either fixed-length patterns (both fuzzy patterns and fixed-length wildcard patterns) or interval-length patterns. The result is adjacent chunks with alternating general pattern types. The general idea of the overarching algorithm is to search for candidate match locations for each entire chunk of fixed-length patterns and checking whether the interval-length patterns immediately before that chunk matches the text. We show a sketch of the algorithm in Algorithm 1. Note that the first and last patterns must be separately handled to take into account start and the end of the line of input text in the actual implementation.

To search for a matching location with an entire fixed-length pattern chunk, Fuzzysplit goes through each pattern in the chunk and uses each of them as the starting pattern until there is no more patterns left or if a valid match is found. The starting pattern is first searched by itself throughout the input text. For each match index of the starting pattern, Fuzzysplit first attempts to match the interval-length chunk before the current fixed-length chunk, and then it matches rest of the patterns in the fixed-length chunk that are after the starting index. If valid matches are found for both the interval-length chunk and the fixed-length chunk, then Fuzzysplit repeats the same algorithm with the next pair of interval-length and fixed-length pattern chunks. Note that if all the patterns in a fixed-length chunk are not required and none of the patterns actually match the text, then the interval-length chunks before and after the fixed-length chunk are merged into a larger interval-length chunk.

Unfortunately, the algorithm we describe is not efficient, as the matching operation for the interval-length chunk and the fixed-length chunk is repeated for each starting pattern and for each match location. We optimize the algorithm with two simple strategies. First, we cache the result of matching each fixed-length pattern at some location in the input text so they can be accessed in O(1) time without having to be matched again at that location. Then, we cache the dynamic programming matrix used for matching interval-length patterns, which can be extended without recalculating previously calculated dynamic programming values. This is possible because the start location in the input text to start matching the interval-length chunk is always the same no matter where the current fixed-length chunk matches, as the algorithm never backtracks on previous chunk matches. Only the end location of the interval-length pattern changes as different candidate match locations of the current fixed-length chunk are examined. Overall, the worst-case amortized time complexity for matching the interval-length chunk across all of the candidate match locations of the current fixed-length chunk is linear to the length of the input text.

The algorithm we describe greedily picks the first encountered valid match for the current interval-length and fixed-length chunks that is valid. This is efficient, as it can terminate early if a required pattern is not matched, but it does not lead to the most optimal match. In some cases, previous greedy matches of fixed-length chunks may cause later chunks to run out of space in the text. We choose to implement this greedy algorithm instead of an backtracking exhaustive search algorithm that finds the best possible match due to time complexity concerns.

We also use a greedy approach for matching within a chunk of fixed-length patterns after a starting pattern match candidate has been found for the chunk. The fixed-length patterns are matched one by one, without any backtracking. For fuzzy patterns, Fuzzysplit greedily selects the match with the lowest edit distance.

Multithreading

As reads can be handled independently of each other, parallelizing the matching process is very feasible with the producer-consumer paradigm. Fuzzysplit uses the main thread to read in batches of reads from all input files and place them in a queue. Meanwhile, it spawns multiple worker threads that handle batches of reads from the queue in parallel. Each thread handles one batch of reads (each read is handled independently of other reads, even within a batch) and then outputs them as soon as it is done. The output step is synchronized to allow only one thread to write out to some file at a time. Batching reads allow a reduction in output synchronization at the expense of increase memory usage. Furthermore, the main thread may block to wait for the worker threads by using a semaphore, in order to constrain the amount of reads stored in memory at one time.

Adapter trimming

We test Fuzzysplit on the simple task of searching for 3′ adapters. The input data is one FASTQ file with 10 million reads. Around 75% of the reads have a randomly generated 3′ adapter that is 30 nucleotides long. Each read’s adapter has a 50% chance of having no edits, and 50% chance of having one edit. In total, each DNA sequence is padded with random nucleotides to around 130 nucleotides long, with or without the adapter. We generate two sets of data, one with substitution edits only, and another with insertions, deletions, and substitutions. The data simulation scripts are available here: https://github.com/Daniel-Liu-c0deb0t/Java-Fuzzy-Search/tree/master/scripts.

We measure the run times for a few different settings and we show them in Fig. 1. Each run time is relative to the baseline at 100%. The baseline is mostly the default settings in Fuzzysplit: allowing up to one edit based on the Levenshtein (1966) metric, using a single thread, and using the default adaptable n-gram setting (in this case, it uses an n-gram size of 15 due to the adapter length being 30 nucleotides long). Each setting shown in Fig. 1 is a modification relative to the baseline. Also, all tests were ran on the data set with insertions, deletions, and substitutions except for the run that uses the Hamming (1950) metric. Note that Fuzzysplit perfectly found all reads with adapters.

We find that increasing n-gram sizes does indeed lower the run time, as it is better at eliminating subpatterns before the Levenshtein (1966) edit distance is calculated. Running the program with two threads does improve the run time. We do not show more than two threads because the run time gets worse due the testing machine being limited to only two CPU cores. However, we expect the run time to continue to improve with multithreading if more cores are available. Finally, as expected, matching with the Hamming (1950) metric drastically improves the run time.

Demultiplexing

We first generate 48 random barcodes of varying length between 8 and 15 nucleotides long. Then we generate 10 million reads, where each read’s DNA sequence has a 50% chance of having one of the 48 5′ barcodes. Overall, each DNA sequence consists of a possible barcode sequence and 100 random nucleotides. Each barcode can have up to one insertion, deletion, or substitution edit, similar to how the adapters were generated.

We test the run time for Fuzzysplit, Cutadapt (Martin, 2011), and GBSX (Herten et al., 2015), and we show the results in Fig. 2. Fuzzysplit was ran using two threads, the default n-gram setting, and allowing one edit with the Levenshtein metric. For Cutadapt (Martin, 2011), each barcode was anchored to the 5′ end and we used 0.125 as the edit threshold for Levenshtein distance. GBSX (Herten et al., 2015) was also ran with allowing one edit with Levenshtein (1966) distance, but checking for enzymes and 3′ adapters were disabled to match the other tools. Overall, Fuzzysplit is slower than the other tools since it trades speed for much greater flexibility.

The difference percentages for all three tools are shown in Fig. 3, where each difference D is defined as (13) $D=\frac{1}{\text{with}\_\text{barcode}}{\displaystyle \sum _b}{\text{wrong}}_b+{\text{missing}}_b$

The difference metric counts both the number of reads that are missing from a file and the number of reads that are not supposed to be in that file, and divides the sum by the total number of reads that have a barcode. Overall, the difference percentage of Fuzzysplit is comparable to that of Cutadapt (Martin, 2011), and GBSX (Herten et al., 2015) is less accurate than either tools. Note that Cutadapt and Fuzzysplit both resulted in more false positives than false negatives over all barcodes, while the opposite was true for GBSX.