[wdecoster]

Closes #79. chopper now reads unaligned BAM (uBAM), the default output of ONT basecallers, and can write it. The format is detected automatically. This PR changes I/O only: filtering, trimming, thresholds and FASTQ output are untouched.

Design

Crate: noodles (noodles-bam, -sam, -bgzf). It is pure Rust and MIT-licensed. I didn't use rust-htslib because it would pull in a C build of htslib (clang/bindgen), which complicates the static musl release binaries.

One pipeline for both record types. Filters and trimmers keep taking &fastq::Record. A BAM record gets a FASTQ view and keeps its original record alongside it:

pub enum RawRecord { Fastq(fastq::Record), Bam(RecordBuf) }   // off the reader
pub struct InputRecord {
    pub fastq: fastq::Record,   // what filters/trimmers see (phred+33, tags as SAM text after a tab)
    pub bam: Option<RecordBuf>, // original, for BAM output
}
pub enum WritableRecord { Fastq(String), Bam(RecordBuf) }      // built in the output format

get_valid_segments and every TrimStrategy are unchanged. --update-mods works for BAM without changes because modtags already operates on the SAM text of MM/ML/MN.

Detection. BAM is BGZF, i.e. gzip, so it goes through the existing gzip decoder and is recognised by the decompressed BAM\1 magic. stdin is now peeked for gzip so that cat x.bam | chopper works. Side effect: gzipped FASTQ on stdin, which used to fail to parse, now works too.

Output format defaults to the input format. A new flag, --output-format fastq|bam, converts between them; it is needed for the two cross-format cases.

Tags and header on a round trip

Input → output	Header	Tags
uBAM → uBAM	Passed through (@HD/@RG/@PG/@CO…) plus a @PG ID:chopper chained via PP after the last program	All kept, with exact type and original order (the test asserts the records are identical). Only name (_segment_N), seq and qual change, plus MM/ML/MN with --update-mods, which are replaced in place. Flags kept.
uBAM → FASTQ	n/a	All tags written as SAM fields after the name, like samtools fastq -T '*'. Byte-identical to htslib's formatting for the fixture.
FASTQ → uBAM	Rebuilt: @HD VN:1.6 SO:unknown + @PG	Header SAM fields become typed tags (integers narrowed to the smallest type, like htslib). Free text after the name goes into CO:Z:. Flag 4.
FASTQ → FASTQ	n/a	Unchanged, byte for byte.

Nothing is dropped silently. The known limitations, all documented in the README:

• Integer widths aren't kept through FASTQ, because SAM text has a single i type.
• A FASTQ RG:Z: tag gets no @RG header line.
• Position-dependent tags other than MM/ML/MN (qs, ns, ts, du, mv) go stale, as they already did for FASTQ.

The following are fatal errors rather than lossy output:

• aligned records: CIGAR/position would be invalidated, and reverse-strand reads are stored reverse complemented
• BAM records without qualities
• FASTQ header fields after a tab that aren't valid SAM tags, duplicate tags, and free text alongside an existing CO tag

Bug caught along the way

With --update-mods, a split segment that keeps no modification calls gets an empty ML:B:C. The SAM spec allows this, but noodles' SAM parser rejects it. text_to_tags now builds empty arrays itself; there is a unit test for it.

Validation

• FASTQ unchanged: test-data/regression/v0.14.1.sha256 holds digests of the v0.14.1 release binary's output for 14 command lines: every trim approach, the GC/quality/contam filters, --inverse, --update-mods, a multi-threaded run (order-insensitive) and stdin input. fastq_output_is_identical_to_v0_14_1 checks them. I also compared a wider set by hand, at -t 1 and -t 4 and including gzip input and stderr. All identical.
• BAM against htslib (pysam), outside the test suite: for uBAM→uBAM passthrough, fixed-crop and split-by-low-quality with --update-mods (28 segments), I decoded the original and the output reads independently with htslib's modified_bases. Every surviving call lands at pos - start with the same probability, MN equals the segment length, every other tag is unchanged in value, type and order, and sequence/qualities equal the slice.
• Fixtures in test-data/ubam/ are synthetic: 6 seeded random reads with dorado-like tags and header. They are written by htslib through make_fixtures.py, so tests compare against an independent implementation.
• make ci, the 4 #[ignore]d smoke tests and clippy -D warnings are clean. I did not run make musl locally (no musl-gcc here); the new crates are pure Rust.

Performance (rough, release build, 4 threads, ~200 MB of reads)

FASTQ→FASTQ 0.09 s · BAM→BAM 0.78 s · BAM→FASTQ 0.59 s · FASTQ→BAM 0.71 s. BAM input is decompressed on one thread (the existing MultiGzDecoder). BAM output uses noodles' multithreaded BGZF writer with --threads workers. A multithreaded BGZF reader would be a natural follow-up if BAM input throughput matters.

Test output

test formats::tests::bam_record_becomes_fastq_with_all_tags_as_text ... ok
test formats::tests::aligned_bam_records_are_rejected ... ok
test formats::tests::invalid_tag_text_is_an_error ... ok
test modtags::tests::chebi_codes_are_a_single_modification ... ok
test modtags::tests::headcrop_shifts_the_first_skip_count ... ok
test modtags::tests::mismatched_mm_and_ml_counts_are_an_error ... ok
test formats::tests::empty_arrays_keep_their_place ... ok
test formats::tests::mods_only_keeps_just_the_modification_tags ... ok
test formats::tests::output_header_keeps_input_and_chains_a_program_line ... ok
test formats::tests::tags_round_trip_through_text ... ok
test modtags::tests::malformed_mm_is_an_error ... ok
test modtags::tests::ml_without_the_separating_comma_is_an_error ... ok
test modtags::tests::mm_without_ml_is_valid ... ok
test modtags::tests::malformed_ml_is_an_error ... ok
test formats::tests::bam_records_without_qualities_are_rejected ... ok
test modtags::tests::ml_without_mm_is_an_error ... ok
test modtags::tests::malformed_mn_is_reported_as_an_mn_error ... ok
test modtags::tests::n_lists_do_not_fail_on_reads_without_literal_n ... ok
test modtags::tests::rewrite_tags_accepts_the_old_lowercase_tag_names ... ok
test modtags::tests::rewrite_tags_passes_through_reads_without_mod_tags ... ok
test modtags::tests::rewrite_tags_empties_ml_when_no_call_survives ... ok
test modtags::tests::rewrite_tags_updates_mm_ml_and_mn_and_leaves_others_alone ... ok
test modtags::tests::segment_keeps_only_the_calls_it_contains ... ok
test modtags::tests::stale_mn_is_an_error ... ok
test modtags::tests::segment_without_calls_keeps_an_empty_list ... ok
test modtags::tests::tags_describing_a_longer_sequence_are_an_error ... ok
test modtags::tests::multiple_codes_keep_their_probabilities_interleaved ... ok
test records::tests::test_record_to_string_keeps_sam_tags_on_single_segment ... ok
test records::tests::test_record_to_string_finds_tags_after_a_tag_value_with_a_space ... ok
test modtags::tests::multiple_lists_are_rewritten_independently ... ok
test records::tests::test_record_to_string_multiple_segments ... ok
test records::tests::test_record_to_string_leaves_mod_tags_alone_by_default ... ok
test modtags::tests::tailcrop_drops_calls_and_their_probabilities ... ok
test modtags::tests::u_is_treated_as_t ... ok
test records::tests::test_record_to_string_multiple_segments_with_desc ... ok
test records::tests::test_record_to_string_one_segment ... ok
test records::tests::test_record_to_string_segment_suffix_goes_on_name_not_tags ... ok
test records::tests::test_record_to_string_reports_broken_mod_tags ... ok
test records::tests::test_record_to_string_segment_suffix_with_tags_and_desc ... ok
test modtags::tests::untrimmed_range_is_unchanged ... ok
test records::tests::test_completed_record_to_string ... ok
test modtags::tests::n_counts_every_base_not_just_literal_n ... ok
test records::tests::test_record_to_string_suffixes_the_name_not_the_description ... ok
test records::tests::test_record_to_string_updates_mod_tags_when_requested ... ok
test tests::phred_score_to_probability_test ... ok
test tests::test_ave_qual ... ok
test records::tests::test_record_to_string_with_no_valid_segment - should panic ... ok
test tests::test_filter ... ignored
test tests::test_filter_with_best_read_segment_approach ... ignored
test tests::test_filter_with_contam ... ignored
test tests::test_empty_record_returns_empty_vec ... ok
test tests::test_filter_with_trim_by_quality_approach ... ignored
test tests::test_multiple_segments_only_some_valid ... ok
test tests::test_valid_before_and_after_trimming ... ok
test tests::test_too_short_record_filtered_before_trimming ... ok
test trimmers::tests::split_by_low_quality_edge_cases_test ... ok
test trimmers::tests::highest_quality_read_test ... ok
test trimmers::tests::fixed_crop_strategy_test ... ok
test trimmers::tests::split_by_low_quality_multiple_segments_test ... ok
test trimmers::tests::split_by_low_quality_strategy_test ... ok
test trimmers::tests::split_by_low_quality_window_test ... ok
test tests::test_quals ... ok
test trimmers::tests::trim_by_quality_strategy_test ... ok
test tests::test_record_qual_len ... ok
test tests::test_no_contam ... ok
test tests::test_contam ... ok
test result: ok. 62 passed; 0 failed; 4 ignored; 0 measured; 0 filtered out; finished in 0.02s
test fastq_description_is_kept_as_a_comment_tag ... ok
test aligned_bam_is_rejected ... ok
test ubam_to_fastq_writes_all_tags_like_samtools_fastq ... ok
test ubam_to_ubam_without_trimming_keeps_every_record_unchanged ... ok
test fastq_to_ubam_types_the_header_tags ... ok
test ubam_on_stdin_is_detected ... ok
test ubam_output_without_reads_is_a_valid_empty_bam ... ok
test ubam_to_fastq_matches_fastq_to_fastq_when_trimming ... ok
test multithreaded_ubam_output_holds_the_same_records ... ok
test ubam_to_ubam_trimming_rewrites_only_sequence_and_mod_tags ... ok
test split_segments_are_named_like_fastq_output ... ok
test fastq_output_is_identical_to_v0_14_1 ... ok
test result: ok. 12 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 4.86s

🤖 Generated with Claude Code