Chunking for Fasta Parsing

This is a quick note on some thoughts & experiments on fasta parsing, alongside github:RagnarGrootKoerkamp/fasta-parsing-playground.

It’s a common complaint these days that Fasta parsing is slow. A common parser is Needletail (github), which builds on seq_io (github). Another recent one is paraseq (github).

Paraseq helps the user by providing an interface for parallel processing of reads, see eg this example. Unfortunately, this still has a bottleneck: while the users processing of reads is multi threaded, the parsing itself is still single threaded.

1 Minimizing critical sections Link to heading

In Langmead et al. (2018), some approaches are discussed for FASTQ parsing.

1. Deferred parsing: search the end of the current block of 4 lines inside a critical section, and then parse those in separate threads.
2. Batch deferred: read \(4N\) lines, and then parse those. This reduces contention on the locks.
3. Block deferred: read \(B\) bytes, assuming that this exactly matches a single record. This can be done by padding records with spaces to give them all the same size, but is more of a hypothetical best case.

2 RabbitFx: chunking Link to heading

RabbitFx (Zhang et al. 2023) (github) is a recent parser written in C++ that does something new: it splits the entire input file into fixed size chunks (of 8 MB by default) that are distributed those over the threads. Reads that cross chunk boundaries store a link to the next chunk:

The RabbitFx paper then reimplements some basic algorithms on top of this parser and should significantly higher throughput, reaching up to 2.4 GB/s¹ for counting bases on a server with two 20-core CPUs, compared to FQFeeder’s 740 MS/s.

Unfortunately, all of this doesn’t work well with gzipped input which requires sequential unzipping. But anyway, I was curious to see how much perf this would give in practice. This can be used as an argument why single-block gzipping is bad.

3 Different chunking Link to heading

In my 1brc post, I took a different approach. When there are \(t\) threads, simply partition the input file into \(t\) non-overlapping chunks, and give one chunk to each thread. Then there might be some reads that cross chunk boundaries, but this is not a problem: each thread first finds the start of the first record (the first > character), and then processes things from there. Each thread ends processing when the next > character it runs into is in the next chunk (or the end of the file is reached).

In this model, the parsing does not require any synchronization at all!

In the code, I simply scan for all > characters and then collect a batch of reads spanning at least 1 MB. Then I send this to a queue to be picked up by a worker thread.² In the worker thread, I search the first newline inside each record for the end of the header, and the rest is the read itself.

Note that the underlying file is memory mapped, so that accessing it in multiple streams is trivial from the code perspective.

4 Experiments Link to heading

I did some experiments both on the CHM13v2.0.fa human genome, and a set of simulated reads³ from zenodo by Bede Constantinides that are used for benchmarking Deacon (github).

See the github code for the implementation.

rsviruses17900.fa is 5.4 GB and has reads of average length 2000 bp, and each read is a single line (ie no linebreaks inside records).

Adding RabbitFx to these benchmarks is still TODO, since I don’t feel like messing with C/C++ bindings at the moment.⁴

Needletail nor paraseq only allocate a newline-free record on request, which I don’t do in this benchmark, but I believe at least needletail does unconditionally store the positions of all newlines, and so does slightly more work.

The new method caps out at 17 GB/s, while the peak memory transfer rate of my DDR4 RAM⁵ is around 23 GB/s, so this is not too bad. If I drop the fixed 2.6 GHz CPU frequency that I used for the experiments, instead I get speeds up to 22.7 GB/s, which is actually very close!

Thus, if the throughput of your code on all threads is more than 2-4 GB/s, you might try this technique to avoid file parsing being the point of contention.

5 Outlook Link to heading

I think this is a neat idea, but probably with limited practical usage, since most large inputs will be zipped. And in that case, you may have many input files, and it’s better to spin up one thread/program per input file, so that decompression can be done in parallel.

6 Some helicase benchmarks Link to heading

See https://github.com/imartayan/helicase

Running the helicase benchmarks: cd bench && cargo run -r <path>:

Input files are taken/modified from https://zenodo.org/records/15411280:

Benchmarks take different inputs:

• path takes a file path and reads the file.
• mmap takes a file path and mmaps the file.
• slice and reader take a &[u8] from an already-read file.

Observations:

• Helicase is much faster for fastq than for fasta, surprisingly, both for &str and columnar/packed output.
• For fastq, helicase is best.
• For short reads, helicase is best (less per-record overhead).
• For fasta with newlines, helicase is slower than on no-newline data, but faster than needletail/paraseq.
• Needletail is quite good, but slow on short reads and human data with newlines.
• Paraseq is only comparable in speed for long chromosomes.
• In all cases, processing an in-memory &[u8] is faster than file input.
• mmap is faster than usual file parsing, but usually (not always) slower than bytes from memory.

We can modify our mmap-based multithreading to call helicase (slice input, string output) instead of our stub parser in each thread, as in this code. Note that we simply split the input into 6 chunks (one for each thread) and ignore cross-chunk reads in helicase.

Recall that our baseline simply searches for > characters and then for the first following \n to determine the end of the header, and does nothing else.

Observations:

• Our trivial method is already saturated at 3 threads.
• Helicase needs the full 6 threads for reads and viruses.
• Helicase is quite a bit slower on the line-breaked human data, probably because it does a lot of additional (and mostly inevitable) I/O.
• Needletail is quite fast in all single-threaded cases (but note that we do not call .seq()).

References Link to heading