Errors on particular BAM file

[cmdcolin]

I saw this on the full bam file but selecting a small subset appears to reproduce

❯❯❯ samtools view https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data_somatic/HG008/Liss_lab/analysis/NIH_HiFi_mm2_mapping/HG008-T_PacBio-HiFi-Revio_20240125_116x_GRCh38-GIABv3_mm2.bam chr1:1-100000 -o subset.bam
❯❯❯ samtools index subset.bam
❯❯❯ RUST_BACKTRACE=1 ./perbase-linux-amd64 base-depth subset.bam
[2025-09-16T18:43:07Z INFO  perbase::commands::base_depth] Running base-depth on: "subset.bam"
[2025-09-16T18:43:07Z INFO  perbase_lib::par_granges] Using 14 worker threads.
[2025-09-16T18:43:07Z INFO  perbase_lib::par_granges] Creating channel of length 21681324 (* 120 bytes to get mem)
[2025-09-16T18:43:08Z INFO  perbase_lib::par_granges] Reading from "subset.bam"
[2025-09-16T18:43:08Z INFO  perbase_lib::par_granges] Processing TID 0:0-248956422

thread '<unnamed>' panicked at /home/runner/.cargo/registry/src/index.crates.io-1949cf8c6b5b557f/rust-htslib-0.39.5/src/bam/record.rs:1540:6:
index out of bounds: the len is 0 but the index is 148
stack backtrace:
   0: __rustc::rust_begin_unwind
   1: core::panicking::panic_fmt
   2: core::panicking::panic_bounds_check
   3: <core::iter::adapters::flatten::FlatMap<I,U,F> as core::iter::traits::iterator::Iterator>::next
   4: perbase_lib::par_granges::ParGranges<R>::process::{{closure}}::{{closure}}::{{closure}}::{{closure}}
   5: rayon::iter::plumbing::bridge_producer_consumer::helper
   6: rayon_core::join::join_context::{{closure}}
   7: rayon::iter::plumbing::bridge_producer_consumer::helper
   8: rayon_core::join::join_context::{{closure}}
   9: rayon::iter::plumbing::bridge_producer_consumer::helper
  10: rayon_core::join::join_context::{{closure}}
  11: rayon::iter::plumbing::bridge_producer_consumer::helper
  12: rayon_core::join::join_context::{{closure}}
  13: rayon_core::thread_pool::ThreadPool::install::{{closure}}
  14: <rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute
  15: rayon_core::registry::WorkerThread::wait_until_cold
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.
❯❯❯ ./perbase-linux-amd64 --version
perbase 1.1.0

here is backtrace full

 RUST_BACKTRACE=full ./perbase-linux-amd64 base-depth subset.bam
[2025-09-16T18:44:16Z INFO  perbase::commands::base_depth] Running base-depth on: "subset.bam"
[2025-09-16T18:44:16Z INFO  perbase_lib::par_granges] Using 14 worker threads.
[2025-09-16T18:44:16Z INFO  perbase_lib::par_granges] Creating channel of length 21681324 (* 120 bytes to get mem)
[2025-09-16T18:44:17Z INFO  perbase_lib::par_granges] Reading from "subset.bam"
[2025-09-16T18:44:17Z INFO  perbase_lib::par_granges] Processing TID 0:0-248956422

thread '<unnamed>' panicked at /home/runner/.cargo/registry/src/index.crates.io-1949cf8c6b5b557f/rust-htslib-0.39.5/src/bam/record.rs:1540:6:
index out of bounds: the len is 0 but the index is 148
stack backtrace:
   0:     0x555d93ee4c26 - <std::sys::backtrace::BacktraceLock::print::DisplayBacktrace as core::fmt::Display>::fmt::h46a716bba2450163
   1:     0x555d93e2fbb3 - core::fmt::write::h275e5980d7008551
   2:     0x555d93ee46ef - std::io::Write::write_fmt::h561a66a0340b6995
   3:     0x555d93ee4963 - std::sys::backtrace::BacktraceLock::print::hafb9d5969adc39a0
   4:     0x555d93ee43ce - std::panicking::rust_panic_with_hook::h409da73ddef13937
   5:     0x555d93f175d8 - std::panicking::begin_panic_handler::{{closure}}::h159b61b27f96a9c2
   6:     0x555d93f17539 - std::sys::backtrace::__rust_end_short_backtrace::h5b56844d75e766fc
   7:     0x555d93f17b8c - __rustc[4794b31dd7191200]::rust_begin_unwind
   8:     0x555d93d357ff - core::panicking::panic_fmt::hc8737e8cca20a7c8
   9:     0x555d93d357d9 - core::panicking::panic_bounds_check::hda0827d94e974e71
  10:     0x555d93d552d8 - <core::iter::adapters::flatten::FlatMap<I,U,F> as core::iter::traits::iterator::Iterator>::next::h1b63daffda1becfc
  11:     0x555d93d56115 - perbase_lib::par_granges::ParGranges<R>::process::{{closure}}::{{closure}}::{{closure}}::{{closure}}::h2f9776ee368b35b5
  12:     0x555d93d8f733 - rayon::iter::plumbing::bridge_producer_consumer::helper::he366852abe0f928e
  13:     0x555d93d4c8ef - rayon_core::join::join_context::{{closure}}::h173a48138f4b378e
  14:     0x555d93d8fb5b - rayon::iter::plumbing::bridge_producer_consumer::helper::he366852abe0f928e
  15:     0x555d93d4c8ef - rayon_core::join::join_context::{{closure}}::h173a48138f4b378e
  16:     0x555d93d8fb5b - rayon::iter::plumbing::bridge_producer_consumer::helper::he366852abe0f928e
  17:     0x555d93d4c8ef - rayon_core::join::join_context::{{closure}}::h173a48138f4b378e
  18:     0x555d93d8fb5b - rayon::iter::plumbing::bridge_producer_consumer::helper::he366852abe0f928e
  19:     0x555d93d4e398 - rayon_core::join::join_context::{{closure}}::h4b3a052ea4f34d43
  20:     0x555d93d47b61 - rayon_core::thread_pool::ThreadPool::install::{{closure}}::h2e89851a859aea23
  21:     0x555d93d99134 - <rayon_core::job::StackJob<L,F,R> as rayon_core::job::Job>::execute::h8ae08bd33531c99d
  22:     0x555d93d37f7a - rayon_core::registry::WorkerThread::wait_until_cold::he08018336a945cd2
  23:     0x555d93e50345 - std::sys::backtrace::__rust_begin_short_backtrace::h8b80dd194ccddbdf
  24:     0x555d93e52987 - core::ops::function::FnOnce::call_once{{vtable.shim}}::h9309d2923e070db7
  25:     0x555d93f17e8b - std::sys::pal::unix::thread::Thread::new::thread_start::h1822d22fde68314f
  26:     0x7941ac4a27f1 - start_thread
                               at ./nptl/pthread_create.c:448:8
  27:     0x7941ac533c9c - __GI___clone3
                               at ./misc/../sysdeps/unix/sysv/linux/x86_64/clone3.S:78:0
  28:                0x0 - <unknown>

[cmdcolin]

looks like the error message refers to htslib...not sure if that means it is upstream bug

[ghuls]

It happens because you don't have a sequence in your BAM file (but *):

m84039_240113_032943_s4/114821634/ccs   256     chr1    10001   0       297S15=1I23=1X5=1X5=1X5=1X16=1D34=2I6=1I18=2I12=1I2=1X10=2I20=1I19=4D33=2D1X4=1X6=1D1=1X5=2X11=1X15=1X48=111D5=1X18=1X23=1X25=1X12=2X58=1X8=1X1=1X18=2X2=8I1=1X173=7D1X6=1X2=3I1X2=1X1=1X5=162I1X54=1X5=47I7=1X29=1X40=1D16=1I17=1X32=1X31=8I1=1X6=1X37=1X20=1X19=1X12=1X60=1D67=1X112=1X47=1X50=1X20=1X12=1X22=1X60=1I51=1X14=1X15=1D21=1X65=1I278=1X95=1D43=1X13=1X63=1X63=1X3=4D15=1X53=1X46=1X8=1X16=1X26=1X159=1X43=1X1=1X14=1X17=1X1=1X4=1X19=1X7=1X30=1X16=1X12=1X19=1X21=1X28=1X8=1X37=1X75=1X16=1X36=1X1=1X2=1X16=1X19=1X1=1X141=1X24=1X72=1X23=1X83=1X158=1X131=1X38=1X73=1X33=1X25=1X7=1X217=1I51=1X8=12D79=1X28=1X120=1I105=1X5=1X67=1D20=1X11=1X4=1X6=1X3=1X209=1X37=1X118=1X50=1D4=1X53=1X69=1X150=1X97=1X20=2I95=1X71=1X65=1X42=1X83=1I18=1X145=1X34=1D60=1D78=1X35=1X18=1X29=1I51=1I75=1X164=1X69=1X117=1X1=1X40=1X5=1X117=1X197=1X306=1X13=2I202=1X32=1X17=1X58=1X89=5I3=1X5=1X63=1X58=1X2=1X63=1X44=1I25=1X29=1X5=1X221=1X76=1X131=1X114=1X21=1X68=1D11=1X23=1X7=1X73=1X48=1X86=1D80=1X53=1X22=1X4=1X24=1X75=1X35=1X114=1X1=1X91=1X3=1X55=1I107=1D8=1X81=1X20=1D72=1X84=1X91=1X105=1X92=1I535=1I53=1X107=1X63=1I5=1X27=1D114=1D45=1X69=1X58=1X68=1X22=1X58=1X17=1X106=1D27=1X47=1X11=1D97=2X46=1X12=1X42=1D112=1I23=1X34=1X9=1X59=1D1=1X84=1X88=1D72=1X259=1I98=1X3=1D165=1X75=1X271=1X15=1X35=1X2=1X107=1X31=1X20=1X24=1X62=1X5=1X46=1X9=1X136=1X8=1X118=1X156=1X87=1X65=1X21=1X26=1I16=1I24=1D49=1X88=1X54=1D69=1X192=1X49=1X3=1X25=1X26=1X4=1X170=1X197=1X33=1X27=1I109=1I15=1X4=2X27=1X5=38D34=1I78=1D43=1X40=1X101=1X57=1X12=1X13=1X22=1X29=2X23=1X69=1D51=1X63=1X5=1X34=1X17=1I20=4D96=1I41=3D95=1D57=1X109=1X82=1X162=1X38=1X2=1X54=1X2=1X19=1D78=1I24=1X68=1D145=1X46=1X16=1X29=1X333=1X19=1X32=1X92=1X41=1I126=1X136=1D15=1X47=1I38=1X43=1X67=1X12=1X61=1X25=5I27=1D46=1X5=1I23=2X105=1D35=1X69=1X53=1X193=1X1=1X34=1X54=1I13=1D134=1X3=1X1=1X43=1X97=1I30=1X57=1D29=1X36=1X28=1X31=1X32=1X9=1D93=1X15=1X28=1X18=1X52=1D16=1X8=1X51=1X19=1I38=1X93=9I74=1X46=1X13=1X34=1X12=1X64=1X95=1X50=1X18=1D16=1X49=1X5=1X50=1X80=1X32=3I1X1=1X22=6D8=2D51=1X43=1X237=7I13=1X3=1I11=1X158=1X12=1X222=1X29=2D142=1X114=1D24=1X62=1X5=1X31=1X9=1I232=1X37=1X9=1X61=1D154=1D61=1X81=1X61=1X165=1X4=1X160=1X116=1I37=1D15=1D36=1X53=1X20=1X13=1X3=1D57=1X9=1X2=1X12=5I89=1D32=1X10=1D194=1I72=1X370=1X150=1X10=1X130=1X156=1X7=1X255=1I60=1X92=1X33=1X1=1X23=1X54=1X28=1X99=1X93=1X247=1X150=1X6=1D353=1D83=1X101=1D29=1X37=1X5=1X6=1D1X186=1X63=1X35=1X24=2X13=1X25=1D21=1D9=1X57=1X17=1X24=1X30=1X2=1X46=1I38=1X5=1D91=2D83=1X5=1X176=1X79=1X93=1I34=1X119=1X5=1D68=1D33=1D14=1X20=1X15=1D54=1X34=1I144=1X154=1X10=1X4=1X202=1X135=1X43=1X102=1X17=1X107=1X164=1X4=1X61=1X240=1I114=1X1=1X26=1X9=1X20=1X12=2X45=1X3=1X7=1X291=1I67=1X16=1X38=1X94=1X6=1X137=1X39=1X114=1X5=1X6=  *       0       0       *       *       NM:i:990        ms:i:27062      AS:i:26567      nn:i:0  tp:A:S  cm:i:2083       s1:i:22419      de:f:0.0181     rl:i:1638       ML:B:C,114,56,150,190,156,161,126,78,190,65,220,123,37,30,92,155,56,175,145,134,192,218,204,100,59,61,18,50,108,136,175,57,54,94,130,147,184,158,90,104,89,55,120,135,145,228,234,18,135,228,27,126,156,217,197,158,220,248,59,187,71,53,81,208,198,142,176,78,73,60,113,191,76,15,68,24,72,156,49,80,30,59,58,13,53,80,167,12,102,60,35,57,158,78,87,37,44,63,112,134,143,250,238,237,170,250,243,147,84,101,23,25,28,35,102,163,186,139,214,145,219,228,110,170,165,84,170,221,136,226,184,95,186,150,144,128,76,175,139,193,235,250,141,112,191,168,178,100,78,178,33,17,29,117,231,254,211,219,208,153,145,150,211,251,235,58,189,74,47,106,135,117,74,58,177,83,109,242,172,81,46,70,189,149,165,113,135,140,236,241,153,158,55,110,110,70,163,160,149,188,225,251,246,208,91,190,53,169,139,94,121,151,186,82,204,199,25,83,243,219,157,49,96,130,79,76,151,164,155,127,71,74,131,134,156,202,221,251,252,115,132,206,190,200,217,164,138,45,56,183,101,173,143,176,165,199,143,59,212,239,199,130,87,139,230,195,233,143,143,59,213,148,180,96,248,226,185,132,117,54,113,178,130,150,216,150,111,63,197,113,170,242,70,163,126,117,176,183,46,72,130,128,194,225,121,167,248,252,149,195,215,192,173,216,129,105,183,146,89,222,87,209,46,98,158,135,202,126,45,67,203,167,203,88,129,178,155,247,140,250,248,165,68,122,232,163,181,208,77,185,173,100,202,171,192,182,145,41,232,73,83,239,164,159,124,73,155,99,72,162,139,103,95,136,225,208,68,154,137,75,176,56,64,174,138,178,106,180,61,150,57,77,122,173,251,248,88,130,137,103,143,197,96,123,235,131,91,178,141,174,157,206,124,123,130,52,137,166,169,199,75,200,248,115,185,229,126,137,55,114,146,171,77,152,211,236,143,142,227,177,106,158,57,167,218,188,148,161,172,119,120,91,163,105,126,222,76,194,57,165,243,57,127,104,141,131,136,177,111,112,151,97,153,69,165,193,85,157,164,20,117,135,116,145,181,52,52,6,39,41,28,29,41,55,124,142,78,58,148,97,206,53,114,28,122,119,31,116,46,56,40,126,168,102,40,58,77,141,54,48,158,40,53,59,69,37,172,46,65,65,48,142,50,36,29,62,166,26,43,28,42,34,102,126,27,80,19,64,138,90,93,37,25,88,118,37,60,27,78,48,45,71,76,144,37,20,154,112,141,89,106,34,52,52,49,75,165,84,76,106,126,37,72,117,66,131,62,73,52,139,96,49,249,191,197,43,201,145,243,44,205,56,161,183,163,46,221,204,230,188,47,212,235,116,233,129,155,121,121,47,158,133,209,142,128,215,155,100,149,163,231,103,200,111,247,103,123,205,142,63,213,253,41,243,146,146,162,218,137,124,174,229,98,90,171,66,189,162,98,148,163,66,41,47,94,110,161,158,69,185,60,187,47,127,176,40,26,110,42,80,26,152,41,154,101,64,92,48,49,183,40,52,123,73,41,201,87,74       MM:Z:C+m?,349,5,2,4,7,2,4,3,3,3,0,3,2,3,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,1,0,1,0,1,0,1,4,0,2,0,1,4,0,3,0,1,1,0,1,0,1,0,1,4,0,1,0,0,1,1,1,2,0,0,2,1,0,1,1,0,0,1,1,0,3,1,0,0,0,0,1,4,0,1,0,8,5,3,3,5,10,7,1,10,1,14,2,8,8,3,51,1,2,6,4,5,2,19,5,68,1,3,13,4,4,23,2,18,3,1,6,31,0,1,9,0,8,22,30,74,19,14,20,9,0,27,1,9,1,24,8,19,1,21,0,84,41,52,16,7,3,13,13,11,20,19,5,0,3,0,5,31,6,3,1,17,8,5,7,1,34,8,2,12,9,13,12,2,2,6,27,3,28,20,3,8,7,18,5,5,6,9,0,12,22,3,0,2,14,0,99,3,13,4,4,6,0,9,14,21,11,3,2,48,74,7,17,11,32,4,8,11,17,5,3,41,15,9,1,2,5,0,7,2,28,7,4,40,55,13,0,47,8,6,4,6,13,18,3,13,2,16,18,3,7,1,8,40,1,2,2,0,5,13,12,16,26,17,7,27,5,2,4,22,15,8,10,9,27,19,18,15,10,17,58,4,5,18,16,11,11,12,9,28,26,0,40,9,22,2,23,8,2,22,2,0,8,51,5,23,34,21,26,9,6,12,0,7,22,27,11,42,8,10,22,3,18,32,33,14,8,3,7,22,27,14,10,60,20,25,10,18,6,10,5,17,68,4,5,1,4,16,1,7,3,11,24,53,16,79,34,14,3,25,124,11,19,32,1,1,2,7,5,0,15,5,2,4,3,0,23,1,67,19,19,76,3,28,9,1,19,23,3,6,9,4,0,1,1,8,7,7,3,4,3,0,5,3,0,3,1,0,0,6,1,0,0,1,1,3,3,3,13,2,1,0,5,0,3,4,0,0,3,0,0,0,2,2,1,1,0,5,3,1,2,4,2,12,4,3,7,6,0,3,0,7,1,5,3,1,1,2,3,1,6,5,1,2,0,3,2,1,1,1,4,1,6,0,0,6,0,2,0,3,1,0,0,3,2,5,3,0,2,5,2,0,2,6,14,5,1,13,4,4,4,3,3,14,21,34,10,3,3,3,0,11,20,5,13,39,10,10,26,70,16,77,1,8,0,28,3,11,16,63,21,49,16,51,15,49,11,9,50,1,48,44,55,7,14,23,9,98,4,1,5,13,41,2,9,3,21,5,11,12,45,25,9,15,58,3,15,3,19,44,6,2,56,32,4,5,1,1,4,3,11,46,25,41,10,36,4,29,14,58,1,6,8,77,40,1,22,25,25,0,47,7,15,19,15,11,13,13,8,17,52,35,2,7,43;

causing out of bound access in : https://github.com/rust-bio/rust-htslib/blob/229d5a265685fd25e59e807524a5df27d4057c83/src/bam/record.rs#L1806-L1809

[cmdcolin]

interesting. it is not necessarily a rare thing for BAM files to have * in the SEQ field though. Secondary reads (e.g. flag 256) frequently do not have SEQ due to things like lh3/minimap2#458

[sstadick]

@ghuls thanks for commenting on this!
@cmdcolin sorry I missed it!

TIL * is a valid SEQ character per the spec in the SEQ field. Could you try running with only-depth? base-depth needs to know what base is at each position, and if the sequence isn't there that won't work.

I guess that only way to guard against that is to check the first byte of every seq to see if it's a *. I'll work on a fix for this in the next few days.

[ghuls]

@sstadick In the BAM file, you will not see *, but you can check record.seq().is_empty() first to see if you have a stored sequence or not. See my pull request.

For accessing the quality per base, you need a similar check for the sequence length.

I still have the feeling that there is a bug in rust_htslib for cases where you have a sequence, but no qual set and you try to access the qual at a position:
https://github.com/rust-bio/rust-htslib/blob/master/src/bam/record.rs#L634-L640

10. SEQ: segment SEQuence. This field can be a ‘*’ when the sequence is not stored. If not a ‘*’, the length
of the sequence must equal the sum of lengths of M/I/S/=/X operations in CIGAR. An ‘=’ denotes the
base is identical to the reference base. No assumptions can be made on the letter cases.
11. QUAL: ASCII of base QUALity plus 33 (same as the quality string in the Sanger FASTQ format). A
base quality is the phred-scaled base error probability which equals −10 log 10 Pr{base is wrong}. This
field can be a ‘*’ when quality is not stored. If not a ‘*’, SEQ must not be a ‘*’ and the length of the
quality string ought to equal the length of SEQ.

It would be nice if some/all perbase dependencies are updated:

gzp = "2.0"

rust-htslib = { version = "0.51", default-features = false, features = ["libdeflate"] }

[sstadick]

Release 1.2 should be out sometime today. Also including #93

Thanks @ghuls for the fix and @cmdcolin for the report!