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

[SchwarzMarek]

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>

[ArtRand]

@SchwarzMarek Thanks for the report, I'll take a look into fixing this.

Could you send me the BAM (with the header) and the reference, or just the CP171645 contig? I'll send you an email. Based on the result from check-tags there is one mod on each A and T primary sequence bases, so --combine-mods wouldn't do anything - it still shouldn't crash like this, however.

Thanks.

[ArtRand]

@SchwarzMarek Looks like my email to you got bounced back. Could you try me: art.rand [at] nanoporetech [dot] com

email used to be easy...