unwrap error on bam_sa_parser

[mrvollger]

Hi,

I was testing out sawfish on a small bam and am running into this error:

sawfish discover --ref test.fa  --clobber --threads 16 --bam test.bam --output-dir temp/test/sawfish
[2024-12-11][15:21:29][sawfish][INFO] Starting sawfish 0.12.7
[2024-12-11][15:21:29][sawfish][INFO] cmdline: sawfish discover --ref test.fa --clobber --threads 16 --bam test.bam --output-dir temp/test/sawfish
[2024-12-11][15:21:29][sawfish][INFO] Running on 16 threads
[2024-12-11][15:21:29][sawfish][INFO] Writing discover settings to file: 'temp/test/sawfish/discover.settings.json'
[2024-12-11][15:21:29][sawfish][INFO] Reading reference genome from file '/mmfs1/gscratch/stergachislab/mvollger/projects/dev-k-mer-variant-phasing/k-mer-test-data/test.fa'
[2024-12-11][15:21:29][sawfish][INFO] Processing alignment file '/mmfs1/gscratch/stergachislab/mvollger/projects/dev-k-mer-variant-phasing/k-mer-test-data/test.bam'
[00:00:00] [==============================>         ] Processed alignments on 48,333 of 64,444 ref genome kb (75%)
thread '<unnamed>' panicked at src/bam_sa_parser.rs:156:77:
called `Option::unwrap()` on a `None` value
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
Aborted (core dumped)

Uploaded the failing test data here:
https://s3-us-west-2.amazonaws.com/stergachis-public1/index.html?prefix=Mitchell/k-mer-test-data/

Thanks!
Mitchell

[ctsa]

Thanks for sharing the test data for this one Mitchell,

In this case the split read SA tag includes split mappings to "chr4", but it looks like this has been removed from the reduced "test.fa" reference sequence.

I'll followup tomorrow with a proper error message for this case.

[ctsa]

Fix will have this behavior:

thread '<unnamed>' panicked at src/bam_sa_parser.rs:160:17:
In read 'm54329U_210328_013807/592310/ccs', the SA aux tag desribes a split read mapped to chr4:49100914 (in segment 0), which is not found in the input reference fasta

[mrvollger]

Great, thanks Chris!

I think this is a small bug in pbmm2 then, not dropping SA tags when remapping. My process here was:
(pbmm2 aligned bam) -> (subset bam to small region) -> (realign bam to a subsetted reference without any chr4)
Does that sound like a small pbmm2 bug to you since it accepts aligned bam as input?

[ctsa]

Agreed it sounds like a stale aux tag that the mapper would ideally have reset, though I'm less familiar with the support level for mapped input. For a quick workaround you could use samtools reset prior to re-mapping. I can check into this a bit more on pbmm2 side.

[armintoepfer]

We never had that use case of a reduced contig set. It's fixed now and will be part of the next pbmm2 release

[mrvollger]

Awesome, thanks @armintoepfer and @ctsa!

[ctsa]

Included in https://github.com/PacificBiosciences/sawfish/releases/tag/v0.12.8, closing as completed.

[armintoepfer]

We've released a new pbmm2 version on bioconda.