Keep base modification tags in step with trimming (--update-mods)

[wdecoster]

Closes #69.

--update-mods

MM encodes each modified base as a count of canonical bases to skip, so trimming or splitting a read silently invalidates it — the tags end up describing a sequence the read no longer has, and htslib rejects them outright (MM tag refers to bases beyond sequence length), losing the modification data entirely.

The new opt-in --update-mods recomputes MM, ML and MN for whatever part of the read is kept. It works with all four trimming approaches, since they all express their result as a (start, end) range over the original record. Calls outside the kept range are dropped along with their ML probabilities; lists that lose every call are kept as empty lists so that "called here, nothing found" survives.

Tags that cannot be rewritten are fatal rather than a warning — writing them through would attach modification coordinates to a sequence they do not describe, which is the bug the flag exists to prevent. That covers malformed MM/ML, MM and ML disagreeing on the number of calls, MM skipping past the end of the sequence, a stale MN, and ML without MM. Rerunning without the flag writes them through unchanged.

Only MM, ML and MN are corrected; qs, ns, ts and du are passed through and documented as going stale.

Header bug fix

The _segment_N suffix added by split-by-low-quality was appended to the end of the header line rather than to the read name, so with tab-separated SAM tags it landed on the last tag:

@read1  MM:Z:C+m?,1,1;  ML:B:C,200,10   qs:f:20.1_segment_1

samtools import parsed qs:f:20.1_segment_1 back as qs:f:20.1 and dropped the suffix, giving every segment of a read the same name. The header is now reassembled from id and description before being split on tabs, so tags are found wherever the parser's first-space split happened to fall (a tag value containing a space would otherwise hide every later tag), and the suffix goes on the name.

write_all fix

Write::write may consume only part of the buffer, and its result was discarded with let _ =. Records at or above the BufWriter's 8 KiB capacity bypass the buffer and reach stdout in a single write call, so any read over about 4 kb — routine here — could be silently truncated. The same discarded Result also swallowed write errors on that path, so a read that failed to write produced no message and no non-zero exit.

Also fixes chopper ... | head, which previously died with a Rust panic and exit 101 on the broken pipe.

Validation

Cross-checked against htslib (via pysam) rather than only against itself:

• 1000 reads under fixed-crop: all 35,050 surviving calls decode at the expected position with the expected probability.
• Same reads under split-by-low-quality: 37,404 calls across 8,144 segments, zero mismatches.
• Output round-trips through samtools import -T '*' with no htslib warnings.
• Untrimmed input is a byte-exact identity transform; untagged input is byte-identical to the previous binary.
• N+ lists are counted over every base, matching htslib's freq[15] = l_qseq.

54 unit tests, clippy and fmt clean. Cost on 100k reads / 362 MB: 0.27 s without the flag, 0.57 s with; zero when unset.

🤖 Generated with Claude Code

https://claude.ai/code/session_01B6KA7ghXr29T6ghPxdM9v9