Bug in modkit pileup with `--combine-mods` v0.6.4

First reproduce the reported command with the attached ak.sam.zip and --combine-mods, then compare it with the working single-modification commands. Inspect modkit-core/src/pileup/pileup_processor.rs at line 1894 and the worker failure reported at subcommand.rs:1659. Done means the combined-modifications run completes without hanging or panicking and produces pileup rows.

Dear ONT, I believe

that I've stumbled upon an bug in modkit pileup when --combine-mods flag is used.

On larger bam, the process froze totally doing nothing. On smaller bam, the process finishes with the second backtrace.

dist_modkit_v0.6.4/modkit pileup --modified-bases T:17802 A:a --combine-mods --ref ${reffasta} ak.bam ak.bed
> discarded 0 contigs with zero aligned reads
> parsed 2 base modification(s), combining modification counts together for primary sequence bases 'A,T'
> adding single-base motif: 'A 0'
> adding single-base motif: 'T 0'
> attempting to sample at least 10042 reads
> calculating threshold value with probabilities aligned to reference base(s) A,T
> setting threshold 0.5605469 (qual: 143) for base A, percentile 10, 597 total probabilites
> setting threshold 0.5761719 (qual: 147) for base T, percentile 10, 440 total probabilites
> using optimized workers for all-context
> 0 B written to output: ak.bed [0 B/s]
[00:00:00] ----------------------------------------       0/4215502 genome positions 0/s 0s
> 0 rows written
> 0 ~records errored                                                                                                                                                                      
thread '<unnamed>' panicked at modkit-core/src/pileup/pileup_processor.rs:1894:9:
STRANDS 2
stack backtrace:
   0:     0x5578e4748b32 - <std::sys::backtrace::BacktraceLock::print::DisplayBacktrace as core::fmt::Display>::fmt::hdcfcb6d4c8489523
   1:     0x5578e476cbb3 - core::fmt::write::h8a494366950f23bb
   2:     0x5578e4745953 - std::io::Write::write_fmt::h6556609fca33d0b1
   3:     0x5578e4748982 - std::sys::backtrace::BacktraceLock::print::hb2a626a81e06b2dc
   4:     0x5578e4749e93 - std::panicking::default_hook::{{closure}}::h4f78485264f12d10
   5:     0x5578e4749c75 - std::panicking::default_hook::h2c66fc99e962531d
   6:     0x5578e474a825 - std::panicking::rust_panic_with_hook::h33ac55f64bbd807d
   7:     0x5578e474a5da - std::panicking::begin_panic_handler::{{closure}}::h30e7cb89678a57fe
   8:     0x5578e4749039 - std::sys::backtrace::__rust_end_short_backtrace::hed60f27456c16ced
   9:     0x5578e474a26d - __rustc[de2ca18b4c54d5b8]::rust_begin_unwind
  10:     0x5578e3802eb0 - core::panicking::panic_fmt::h62f63d096dd276af
  11:     0x5578e3d08682 - mod_kit::pileup::pileup_processor::ACountsMatrix::incr_call::hb767530f0ff78e6b
  12:     0x5578e3cfee89 - <mod_kit::pileup::pileup_processor::DnaPileupWorker<T,M,_,_> as mod_kit::pileup::pileup_processor::PileupWorker>::process::hbead0a384f097cd2
  13:     0x5578e3bade03 - std::sys::backtrace::__rust_begin_short_backtrace::hb509df9362b458d6
  14:     0x5578e3bbef98 - core::ops::function::FnOnce::call_once{{vtable.shim}}::hbefd88c46c41be61
  15:     0x5578e474db9f - std::sys::pal::unix::thread::Thread::new::thread_start::hb6e99e73da4d28f8
  16:     0x150e0e82fb7b - start_thread
                               at ./nptl/pthread_create.c:448:8
> Done, processed 0 rows.
> 0 B written to output: ak.bed [0 B/s]                                                                                                                                                   
thread 'main' panicked at modkit-core/src/pileup/subcommand.rs:1659:34:
worker thread 2 paniced: Any { .. }
stack backtrace:
   0:     0x5578e4748b32 - <std::sys::backtrace::BacktraceLock::print::DisplayBacktrace as core::fmt::Display>::fmt::hdcfcb6d4c8489523
   1:     0x5578e476cbb3 - core::fmt::write::h8a494366950f23bb
   2:     0x5578e4745953 - std::io::Write::write_fmt::h6556609fca33d0b1
   3:     0x5578e4748982 - std::sys::backtrace::BacktraceLock::print::hb2a626a81e06b2dc
   4:     0x5578e4749e93 - std::panicking::default_hook::{{closure}}::h4f78485264f12d10
   5:     0x5578e4749c75 - std::panicking::default_hook::h2c66fc99e962531d
   6:     0x5578e474a825 - std::panicking::rust_panic_with_hook::h33ac55f64bbd807d
   7:     0x5578e474a5da - std::panicking::begin_panic_handler::{{closure}}::h30e7cb89678a57fe
   8:     0x5578e4749039 - std::sys::backtrace::__rust_end_short_backtrace::hed60f27456c16ced
   9:     0x5578e474a26d - __rustc[de2ca18b4c54d5b8]::rust_begin_unwind
  10:     0x5578e3802eb0 - core::panicking::panic_fmt::h62f63d096dd276af
  11:     0x5578e38032d6 - core::result::unwrap_failed::h95bc3f5a607b2c95
  12:     0x5578e3bcc335 - core::result::Result<T,E>::expect::hb714fbacf74e1ed4
  13:     0x5578e3beebb7 - mod_kit::pileup::subcommand::ModBamPileup::run::hf7ff9eb12f1f321d
  14:     0x5578e380fc1a - modkit::main::h132b333ca12c63c2
  15:     0x5578e3815766 - std::sys::backtrace::__rust_begin_short_backtrace::h9bb5aca9ca64d7de
  16:     0x5578e3813d22 - std::rt::lang_start::{{closure}}::h5a4f1a2efbc5fbb6
  17:     0x5578e473b160 - std::rt::lang_start_internal::h34f9328d113fd60a
  18:     0x5578e3810295 - main
  19:     0x150e0e7c6ca8 - __libc_start_call_main
                               at ./csu/../sysdeps/nptl/libc_start_call_main.h:58:16
  20:     0x150e0e7c6d65 - __libc_start_main_impl
                               at ./csu/../csu/libc-start.c:360:3
  21:     0x5578e38037a9 - _start
  22:                0x0 - <unknown>

There is no issue if the flag is not used

dist_modkit_v0.6.4/modkit pileup --modified-bases T:17802 --ref ${reffasta} ak.bam ak.bed
> discarded 0 contigs with zero aligned reads
> parsed 1 base modification(s). Base modifications other than 'T:17802' will be counted as 'N_other'.
> adding single-base motif: 'T 0'
> attempting to sample at least 10042 reads
> calculating threshold value with probabilities aligned to reference base(s) T
> setting threshold 0.5761719 (qual: 147) for base T, percentile 10, 440 total probabilites
> using optimized workers for all-context
> Done, processed 397 rows.

dist_modkit_v0.6.4/modkit pileup --modified-bases A:a --ref ${reffasta} ak.bam ak.bed
> discarded 0 contigs with zero aligned reads
> parsed 1 base modification(s). Base modifications other than 'A:a' will be counted as 'N_other'.
> adding single-base motif: 'A 0'
> attempting to sample at least 10042 reads
> calculating threshold value with probabilities aligned to reference base(s) A
> setting threshold 0.5605469 (qual: 143) for base A, percentile 10, 597 total probabilites
> using optimized workers for all-context
> Done, processed 539 rows.

I've drilled this down to 1 read that is causing this issue.

I'm attaching the read
ak.sam.zip

Should you need the .bam and reference sequences please reach out marek.schwarz AT biomed.cas.cz.

./dist_modkit_v0.6.4/modkit modbam check-tags ak.bam
> no errors
> num PASS records: 1 (100.00%)
> num records: 1
> valid record tag headers:
+------------+-------+
| tag_header | count |
+------------+-------+
| A+a.       | 1     |
| T+17802.   | 1     |
+------------+-------+


> modified bases:
+--------+--------------+----------+------+
| strand | primary_base | mod_code | mode |
+--------+--------------+----------+------+
| +      | T            | 17802    | .    |
| +      | A            | a        | .    |
+--------+--------------+----------+------+

Ps.: github suggest #562 as possible duplicate, however adding --use-dynamic does not help.

dist_modkit_v0.6.4/modkit pileup --use-dynamic --modified-bases T:17802 A:a --combine-mods --ref ${reffasta} ak.bam ak.bed
> discarded 0 contigs with zero aligned reads
> parsed 2 base modification(s), combining modification counts together for primary sequence bases 'A,T'
> adding single-base motif: 'A 0'
> adding single-base motif: 'T 0'
> attempting to sample at least 10042 reads
> calculating threshold value with probabilities aligned to reference base(s) A,T
> setting threshold 0.5605469 (qual: 143) for base A, percentile 10, 597 total probabilites
> setting threshold 0.5761719 (qual: 147) for base T, percentile 10, 440 total probabilites
> using optimized workers for all-context
> 0 B written to output: ak.bed [0 B/s]
[00:00:00] ----------------------------------------       0/4215502 genome positions 0/s 0s
> 0 rows written
> 0 ~records errored                                                                                                                                                                      
thread '<unnamed>' panicked at modkit-core/src/pileup/pileup_processor.rs:1894:9:
STRANDS 2
stack backtrace:
   0:     0x5562c5548b32 - <std::sys::backtrace::BacktraceLock::print::DisplayBacktrace as core::fmt::Display>::fmt::hdcfcb6d4c8489523
   1:     0x5562c556cbb3 - core::fmt::write::h8a494366950f23bb
   2:     0x5562c5545953 - std::io::Write::write_fmt::h6556609fca33d0b1
   3:     0x5562c5548982 - std::sys::backtrace::BacktraceLock::print::hb2a626a81e06b2dc
   4:     0x5562c5549e93 - std::panicking::default_hook::{{closure}}::h4f78485264f12d10
   5:     0x5562c5549c75 - std::panicking::default_hook::h2c66fc99e962531d
   6:     0x5562c554a825 - std::panicking::rust_panic_with_hook::h33ac55f64bbd807d
   7:     0x5562c554a5da - std::panicking::begin_panic_handler::{{closure}}::h30e7cb89678a57fe
   8:     0x5562c5549039 - std::sys::backtrace::__rust_end_short_backtrace::hed60f27456c16ced
   9:     0x5562c554a26d - __rustc[de2ca18b4c54d5b8]::rust_begin_unwind
  10:     0x5562c4602eb0 - core::panicking::panic_fmt::h62f63d096dd276af
  11:     0x5562c4b08682 - mod_kit::pileup::pileup_processor::ACountsMatrix::incr_call::hb767530f0ff78e6b
  12:     0x5562c4afee89 - <mod_kit::pileup::pileup_processor::DnaPileupWorker<T,M,_,_> as mod_kit::pileup::pileup_processor::PileupWorker>::process::hbead0a384f097cd2
  13:     0x5562c49ade03 - std::sys::backtrace::__rust_begin_short_backtrace::hb509df9362b458d6
  14:     0x5562c49bef98 - core::ops::function::FnOnce::call_once{{vtable.shim}}::hbefd88c46c41be61
  15:     0x5562c554db9f - std::sys::pal::unix::thread::Thread::new::thread_start::hb6e99e73da4d28f8
  16:     0x14cafdb20b7b - start_thread
                               at ./nptl/pthread_create.c:448:8
> Done, processed 0 rows.
> 0 B written to output: ak.bed [0 B/s]                                                                                                                                                   
thread 'main' panicked at modkit-core/src/pileup/subcommand.rs:1659:34:
worker thread 1 paniced: Any { .. }
stack backtrace:
   0:     0x5562c5548b32 - <std::sys::backtrace::BacktraceLock::print::DisplayBacktrace as core::fmt::Display>::fmt::hdcfcb6d4c8489523
   1:     0x5562c556cbb3 - core::fmt::write::h8a494366950f23bb
   2:     0x5562c5545953 - std::io::Write::write_fmt::h6556609fca33d0b1
   3:     0x5562c5548982 - std::sys::backtrace::BacktraceLock::print::hb2a626a81e06b2dc
   4:     0x5562c5549e93 - std::panicking::default_hook::{{closure}}::h4f78485264f12d10
   5:     0x5562c5549c75 - std::panicking::default_hook::h2c66fc99e962531d
   6:     0x5562c554a825 - std::panicking::rust_panic_with_hook::h33ac55f64bbd807d
   7:     0x5562c554a5da - std::panicking::begin_panic_handler::{{closure}}::h30e7cb89678a57fe
   8:     0x5562c5549039 - std::sys::backtrace::__rust_end_short_backtrace::hed60f27456c16ced
   9:     0x5562c554a26d - __rustc[de2ca18b4c54d5b8]::rust_begin_unwind
  10:     0x5562c4602eb0 - core::panicking::panic_fmt::h62f63d096dd276af
  11:     0x5562c46032d6 - core::result::unwrap_failed::h95bc3f5a607b2c95
  12:     0x5562c49cc335 - core::result::Result<T,E>::expect::hb714fbacf74e1ed4
  13:     0x5562c49eebb7 - mod_kit::pileup::subcommand::ModBamPileup::run::hf7ff9eb12f1f321d
  14:     0x5562c460fc1a - modkit::main::h132b333ca12c63c2
  15:     0x5562c4615766 - std::sys::backtrace::__rust_begin_short_backtrace::h9bb5aca9ca64d7de
  16:     0x5562c4613d22 - std::rt::lang_start::{{closure}}::h5a4f1a2efbc5fbb6
  17:     0x5562c553b160 - std::rt::lang_start_internal::h34f9328d113fd60a
  18:     0x5562c4610295 - main
  19:     0x14cafdab7ca8 - __libc_start_call_main
                               at ./csu/../sysdeps/nptl/libc_start_call_main.h:58:16
  20:     0x14cafdab7d65 - __libc_start_main_impl
                               at ./csu/../csu/libc-start.c:360:3
  21:     0x5562c46037a9 - _start
  22:                0x0 - <unknown>