Performance question

[anderspitman]

First of all thanks for the excellent library.

I'm working on a program for copying certain references into a new bam file. Source here.

It's running slower than I would expect, only managing ~6MB/s at 100% CPU usage. I suspect I may be using noodles incorrectly. Any suggestions?

[zaeleus]

Hi @anderspitman,

When looking to optimize, I first recommend profiling with performance counters to better understand where time is being spent during execution. For example, on a Linux system using perf:

$ cargo build --release
$ perf record --call-graph dwarf target/release/samslice in.bam out.bam chr9 chrM
$ perf report --hierarchy --disassembler-style intel

-  100.00%        samslice
   -   95.12%        samslice
      +   74.28%        [.] miniz_oxide::deflate::core::compress_inner
      +    3.80%        [.] miniz_oxide::inflate::core::decompress
      +    3.47%        [.] miniz_oxide::deflate::core::compress_block
[...]

This shows that 80+% of time is being spent in the DEFLATE codec. Like htslib, linking to a optimized DEFLATE implementation for BGZF reading and writing is beneficial for this particular application. noodles-bgzf supports linking to libdeflate via a feature flag.

$ cargo add noodles-bgzf --features libdeflate

This is a simple change that provides a fair amount of improvement, assuming there is no requirement preventing you from linking to C libraries.

$ hyperfine --warmup 3 --parameter-list src in.bam --export-markdown results.md "target/release/samslice_miniz_oxide {src} out.bam chr9 chrM" "target/release/samslice_libdeflate {src} out.bam chr9 chrM"

Command	Mean [s]	Min [s]	Max [s]	Relative
samslice_miniz_oxide	5.545 ± 0.101	5.440	5.813	1.95 ± 0.08
samslice_libdeflate	2.846 ± 0.110	2.717	3.050	1.00

Depending on how performance-sensitive your application is, the query reader can be improved to reduce allocations. While bam::io::reader::Query does have an internal record buffer, the iterator implementation cannot self-reference it, which requires a copy per record returned. One workaround is to compose a query reader with a BAM reader and reusing an owned bam::Record. See #261 (reply in thread) for an example.

Also note that noodles does more comprehensive validation than other implementations upon serialization, which looks unfavorable in benchmarks at the expense of correctness. If this seems unnecessary, you may be interested in using rust-htslib instead.

[anderspitman]

My curiosity got the better of me so I ported samslice to rust-htslib. It was about the same speed as the noodles version. Interestingly, reusing a pre-allocated Record didn't speed it up. What did speed it up is increasing the number of threads. Went from ~90s to ~6s on my 32-thread machine. Does noodles have any threading functionality?

One downside is the executable is much larger (~10MB vs ~1MB).

Implementation here: https://github.com/anderspitman/samslice-htslib

[zaeleus]

Does noodles have any threading functionality?

Yes, bgzf::MultithreadedReader and bgzf::MultithreadedWriter are available to be used in advanced situations. For BAM, you will also have to read the index yourself and use it when querying. E.g.,

use std::{env, fs::File, io, num::NonZero, thread};

use noodles::{
    bam::{self, bai},
    bgzf,
};

// ...

let worker_count = thread::available_parallelism().unwrap_or(NonZero::<usize>::MIN);

let index_src = format!("{src}.bai");
let index = bai::read(index_src)?;

let mut reader = File::open(src)
    .map(|f| bgzf::MultithreadedReader::with_worker_count(worker_count, f))
    .map(bam::io::Reader::from)?;

let header = reader.read_header()?;

let query = reader.query(&header, &index, &region)?;

let mut writer = File::create(dst)
    .map(|f| bgzf::MultithreadedWriter::with_worker_count(worker_count, f))
    .map(bam::io::Writer::from)?;

for result in query {
    let record = result?;
    writer.write_record(&header, &record)?;
}

[anderspitman]

Awesome, thanks again!

[anderspitman]

Is there a way to use a single worker pool for read/write, similar to how rust-htslib does it, or do you always have to specify the number of workers separately?

[zaeleus]

No, the encode and decode thread pools are separate. They currently cannot be combined.