From 7034a03163db02ffbd6b20389255f4022c7a238f Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Sun, 27 Sep 2026 04:37:04 -0700 Subject: [PATCH 1/3] fix(runall/align): accept bwa's mid-pair -K split and keep both halves' tags With mixed single-end/paired input, bwa's even-read-count -K chunk cut can land between a pair's two reads; with -p smart pairing bwa then aligns them as two unpaired reads and emits them back to back under the pair's name. The subprocess reader failed on that output. It now recognizes the shape (no record paired, exactly two primaries) and splits the pair's unmapped template to match, clearing each half's pairing bits so the zipper merge copies each read's own tags and QC-fail flag. The split is accepted only from the bwa and bwa-mem3 presets, which run `mem -p -K`; a free-form --aligner::command keeps the loud "multiple primaries" error, since the same shape from an aligner that never pairs would otherwise split every pair silently. Also cap a preset's --aligner::chunk-size at i32::MAX: bwa parses -K with atoi into an int. And add a subprocess align/merge guard (tests/align_subprocess_split_parity.rs) that runs the preset over a fixture of substitutions, indels, unmappable, chimeric, discordant, repeat and zero-length reads with RX tags, checks every read keeps its tags (including across a forced mid-pair split), and that the output is identical across --threads values; it runs where a bwa-mem3 binary is available and skips otherwise. --- src/lib/aligner.rs | 72 ++- src/lib/pipeline/steps/align/merge.rs | 68 ++ src/lib/pipeline/steps/align/mod.rs | 60 +- src/lib/pipeline/steps/align/subprocess.rs | 328 ++++++++-- tests/align_common/mod.rs | 718 +++++++++++++++++++++ tests/align_subprocess_split_parity.rs | 179 +++++ 6 files changed, 1375 insertions(+), 50 deletions(-) create mode 100644 tests/align_common/mod.rs create mode 100644 tests/align_subprocess_split_parity.rs diff --git a/src/lib/aligner.rs b/src/lib/aligner.rs index 14c8d14c7..3095757e0 100644 --- a/src/lib/aligner.rs +++ b/src/lib/aligner.rs @@ -707,7 +707,17 @@ impl Default for AlignerOptions { pub(crate) enum ResolvedBackend { /// `--aligner::preset` or `--aligner::command`: the shell command to spawn /// via [`AlignerProcess::spawn`] (already substituted and validated). - Subprocess { command: String }, + Subprocess { + command: String, + /// Whether the aligner's output may carry bwa's mid-pair split (a + /// pair aligned as two unpaired reads because a `-K` chunk cut fell + /// between them). `true` only for the presets, which run `mem -p -K` + /// and so produce exactly that shape; a free-form `--aligner::command` + /// gets the loud "multiple primaries" error instead, because the same + /// shape from an aligner that never pairs would otherwise split every + /// pair silently. + accept_mid_pair_split: bool, + }, } /// Result of [`AlignerOptions::resolve`] — a ready-to-construct aligner @@ -751,6 +761,10 @@ pub(crate) enum ResolvedAlignerMode { Command, } +/// The largest `--aligner::chunk-size` a preset accepts: `i32::MAX`, since bwa +/// and bwa-mem3 parse `-K` with `atoi` into an `int`. +const MAX_PRESET_CHUNK_SIZE: u64 = 2_147_483_647; + impl AlignerOptions { /// Validate the option combination and produce a [`ResolvedAligner`] /// ready for the align stage's `backend_for` to construct. @@ -766,6 +780,7 @@ impl AlignerOptions { /// /// # Errors /// + /// - `--aligner::chunk-size` is zero, or exceeds `i32::MAX` with a preset. /// - Neither `--aligner::preset` nor `--aligner::command` was set. /// - Both were set (mutual exclusion violation). /// - Command mode + a preset-only flag (`--aligner-bin` or @@ -789,6 +804,17 @@ impl AlignerOptions { aligner's -K batch size and the pipeline's in-flight unmapped budget" ); } + // Every preset hands the chunk size to bwa's `-K`, which bwa and + // bwa-mem3 parse with `atoi` into an `int`: a larger value overflows + // there instead of reaching the aligner as given. Command mode sets its + // own `-K` (this flag only sizes its in-flight budget), so it is exempt. + if self.preset.is_some() && self.chunk_size > MAX_PRESET_CHUNK_SIZE { + bail!( + "--aligner::chunk-size must be at most {MAX_PRESET_CHUNK_SIZE} (got {}) with \ + --aligner::preset: bwa parses -K as a 32-bit int", + self.chunk_size + ); + } match (self.preset, self.command) { (None, None) => bail!( "--start-from align requires one of `--aligner::preset` \ @@ -811,7 +837,10 @@ impl AlignerOptions { let command = preset.build_command(reference, threads, self.chunk_size, aligner_bin); Ok(ResolvedAligner { - backend: ResolvedBackend::Subprocess { command }, + // Both presets run `mem -p -K`, whose smart pairing splits + // a pair across a chunk cut (bwa and bwa-mem3 share + // `bseq_read`'s even-count cut and `bseq_classify`). + backend: ResolvedBackend::Subprocess { command, accept_mid_pair_split: true }, chunk_size: self.chunk_size, threads: Some(threads), mode: ResolvedAlignerMode::Preset(preset), @@ -846,7 +875,7 @@ impl AlignerOptions { // hardcode any thread count). let command = substitute_template(&template, reference, top_threads)?; Ok(ResolvedAligner { - backend: ResolvedBackend::Subprocess { command }, + backend: ResolvedBackend::Subprocess { command, accept_mid_pair_split: false }, chunk_size: self.chunk_size, threads: None, mode: ResolvedAlignerMode::Command, @@ -1149,7 +1178,8 @@ mod tests { chunk_size: DEFAULT_ALIGNER_CHUNK_SIZE, }; let resolved = opts.resolve(&ref_path, 4, None).unwrap(); - let ResolvedBackend::Subprocess { command } = resolved.backend; + let ResolvedBackend::Subprocess { command, accept_mid_pair_split } = resolved.backend; + assert!(!accept_mid_pair_split, "command mode must not accept bwa's mid-pair split"); assert!(command.contains(&ref_path.display().to_string())); assert!(command.contains("-t 4")); assert!(matches!(resolved.mode, ResolvedAlignerMode::Command)); @@ -1327,6 +1357,40 @@ mod tests { assert_eq!(result, "bwa-mem3 mem -t 8 /data/{threads}/genome.fa /dev/stdin"); } + /// `resolve` caps a preset's `--aligner::chunk-size` at `i32::MAX`, since bwa + /// and bwa-mem3 parse `-K` with `atoi` into an `int`; the check runs before + /// preset validation, so the reference needs no index. + #[rstest] + #[case::bwa_mem3(AlignerPreset::BwaMem3, 2_147_483_648)] + #[case::bwa(AlignerPreset::Bwa, 4_294_967_296)] + fn resolve_rejects_preset_chunk_size_past_i32( + #[case] preset: AlignerPreset, + #[case] chunk_size: u64, + ) { + let tmp = tempfile::tempdir().unwrap(); + let ref_path = tmp.path().join("ref.fa"); + std::fs::write(&ref_path, b">chr1\nACGT\n").unwrap(); + let opts = AlignerOptions { preset: Some(preset), chunk_size, ..AlignerOptions::default() }; + let err = opts.resolve(&ref_path, 4, None).unwrap_err().to_string(); + let needle = format!("--aligner::chunk-size must be at most 2147483647 (got {chunk_size})"); + assert!(err.contains(&needle), "got: {err}"); + } + + /// Command mode sets its own `-K`, so a chunk size past `i32::MAX` (which + /// only sizes its in-flight budget) is accepted. + #[test] + fn resolve_accepts_command_mode_chunk_size_past_i32() { + let tmp = tempfile::tempdir().unwrap(); + let ref_path = tmp.path().join("ref.fa"); + std::fs::write(&ref_path, b">chr1\nACGT\n").unwrap(); + let opts = AlignerOptions { + command: Some("bwa-mem3 mem -t {threads} {ref} /dev/stdin".to_string()), + chunk_size: 2_147_483_648, + ..AlignerOptions::default() + }; + assert!(opts.resolve(&ref_path, 4, None).is_ok()); + } + /// `resolve` rejects `--aligner::chunk-size 0` (would reach the aligner as /// `-K 0` and zero the in-flight budget) with a flag-attributed error. #[test] diff --git a/src/lib/pipeline/steps/align/merge.rs b/src/lib/pipeline/steps/align/merge.rs index 33029d60b..68e78942f 100644 --- a/src/lib/pipeline/steps/align/merge.rs +++ b/src/lib/pipeline/steps/align/merge.rs @@ -299,6 +299,74 @@ mod tests { ); } + /// A pair split across a mid-pair `-K` cut zips each half with the aligner's + /// unpaired record for that read. Both halves (and the second read's + /// supplementary) must receive their own unmapped read's tags, and the + /// second read's QC-fail flag must transfer. Before the split halves' + /// pairing bits were cleared, the second half (still `PAIRED | + /// LAST_SEGMENT`) looked for a mapped R2, found none, and silently copied + /// nothing. + #[test] + fn split_pair_halves_receive_their_own_tags_and_qc_flag() { + use fgumi_raw_bam::flags::{ + FIRST_SEGMENT, LAST_SEGMENT, MATE_UNMAPPED, PAIRED, QC_FAIL, REVERSE, SUPPLEMENTARY, + UNMAPPED, + }; + let cfg = make_test_cfg(); + let tags = ZipperTags::from_tag_info(&cfg.tag_info); + let r1 = make_record_with_string_tag( + b"pe10", + PAIRED | FIRST_SEGMENT | UNMAPPED | MATE_UNMAPPED, + crate::sam::SamTag::RX, + b"AAAA", + ); + let r2 = make_record_with_string_tag( + b"pe10", + PAIRED | LAST_SEGMENT | UNMAPPED | MATE_UNMAPPED | QC_FAIL, + crate::sam::SamTag::RX, + b"CCCC", + ); + let pair = Template::from_records(vec![r1, r2]).expect("unmapped pair"); + let (first, second) = + crate::pipeline::steps::align::split_pair_into_singles(pair).expect("split"); + let flags_of = |t: &Template| t.records()[0].flags(); + assert_eq!(flags_of(&first), UNMAPPED, "first half keeps only non-pairing bits"); + assert_eq!(flags_of(&second), UNMAPPED | QC_FAIL, "second half keeps QC_FAIL"); + + // The aligner emitted both reads unpaired; the second has a supplementary. + let mapped_first = + Template::from_records(vec![make_record(b"pe10", 0)]).expect("mapped first half"); + let mapped_second = Template::from_records(vec![ + make_record(b"pe10", REVERSE), + make_record(b"pe10", SUPPLEMENTARY), + ]) + .expect("mapped second half"); + + let zb = ZipperBatch { + serial: 0, + mapped: vec![mapped_first, mapped_second], + unmapped: BamTemplateBatch::new(0, vec![first, second]), + }; + let out = merge_zipper_batch(zb, &cfg, &tags).expect("merge ok"); + + let rx_of = |rec: &fgumi_raw_bam::RawRecord| { + fgumi_raw_bam::tags::find_string_tag( + fgumi_raw_bam::fields::aux_data_slice(rec), + crate::sam::SamTag::RX, + ) + .map(<[u8]>::to_vec) + }; + let first_out = &out.templates()[0].records; + let second_out = &out.templates()[1].records; + assert_eq!(rx_of(&first_out[0]), Some(b"AAAA".to_vec()), "first half gets R1's RX"); + assert_eq!(second_out.len(), 2, "second half keeps its supplementary"); + for rec in second_out { + assert_eq!(rx_of(rec), Some(b"CCCC".to_vec()), "second half gets R2's RX"); + assert_ne!(rec.flags() & QC_FAIL, 0, "second half gets R2's QC-fail flag"); + } + assert_eq!(first_out[0].flags() & QC_FAIL, 0, "first half's QC-pass status is unchanged"); + } + /// A `ZipperBatch` whose mapped and unmapped halves differ in length is a /// structural invariant violation: `merge_zipper_batch` must hard-error /// (release-safe) rather than let `zip` silently truncate to the shorter diff --git a/src/lib/pipeline/steps/align/mod.rs b/src/lib/pipeline/steps/align/mod.rs index f97725e85..f1d22ab98 100644 --- a/src/lib/pipeline/steps/align/mod.rs +++ b/src/lib/pipeline/steps/align/mod.rs @@ -17,7 +17,7 @@ //! //! This module holds the trait, its wiring context/result types, the //! backend-agnostic `ZipperBatch`, the subprocess backend's `InFlightGate` -//! byte-budget gate, and the header helpers (`validate_sq_consistency`, +//! byte-budget gate, the mid-pair split helper, and the header helpers (`validate_sq_consistency`, //! `merge_aligner_header`). pub(crate) mod merge; @@ -297,8 +297,12 @@ pub(crate) trait AlignBackend: Send + 'static { /// `aligner` module does not depend on the pipeline's align stage. pub(crate) fn backend_for(resolved: ResolvedAligner) -> Box { match resolved.backend { - ResolvedBackend::Subprocess { command } => { - Box::new(subprocess::SubprocessBackend { command, chunk_size: resolved.chunk_size }) + ResolvedBackend::Subprocess { command, accept_mid_pair_split } => { + Box::new(subprocess::SubprocessBackend { + command, + chunk_size: resolved.chunk_size, + accept_mid_pair_split, + }) } } } @@ -354,6 +358,56 @@ pub(crate) fn validate_sq_consistency(partial: &Header, aligner: &Header) -> io: Ok(()) } +/// The pairing bits [`split_pair_into_singles`] clears on each half's unmapped +/// record: everything that marks it as one segment of a pair. `QC_FAIL`, +/// `REVERSE` and the rest are kept. +const SPLIT_HALF_CLEARED_FLAGS: u16 = fgumi_raw_bam::flags::PAIRED + | fgumi_raw_bam::flags::PROPER_PAIR + | fgumi_raw_bam::flags::MATE_UNMAPPED + | fgumi_raw_bam::flags::MATE_REVERSE + | fgumi_raw_bam::flags::FIRST_SEGMENT + | fgumi_raw_bam::flags::LAST_SEGMENT; + +/// Split a two-primary-read template into two single-record templates, in record +/// order: bwa's mid-pair split. With mixed SE/PE input a `-K` +/// chunk boundary can fall between a pair's two reads, and bwa then aligns them as +/// two unpaired reads. The pair's unmapped template is split the same way so each +/// half zips with its own read's alignment. Errors if the template carries records +/// beyond its two primaries (secondaries make the split ambiguous — out of +/// contract for valid unmapped input). +/// +/// Each half's record has its pairing bits ([`SPLIT_HALF_CLEARED_FLAGS`]) cleared, +/// so it is an unpaired read like the aligner's record for it. The zipper merge +/// picks the mapped segment to copy tags and the QC-fail flag onto from the +/// *unmapped* record's `PAIRED`/`FIRST_SEGMENT` bits; left set, the second half +/// (still `PAIRED | LAST_SEGMENT`) would look for a mapped R2, find none (the +/// aligner emitted it unpaired, which files as R1), and silently copy nothing. +/// The merged record's own flags come from the aligner, so clearing these bits +/// on the unmapped side changes nothing else in the output. +pub(crate) fn split_pair_into_singles(template: Template) -> io::Result<(Template, Template)> { + // Reject anything but exactly two primaries before consuming the template, + // so the error path reads `template.name()` by reference and the success + // path never allocates a name copy. + if template.records().len() != 2 { + return Err(io::Error::other(format!( + "align-and-merge: cannot split template '{name}' across a mid-pair -K chunk \ + boundary: it carries {n} records, not exactly two primaries (secondary/supplementary \ + records are out of contract for mixed single/paired input).", + name = String::from_utf8_lossy(template.name()), + n = template.records().len(), + ))); + } + let mut records = template.into_records(); + let second = records.pop().expect("len checked == 2"); + let first = records.pop().expect("len checked == 2"); + let to_template = |mut record: fgumi_raw_bam::RawRecord| { + record.set_flags(record.flags() & !SPLIT_HALF_CLEARED_FLAGS); + Template::from_records(vec![record]) + .map_err(|e| io::Error::other(format!("align-and-merge: split-half template: {e:#}"))) + }; + Ok((to_template(first)?, to_template(second)?)) +} + /// Merge aligner-emitted header lines into the partial output header. /// /// The aligner contributes: diff --git a/src/lib/pipeline/steps/align/subprocess.rs b/src/lib/pipeline/steps/align/subprocess.rs index 7baf199db..67276279c 100644 --- a/src/lib/pipeline/steps/align/subprocess.rs +++ b/src/lib/pipeline/steps/align/subprocess.rs @@ -83,7 +83,7 @@ use crate::pipeline::steps::align::merge::MergeAlignedStep; use crate::pipeline::steps::align::{ AlignBackend, AlignWired, AlignWiringCtx, ConsumerGoneGuard, InFlightGate, ZipperBatch, in_flight_budget_for_chunk_size, is_primary_for_alignment, merge_aligner_header, - no_primary_records_message, validate_sq_consistency, + no_primary_records_message, split_pair_into_singles, validate_sq_consistency, }; use crate::pipeline::steps::types::BamTemplateBatch; use crate::template::Template; @@ -156,6 +156,13 @@ pub(crate) struct SubprocessConfig { /// unmapped reads in RAM (issue #382). Derive it from the aligner's /// `-K` chunk size via [`in_flight_budget_for_chunk_size`]. pub(crate) in_flight_unmapped_budget: u64, + + /// Whether a same-queryname output group with two unpaired primaries is + /// accepted as bwa's mid-pair split (see [`AlignedGroup::MidPairSplit`]). + /// Set only for the subprocess presets (`mem -p -K`); otherwise the group + /// fails `Template::from_records`'s "multiple primaries" check, as it did + /// before the split was recognized. + pub(crate) accept_mid_pair_split: bool, } // ────────────────────────────────────────────────────────────────────────── @@ -580,8 +587,10 @@ fn reader_loop_inner( .read_header() .map_err(|e| io::Error::other(format!("aligner BAM header: {e}")))?; let bgzf = bam_reader.into_inner(); - let stream = - TemplateStream::Bam(BamTemplateStream::new(fgumi_raw_bam::RawBamReader::new(bgzf))); + let stream = TemplateStream::Bam(BamTemplateStream::new( + fgumi_raw_bam::RawBamReader::new(bgzf), + cfg.accept_mid_pair_split, + )); (stream, aligner_header) } else { let buffered = BufReader::with_capacity(STDOUT_READER_BUF_BYTES, stdout); @@ -590,7 +599,11 @@ fn reader_loop_inner( .read_header() .map_err(|e| io::Error::other(format!("aligner SAM header: {e}")))?; let header_arc = Arc::new(aligner_header.clone()); - let stream = TemplateStream::Sam(SamTemplateStream::new(sam_reader, header_arc)); + let stream = TemplateStream::Sam(SamTemplateStream::new( + sam_reader, + header_arc, + cfg.accept_mid_pair_split, + )); (stream, aligner_header) }; @@ -656,6 +669,9 @@ fn reader_loop_inner( } let mut mapped_templates: Vec