Skip to content

Support reading and writing unaligned BAM #79

Description

@wdecoster

Motivation

chopper currently reads and writes FASTQ only. For modified-basecalled ONT data that means routing everything through samtools fastq -T MM,ML,MN | chopper | samtools import -T '*', which smuggles SAM tags through a FASTQ header — a format that was never designed to carry them.

Two bugs found while investigating #69 are symptoms of exactly that:

  • The _segment_N suffix added by --trim-approach split-by-low-quality was appended to the end of the whole header line, i.e. onto the last SAM tag, producing ...\tqs:f:20.1_segment_1. samtools import silently parsed that back as qs:f:20.1 and dropped the suffix, so both segments came back with the same read name.
  • MM/ML were passed through as opaque text, so after any trimming they described a sequence the read no longer had. htslib rejects these outright (MM tag refers to bases beyond sequence length), so the modification data was simply lost.

Both are now fixed for the FASTQ path (--update-mods), but the underlying mismatch stays: in a uBAM these are typed fields that cannot be passed through unexamined, and the other tags that trimming invalidates (qs, ns, ts, du) become addressable instead of silently going stale. The @RG/@PG header that samtools fastq discards would also survive. Dorado emits uBAM natively, so chopper sitting in the middle of a BAM pipeline is the natural shape now, and it saves two full encode/decode passes over the data.

Proposed approach

Use noodles (noodles-bam + noodles-sam), not rust-htslib.

This is the deciding constraint: chopper ships a statically linked musl binary (make musl), and htslib in a static musl build is a known pain. noodles is pure Rust and sidesteps it entirely. Note that rust-htslib is not currently in the dependency tree — it is an optional feature of the minimap2 crate that we do not enable. There are already C dependencies (minimap2-sys, libz-sys, bzip2, xz2), but adding htslib on top of the musl target is a different order of problem.

Work involved

  • Abstract the record type so filters and trimmers operate on (seq, qual) slices. This should be mostly mechanical: the trimmers only ever touch record.seq() and record.qual(), and every strategy already returns (start, end) offsets into the original record.
  • Detect BAM input (magic bytes) and add a reader.
  • Add a writer, propagating the input header and adding a @PG line recording the chopper invocation — users will expect the provenance.
  • Reject aligned input. Trimming invalidates POS and CIGAR, so anything that is not unmapped has to be an error rather than silently corrupted output.
  • Reuse the subset() core in src/modtags.rs for MM/ML/MN. The position arithmetic carries over unchanged; only the parse/serialize layer differs, since in BAM these arrive as typed fields rather than as text.

Rough estimate: a few days of focused work. It roughly doubles chopper's surface area, hence a deliberate feature rather than a bolt-on.

Notes

  • This does not obsolete --update-mods. Tagged FASTQ stays common and plenty of pipelines will not be restructured, so the FASTQ path still needs to be correct.
  • On the BAM path, updating modification tags should be the default rather than a flag. Once they are typed fields you cannot pass them through as opaque text without deciding what to do, so the opt-in framing stops making sense. --update-mods is therefore a FASTQ-path concept, not a permanent part of the interface.

One useful invariant for testing, established while implementing #69: for FASTQ, records are always in original read orientation, so the canonical base in MM is counted left to right and the +/- strand character is metadata. In BAM that is no longer guaranteed — reverse-complement handling is driven by the record's reverse flag (see htslib sam_mods.c, which keys off BAM_FREVERSE and counts the complement base right to left). Unmapped dorado output has no reverse records, which is another reason to reject aligned input rather than try to handle it.

Refs #69

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    enhancementNew feature or request

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions