Use ISA-L even in the Rust version

[marcelm]

Commit 759398f added support for reading gzip files with ISA-L to the C++ code. This has not been ported to Rust, yet.

We currently use the flate2 crate for decompression. We use their zlib-rs backend, which is supposed to be the fastest one.

Decompressing compressed FASTQ files represents only a small fraction of total runtime, so in many cases we don’t actually gain that much, but it is highly relevant when we try to run strobealign on a machine with many cores. Before we had ISA-L in the C++ code, strobealign could only saturate about 32 CPU cores. When more cores are available, the worker threads sit idle part of the time because they have to wait for the (single) decompression thread to provide them with data.

I measured this on the Rust version and can see the same thing: When increasing the number of threads, CPU usage scales accordingly up to about 32 threads and then stays flat.

When I tried to switch our code to use the isal-rs crate, it initially looked straightforward and was essentially just a cargo add isal-rs and switching from flate2::read::MultiGzDecoder to isal::read::GzipDecoder. However, the code does not actually compile because GzipDecoder does not implement the Send trait.

This is a Rust-specific thing: The Send trait marks types that are safe to send from one thread to another. We do this in strobealign because the GzipDecoder is created in the main thread so that we can estimate the read length, and then it is moved into a new thread that reads the entire file.

Here are a couple of ways to solve this:

• Explicitly mark isal::read::GzipDecoder as Send. This requires the use of the unsafe keyword and is very, very likely not the proper way to solve this. We should do this only if we understand much better what the implications are.
• Try cloudflare’s fork of zlib instead, which is also supposed to be a lot faster than regular zlib.
• Restructure the code so that we don’t need to send the GzipDecoder to a different thread. That is, read the initial 500 reads for estimating read length also in the reader thread.

I have tried the first option just to have something that I can perhaps use in benchmarks in order to see whether it actually makes a difference in core usage, but there were some issues compiling it.

I haven’t tried cloudflare’s zlib, but from what I heard, it is not as fast as ISA-L.

Restructuring the code is probably the cleanest way forward.

See also milesgranger/isal-rs#33

[marcelm]

I’ve implemented the first option (the unsafe one that just marks isal::read::GzipDecoder as Send), but the benchmarks show only a small improvement. I have also tested decompressing the input in a separate process using <(igzip -dc input.fastq.gz) as the "filename" for the reads, but even with that, CPU usage does not go above 3200%, so I’m no longer so sure that decompressing the input is the bottleneck.

I’ve done a few measurements in alignment mode on a dataset with 100 million paired-end reads. The runtimes below are reported in microseconds per read (not per read pair).

What	Runtime in µs/read or µs/record
Mapping (alignment)	80
Parsing FASTQ	0.4
Decompression with flate2 plus parsing FASTQ	1.24
Decompression with isal-rs plus parsing FASTQ	0.84 µs

flate2 is what we currently use for decompression.

We have one "reader" thread that decompresses the input, parses the FASTQ records and splits them into chunks of roughly 1 million nucleotides each. Given the above numbers, we should even with the current flate2 solution be able to decompress so quickly that about 80 / 1.24 = 64.5 workers have something to do, but in reality, not more than about 35 cores are used no matter how many worker threads run.

I’ll try to see what is going on by profiling, but this may be a bit hard because I only have 8 cores and the cluster with the 128 core machines doesn’t let me use the perf profiler (needs certain permissions that users typically don’t have).

[marcelm]

I made some progress on this.

First, I temporarily changed the worker thread to not do any mapping. It just receives the data and discards it. This let me get a profile on my 8-core CPU. I was able to see that the writer thread is not the issue.

Then I found out that it is possible to use perf on the cluster after all, but not samply (which I had wanted to use). After some fiddling, I got this profile (just for the reader thread) in the ISA-L version (with 64 workers):

  Children      Self  Command          Shared Object             Symbol                                                                                  
+   30.37%     1.41%  strobealign-rs-  strobealign-rs-isal-prof  [.] alloc::raw_vec::RawVecInner<A>::finish_grow                                         
+   29.73%     0.00%  strobealign-rs-  [unknown]                 [.] 0x0000000000000001                                                                  
+   27.75%     0.22%  strobealign-rs-  strobealign-rs-isal-prof  [.] _mi_malloc_generic                                                                  
+   27.63%     0.07%  strobealign-rs-  strobealign-rs-isal-prof  [.] mi_heap_malloc_zero_aligned_at_generic                                              
+   26.18%     1.41%  strobealign-rs-  strobealign-rs-isal-prof  [.] _mi_heap_delayed_free_partial                                                       
+   25.53%    25.50%  strobealign-rs-  strobealign-rs-isal-prof  [.] _mi_page_free_collect                                                               
+   24.77%     0.11%  strobealign-rs-  strobealign-rs-isal-prof  [.] _mi_free_delayed_block                                                              
+   24.72%     6.77%  strobealign-rs-  strobealign-rs-isal-prof  [.] <strobealign::fastq::FastqReader<R> as core::iter::traits::iterator::Iterator>::next
+   14.25%    14.23%  strobealign-rs-  strobealign-rs-isal-prof  [.] _mi_page_malloc                                                                     
+   12.45%    12.44%  strobealign-rs-  strobealign-rs-isal-prof  [.] ..@42.end                                                                           
     7.35%     7.35%  strobealign-rs-  strobealign-rs-isal-prof  [.] ..@37.end                                                                           
+    5.23%     5.22%  strobealign-rs-  strobealign-rs-isal-prof  [.] std::io::read_until                                                                 
+    4.92%     4.66%  strobealign-rs-  libc-2.31.so              [.] __memcpy_avx_unaligned_erms                                                         
+    4.32%     0.00%  strobealign-rs-  [unknown]                 [.] 0000000000000000                                                                    
+    3.66%     0.00%  strobealign-rs-  [unknown]                 [.] 0xfb89485354415541                                                                  
+    3.57%     0.03%  strobealign-rs-  strobealign-rs-isal-prof  [.] <strobealign::io::UnsafeGzipDecoder<R> as std::io::Read>::read
...

So the bottleneck is actually memory allocation. I cannot tell from the profile, but there are two places where memory is allocated: Within the FASTQ reader and within the reader thread when it creates chunks of SequenceRecords. There are other FASTQ readers that are optimized to do very few memory allocations, we should probably try one of those. They don’t quite have the same features as what we need, but maybe we can make it work.