From c676f1b09a734e49c41661cdc06a5f69a268b2bb Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Thu, 9 Jul 2026 10:14:12 -0700 Subject: [PATCH] fix(filter,clip): require query-grouped input (FILT3-02, CLIP3-05) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit fgbio's FilterConsensusReads and ClipBam both call Bams.requireQueryGrouped before processing; fgumi did no such check. On coordinate-sorted input a template's mates scatter, every read degrades to its own single-read "template", and pair/template logic is silently wrong with a success exit (filter drops all reads; clip no-ops pair clipping). Add a query-grouped guard matching fgbio's isQueryGrouped (SO:queryname || GO:query) — weaker than the consensus template-coordinate guard, so a plain queryname sort is accepted: - fgumi-sam: is_query_grouped() predicate, header_as_queryname(), and a SamBuilder::set_queryname_sort_order() test helper. - commands::common: require_query_grouped() with fgbio-style diagnostics that echo the found SO/GO. - Wire it into filter::execute and both clip::execute paths, guarding the input header before clip rewrites SO for the output. Real-tool parity (adversarial coordinate-sorted fixture with scattered mates): before: fgbio errors; fgumi exits 0 (filter "kept 0, rejected 4"; clip no-op). after: fgumi errors like fgbio ("not queryname sorted or query grouped, found: SO:coordinate GO:none"); both accept query-grouped input (clip 4->4 byte parity). Test fixtures that fed bare-header BAMs updated to stamp SO:queryname. Tracker: reports/2026-07-09-fgumi-final-audit-burndown-tracker.md (W2) --- crates/fgumi-sam/src/builder.rs | 11 ++ crates/fgumi-sam/src/lib.rs | 163 +++++++++++++++++++---- src/lib/commands/clip.rs | 26 ++++ src/lib/commands/common.rs | 103 ++++++++++++++ src/lib/commands/filter.rs | 16 ++- tests/integration/test_clip_command.rs | 94 ++++++++++++- tests/integration/test_filter_command.rs | 75 +++++++++-- 7 files changed, 441 insertions(+), 47 deletions(-) diff --git a/crates/fgumi-sam/src/builder.rs b/crates/fgumi-sam/src/builder.rs index 503ca6103..ba578773e 100644 --- a/crates/fgumi-sam/src/builder.rs +++ b/crates/fgumi-sam/src/builder.rs @@ -256,6 +256,17 @@ impl SamBuilder { self } + /// Stamps queryname sort order (`SO:queryname`) onto the header. + /// + /// Template-based commands (`filter`, `clip`) require query-grouped input + /// (fgbio `Bams.requireQueryGrouped`), so tests that exercise them end-to-end + /// should call this before writing the input BAM. The builder's per-pair + /// insertion order keeps a template's reads adjacent, satisfying the grouping. + pub fn set_queryname_sort_order(&mut self) -> &mut Self { + self.header = crate::header_as_queryname(&self.header); + self + } + /// Returns the next sequential name. fn next_name(&self) -> String { format!("{:04}", self.counter.fetch_add(1, Ordering::SeqCst)) diff --git a/crates/fgumi-sam/src/lib.rs b/crates/fgumi-sam/src/lib.rs index a1f65511f..48e6fe22f 100644 --- a/crates/fgumi-sam/src/lib.rs +++ b/crates/fgumi-sam/src/lib.rs @@ -251,6 +251,42 @@ pub fn is_query_grouped_unsorted(header: &Header) -> bool { } } +/// Checks if a header is queryname sorted or query grouped, matching fgbio's +/// `Bams.isQueryGrouped` (`header.getSortOrder == queryname || header.getGroupOrder == query`). +/// +/// When this holds, all reads with the same query name are adjacent, so +/// template-based commands (filter, clip) can group a template's mates by +/// adjacency. This is *weaker* than [`is_template_coordinate_sorted`] / +/// [`is_query_grouped_unsorted`]: a plain `SO:queryname` file (no `GO`) qualifies +/// here but not there. Template-coordinate input (`SO:unsorted GO:query …`) also +/// qualifies via its `GO:query`. +/// +/// # Arguments +/// +/// * `header` - SAM header to check +/// +/// # Returns +/// +/// `true` if the header advertises `SO:queryname` or `GO:query`, `false` otherwise. +#[must_use] +pub fn is_query_grouped(header: &Header) -> bool { + if let Some(hdr_map) = header.header() { + let other_fields = hdr_map.other_fields(); + + // SO:queryname (a full queryname sort is also query grouped). + let so_queryname = + other_fields.get(b"SO").is_some_and(|so| <_ as AsRef<[u8]>>::as_ref(so) == QUERY_NAME); + + // GO:query (query grouped without a queryname sort — includes template-coordinate). + let go_query = + other_fields.get(b"GO").is_some_and(|go| <_ as AsRef<[u8]>>::as_ref(go) == b"query"); + + so_queryname || go_query + } else { + false + } +} + /// Checks if a BAM file is queryname sorted and logs a warning if not. /// /// This function validates that the input file has the correct sort order for @@ -283,6 +319,33 @@ pub fn check_sort(header: &Header, path: &Path, name: &str) { } } +/// Rebuilds `header` with a replaced `@HD` (header) map, preserving its read +/// groups, reference sequences, programs, and comments. +/// +/// Shared by the `header_as_*` sort-order stampers, which differ only in how +/// they mutate the `@HD` map before rebuilding. +fn rebuild_header_with_hd( + header: &Header, + header_map: noodles::sam::header::record::value::Map< + noodles::sam::header::record::value::map::Header, + >, +) -> Header { + let mut builder = Header::builder(); + for (name, rg) in header.read_groups() { + builder = builder.add_read_group(name.clone(), rg.clone()); + } + for (name, reference) in header.reference_sequences() { + builder = builder.add_reference_sequence(name.clone(), reference.clone()); + } + for (id, pg) in header.programs().as_ref() { + builder = builder.add_program(id.clone(), pg.clone()); + } + for comment in header.comments() { + builder = builder.add_comment(comment.clone()); + } + builder.set_header(header_map).build() +} + /// Returns a copy of `header` with its sort-order metadata cleared. /// /// The `SO` tag is set to `unsorted` and any `GO` (group order) or `SS` @@ -305,20 +368,32 @@ pub fn header_as_unsorted(header: &Header) -> Header { other_fields.shift_remove(b"GO"); other_fields.shift_remove(b"SS"); - let mut builder = Header::builder(); - for (name, rg) in header.read_groups() { - builder = builder.add_read_group(name.clone(), rg.clone()); - } - for (name, reference) in header.reference_sequences() { - builder = builder.add_reference_sequence(name.clone(), reference.clone()); - } - for (id, pg) in header.programs().as_ref() { - builder = builder.add_program(id.clone(), pg.clone()); - } - for comment in header.comments() { - builder = builder.add_comment(comment.clone()); - } - builder.set_header(header_map).build() + rebuild_header_with_hd(header, header_map) +} + +/// Returns a copy of `header` advertising queryname sort order (`SO:queryname`), +/// with any `GO`/`SS` tags removed. +/// +/// This is the weakest order that satisfies [`is_query_grouped`], so it is the +/// natural input order for template-based commands (filter, clip). Used by the +/// test [`crate::builder::SamBuilder`] to stamp query-grouped fixtures. +#[must_use] +pub fn header_as_queryname(header: &Header) -> Header { + use bstr::BString; + use noodles::sam::header::record::value::Map; + use noodles::sam::header::record::value::map::header::tag::SORT_ORDER; + + let mut header_map = header + .header() + .cloned() + .unwrap_or_else(Map::::default); + + let other_fields = header_map.other_fields_mut(); + other_fields.insert(SORT_ORDER, BString::from(QUERY_NAME)); + other_fields.shift_remove(b"GO"); + other_fields.shift_remove(b"SS"); + + rebuild_header_with_hd(header, header_map) } /// Returns a copy of `header` advertising template-coordinate sort order. @@ -344,20 +419,7 @@ pub fn header_as_template_coordinate(header: &Header) -> Header { other_fields.insert(GROUP_ORDER, BString::from("query")); other_fields.insert(SUBSORT_ORDER, BString::from("template-coordinate")); - let mut builder = Header::builder(); - for (name, rg) in header.read_groups() { - builder = builder.add_read_group(name.clone(), rg.clone()); - } - for (name, reference) in header.reference_sequences() { - builder = builder.add_reference_sequence(name.clone(), reference.clone()); - } - for (id, pg) in header.programs().as_ref() { - builder = builder.add_program(id.clone(), pg.clone()); - } - for comment in header.comments() { - builder = builder.add_comment(comment.clone()); - } - builder.set_header(header_map).build() + rebuild_header_with_hd(header, header_map) } /// Reverses a `BufValue` (array or string). @@ -531,6 +593,7 @@ pub fn buf_value_to_smallest_signed_int(value: &BufValue) -> Option { mod tests { use super::*; use noodles::sam::alignment::record_buf::data::field::value::Array; + use rstest::rstest; #[test] fn test_reverse_buf_value_int8_array() { @@ -996,6 +1059,50 @@ mod tests { assert!(is_query_grouped_unsorted(&header)); } + // ========================================================================= + // Tests for is_query_grouped() — matches fgbio Bams.isQueryGrouped. + // The final column contrasts the weaker is_query_grouped_unsorted (which + // additionally requires SO:unsorted + GO:query) so the two are pinned + // side-by-side. + // ========================================================================= + + #[rstest] + // queryname sort qualifies for is_query_grouped but not the unsorted variant. + #[case::queryname("queryname", None, None, true, false)] + // GO:query alone (SO:unsorted) qualifies for both. + #[case::go_query("unsorted", Some("query"), None, true, true)] + // Template-coordinate qualifies for both via GO:query. + #[case::template_coordinate("unsorted", Some("query"), Some("template-coordinate"), true, true)] + // SO:coordinate GO:none — the FILT3-02/CLIP3-05 footgun input. + #[case::coordinate("coordinate", None, None, false, false)] + // Bare SO:unsorted with no GO is not query grouped (mates may scatter). + #[case::bare_unsorted("unsorted", None, None, false, false)] + fn test_is_query_grouped( + #[case] so: &str, + #[case] go: Option<&str>, + #[case] ss: Option<&str>, + #[case] expected_query_grouped: bool, + #[case] expected_query_grouped_unsorted: bool, + ) { + let header = create_template_coord_header(so, go, ss); + assert_eq!(is_query_grouped(&header), expected_query_grouped); + assert_eq!(is_query_grouped_unsorted(&header), expected_query_grouped_unsorted); + } + + #[test] + fn test_header_as_queryname_stamps_so_and_strips_go_ss() { + // Start from a template-coordinate header (SO:unsorted GO:query SS:…) and + // confirm header_as_queryname rewrites SO:queryname and drops GO/SS. + let header = + create_template_coord_header("unsorted", Some("query"), Some("template-coordinate")); + let stamped = header_as_queryname(&header); + assert!(is_query_grouped(&stamped)); + assert!(is_sorted(&stamped, QUERY_NAME)); + let other = stamped.header().expect("has @HD").other_fields(); + assert!(other.get(b"GO").is_none(), "GO must be stripped"); + assert!(other.get(b"SS").is_none(), "SS must be stripped"); + } + // ========================================================================= // Tests for header_as_template_coordinate() // ========================================================================= diff --git a/src/lib/commands/clip.rs b/src/lib/commands/clip.rs index c9b337ab4..7cdad6fb0 100644 --- a/src/lib/commands/clip.rs +++ b/src/lib/commands/clip.rs @@ -237,6 +237,14 @@ impl Command for Clip { // Read header for the 7-step pipeline (supports stdin) let (reader, header) = create_bam_reader_for_pipeline(&self.io.input)?; + // CLIP3-05: fgbio's ClipBam calls Bams.requireQueryGrouped. Clipping is + // template-based (pair clip, overlap, past-mate, mate-fix), so coordinate- + // sorted input silently mis-groups mates. Guard the *input* header. + crate::commands::common::require_query_grouped( + &header, + &self.io.input.display().to_string(), + )?; + // Load reference (always required for clip) let reference = Arc::new(ReferenceReader::new(&self.reference)?); @@ -275,6 +283,12 @@ impl Clip { let reader_threads = self.threading.num_threads(); let (reader, header) = create_raw_bam_reader(&self.io.input, reader_threads)?; + // CLIP3-05: require query-grouped input (fgbio Bams.requireQueryGrouped). + crate::commands::common::require_query_grouped( + &header, + &self.io.input.display().to_string(), + )?; + // Add @PG record with PP chaining to input's last program let header = crate::commands::common::add_pg_record(header, command_line)?; @@ -2003,6 +2017,7 @@ mod tests { // Create overlapping read pair let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2053,6 +2068,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2102,6 +2118,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2151,6 +2168,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2200,6 +2218,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2249,6 +2268,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2299,6 +2319,7 @@ mod tests { let metrics_path = dir.path().join("metrics.txt"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2350,6 +2371,7 @@ mod tests { // Create a fragment (unpaired) read let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_frag() .name("frag1") @@ -2399,6 +2421,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2449,6 +2472,7 @@ mod tests { let metrics_path = dir.path().join("metrics.txt"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2499,6 +2523,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") @@ -2552,6 +2577,7 @@ mod tests { let output_path = dir.path().join("output.bam"); let mut builder = SamBuilder::with_single_ref("chr1", 200); + builder.set_queryname_sort_order(); // clip requires query-grouped input (CLIP3-05) let _ = builder .add_pair() .name("read1") diff --git a/src/lib/commands/common.rs b/src/lib/commands/common.rs index b17c462f9..2e29bb040 100644 --- a/src/lib/commands/common.rs +++ b/src/lib/commands/common.rs @@ -172,6 +172,65 @@ pub fn check_consensus_sort_order(header: &Header, source: &str) -> anyhow::Resu ); } +/// Requires that the input header is queryname sorted or query grouped, matching +/// fgbio's `Bams.requireQueryGrouped` (used by `FilterConsensusReads` and `ClipBam`). +/// +/// Template-based commands group a template's reads by *adjacency*. On coordinate- +/// sorted input the mates scatter, every read degrades to its own single-read +/// "template", and pair-level logic (overlap clipping, mate-info fix, +/// both-primaries-pass) is silently wrong with a success exit. fgbio hard-fails +/// here; fgumi must too. +/// +/// This is the query-grouped guard, weaker than [`check_consensus_sort_order`]: +/// a plain `SO:queryname` file (no template-coordinate sub-sort) is accepted here +/// but rejected there. Use this for filter/clip; use `check_consensus_sort_order` +/// for consensus callers. +/// +/// Reads only the header, so it is safe to call on a stdin-backed reader after the +/// header has been read. +/// +/// # Arguments +/// +/// * `header` - the already-read input BAM header +/// * `source` - a human-readable description of the input (path or ``) +/// +/// # Errors +/// +/// Returns an error if the header is neither queryname sorted (`SO:queryname`) nor +/// query grouped (`GO:query`). +pub fn require_query_grouped(header: &Header, source: &str) -> anyhow::Result<()> { + if fgumi_sam::is_query_grouped(header) { + return Ok(()); + } + let (so, go, ss) = header_sort_and_group_order(header); + // Mirror fgbio's requireQueryGrouped: append " SS:{ss}" only when present. + let ss_suffix = ss.map_or_else(String::new, |ss| format!(" SS:{ss}")); + anyhow::bail!( + "File {source} was not queryname sorted or query grouped, found: SO:{so} GO:{go}{ss_suffix}. \ + A template's reads must be adjacent, so the input must advertise SO:queryname \ + or GO:query.\n\nSort it first, e.g.:\n \ + fgumi sort -i input.bam -o sorted.bam --order queryname" + ); +} + +/// Returns the header's declared sort order (`SO`, default `unsorted` per htsjdk), +/// group order (`GO`, default `none`), and sub-sort (`SS`, `None` when absent) as +/// display strings for diagnostics — mirroring fgbio's `requireQueryGrouped`. +fn header_sort_and_group_order(header: &Header) -> (String, String, Option) { + let Some(hdr_map) = header.header() else { + return ("unsorted".to_string(), "none".to_string(), None); + }; + let other = hdr_map.other_fields(); + let read = |tag: &[u8; 2]| -> Option { + other.get(tag).map(|v| String::from_utf8_lossy(<_ as AsRef<[u8]>>::as_ref(v)).into_owned()) + }; + ( + read(b"SO").unwrap_or_else(|| "unsorted".to_string()), + read(b"GO").unwrap_or_else(|| "none".to_string()), + read(b"SS"), + ) +} + /// EM-Seq methylation-aware consensus calling options. #[derive(Debug, Clone, Default, Args)] pub struct EmSeqOptions { @@ -1054,6 +1113,50 @@ mod tests { } } + /// FILT3-02 / CLIP3-05: `require_query_grouped` mirrors fgbio's + /// `Bams.requireQueryGrouped` (`isQueryGrouped = SO:queryname || GO:query`). + /// Error cases assert the found-order echo and the fgbio-style message. + #[rstest] + // A plain queryname sort qualifies (unlike the stricter consensus guard). + #[case::queryname("@HD\tVN:1.6\tSO:queryname\n", true, &[])] + // GO:query alone qualifies. + #[case::go_query("@HD\tVN:1.6\tSO:unsorted\tGO:query\n", true, &[])] + // Template-coordinate consensus output qualifies via GO:query. + #[case::template_coordinate( + "@HD\tVN:1.6\tSO:unsorted\tGO:query\tSS:template-coordinate\n", + true, + &[] + )] + // The FILT3-02/CLIP3-05 footgun: coordinate-sorted input scatters mates. + #[case::coordinate( + "@HD\tVN:1.6\tSO:coordinate\n", + false, + &["foo.bam", "queryname sorted or query grouped", "SO:coordinate", "GO:none"] + )] + // A rejected header with an SS sub-sort echoes it like fgbio (" SS:{ss}"). + #[case::coordinate_with_ss( + "@HD\tVN:1.6\tSO:coordinate\tSS:coordinate:natural\n", + false, + &["SO:coordinate", "GO:none", "SS:coordinate:natural"] + )] + // No SO/GO at all (htsjdk default unsorted, GO none) → rejected. + #[case::bare_unsorted("@HD\tVN:1.6\n", false, &["queryname sorted or query grouped"])] + fn test_require_query_grouped( + #[case] header_str: &str, + #[case] expect_ok: bool, + #[case] expected_substrings: &[&str], + ) { + let header: Header = header_str.parse().expect("parse"); + let result = require_query_grouped(&header, "foo.bam"); + assert_eq!(result.is_ok(), expect_ok, "for header {header_str:?}"); + if let Err(err) = result { + let msg = err.to_string(); + for sub in expected_substrings { + assert!(msg.contains(sub), "message {msg:?} missing {sub:?}"); + } + } + } + #[test] fn test_none_is_single_threaded() { let opts = ThreadingOptions::none(); diff --git a/src/lib/commands/filter.rs b/src/lib/commands/filter.rs index d6192f4df..b98528de0 100644 --- a/src/lib/commands/filter.rs +++ b/src/lib/commands/filter.rs @@ -337,6 +337,14 @@ impl Command for Filter { self.io.pipeline_reader_opts(), )?; + // FILT3-02: fgbio's FilterConsensusReads calls Bams.requireQueryGrouped. + // Filtering is template-based, so coordinate-sorted input silently scatters + // mates and corrupts the both-primaries-pass logic. Reject it like fgbio. + crate::commands::common::require_query_grouped( + &header, + &self.io.input.display().to_string(), + )?; + // Add @PG record with PP chaining to input's last program let header = crate::commands::common::add_pg_record(header, command_line)?; @@ -2438,16 +2446,20 @@ mod tests { } /// Helper: build a minimal chr1 noodles `Header` for use with the noodles BAM writer. + /// + /// Stamped `SO:queryname` because `filter` now requires query-grouped input + /// (FILT3-02, matching fgbio `Bams.requireQueryGrouped`). fn test_bam_header() -> noodles::sam::Header { use noodles::sam::header::record::value::Map; use noodles::sam::header::record::value::map::ReferenceSequence; use std::num::NonZeroUsize; - noodles::sam::Header::builder() + let header = noodles::sam::Header::builder() .add_reference_sequence( "chr1", Map::::new(NonZeroUsize::new(1000).expect("1000 is non-zero")), ) - .build() + .build(); + fgumi_sam::header_as_queryname(&header) } /// Write a list of raw records as a BAM file with a single chr1 reference. diff --git a/tests/integration/test_clip_command.rs b/tests/integration/test_clip_command.rs index f0f5968f0..11cfe05a4 100644 --- a/tests/integration/test_clip_command.rs +++ b/tests/integration/test_clip_command.rs @@ -13,26 +13,106 @@ use noodles::bam; use noodles::sam::alignment::RecordBuf; use noodles::sam::alignment::io::Write as AlignmentWrite; use noodles::sam::alignment::record::cigar::op::Kind as CigarKind; +use rstest::rstest; use std::fs; use std::path::Path; use tempfile::TempDir; -use crate::helpers::bam_generator::{create_minimal_header, create_test_reference, to_record_buf}; +use crate::helpers::bam_generator::{ + create_coordinate_sorted_header, create_minimal_header, create_test_reference, to_record_buf, +}; -/// Create a BAM with paired reads. -fn create_paired_bam(path: &Path, pairs: Vec<(RawRecord, RawRecord)>) { - let header = create_minimal_header("chr1", 10000); +/// Create a BAM with paired reads using the given header. +fn create_paired_bam_with_header( + path: &Path, + header: &noodles::sam::Header, + pairs: Vec<(RawRecord, RawRecord)>, +) { let mut writer = bam::io::Writer::new(fs::File::create(path).expect("Failed to create BAM file")); - writer.write_header(&header).expect("Failed to write header"); + writer.write_header(header).expect("Failed to write header"); for (r1, r2) in pairs { - writer.write_alignment_record(&header, &to_record_buf(&r1)).expect("Failed to write R1"); - writer.write_alignment_record(&header, &to_record_buf(&r2)).expect("Failed to write R2"); + writer.write_alignment_record(header, &to_record_buf(&r1)).expect("Failed to write R1"); + writer.write_alignment_record(header, &to_record_buf(&r2)).expect("Failed to write R2"); } writer.try_finish().expect("Failed to finish BAM"); } +/// Create a BAM with paired reads and a query-grouped (template-coordinate) header. +fn create_paired_bam(path: &Path, pairs: Vec<(RawRecord, RawRecord)>) { + let header = create_minimal_header("chr1", 10000); + create_paired_bam_with_header(path, &header, pairs); +} + +/// CLIP3-05: clip must reject coordinate-sorted (non-query-grouped) input, +/// matching fgbio's `Bams.requireQueryGrouped`. On coordinate-sorted input mates +/// scatter, so pair clip / overlap / mate-fix silently no-op — hard-fail instead. +/// +/// Exercised on both `Clip::execute` paths (the single-threaded fast path and the +/// `--threads` pipeline) since the guard is wired into each. +#[rstest] +#[case::single_threaded(&[])] +#[case::threaded(&["--threads", "2"])] +fn test_clip_rejects_coordinate_sorted_input(#[case] extra_args: &[&str]) { + let temp_dir = TempDir::new().unwrap(); + let input_bam = temp_dir.path().join("input.bam"); + let output_bam = temp_dir.path().join("output.bam"); + let ref_path = create_test_reference(temp_dir.path()); + + let r1 = { + let mut b = SamBuilder::new(); + b.read_name(b"read1") + .sequence(b"ACGTACGT") + .qualities(&[30; 8]) + .flags(flags::PAIRED | flags::FIRST_SEGMENT | flags::MATE_REVERSE) + .ref_id(0) + .pos(99) + .mapq(60) + .cigar_ops(&[8 << 4]) + .mate_ref_id(0) + .mate_pos(103) + .template_length(12); + b.build() + }; + let r2 = { + let mut b = SamBuilder::new(); + b.read_name(b"read1") + .sequence(b"ACGTACGT") + .qualities(&[30; 8]) + .flags(flags::PAIRED | flags::LAST_SEGMENT | flags::REVERSE) + .ref_id(0) + .pos(103) + .mapq(60) + .cigar_ops(&[8 << 4]) + .mate_ref_id(0) + .mate_pos(99) + .template_length(-12); + b.build() + }; + let header = create_coordinate_sorted_header("chr1", 10000); + create_paired_bam_with_header(&input_bam, &header, vec![(r1, r2)]); + + let mut args = vec![ + "clip", + "--input", + input_bam.to_str().unwrap(), + "--output", + output_bam.to_str().unwrap(), + "--ref", + ref_path.to_str().unwrap(), + "--clip-overlapping-reads", + ]; + args.extend_from_slice(extra_args); + let cmd = Clip::try_parse_from(args).expect("failed to parse clip args"); + + let err = cmd.execute("fgumi clip").expect_err("must reject coordinate-sorted input"); + assert!( + err.to_string().contains("queryname sorted or query grouped"), + "unexpected error message: {err}" + ); +} + /// Test basic clip command with fixed-end clipping. #[test] fn test_clip_command_basic() { diff --git a/tests/integration/test_filter_command.rs b/tests/integration/test_filter_command.rs index be3bf2c27..6080a1799 100644 --- a/tests/integration/test_filter_command.rs +++ b/tests/integration/test_filter_command.rs @@ -16,7 +16,9 @@ use std::fs; use std::path::Path; use tempfile::TempDir; -use crate::helpers::bam_generator::{create_minimal_header, create_test_reference}; +use crate::helpers::bam_generator::{ + create_coordinate_sorted_header, create_minimal_header, create_test_reference, +}; /// Convert a raw BAM record to a `RecordBuf` using the given header. /// @@ -31,24 +33,34 @@ fn to_record_buf( .expect("raw_record_to_record_buf should succeed in test") } -/// Create a consensus BAM with records that have consensus tags. -fn create_consensus_bam(path: &Path, records: Vec) { - let header = create_minimal_header("chr1", 10000); +/// Create a consensus BAM with the given header and consensus-tagged records. +fn create_consensus_bam_with_header( + path: &Path, + header: &noodles::sam::Header, + records: Vec, +) { let mut writer = bam::io::Writer::new(fs::File::create(path).expect("Failed to create BAM file")); - writer.write_header(&header).expect("Failed to write header"); + writer.write_header(header).expect("Failed to write header"); for record in records { writer - .write_alignment_record(&header, &to_record_buf(&record, &header)) + .write_alignment_record(header, &to_record_buf(&record, header)) .expect("Failed to write record"); } writer.try_finish().expect("Failed to finish BAM"); } -/// Write a small mapped-consensus BAM (two reads with CD/CE per-base tags that -/// pass `filter`). Shared by `test_filter_command` and the stdin parity test. -pub(crate) fn write_filter_consensus_bam(path: &Path) { +/// Create a consensus BAM with a query-grouped (template-coordinate) header. +fn create_consensus_bam(path: &Path, records: Vec) { + let header = create_minimal_header("chr1", 10000); + create_consensus_bam_with_header(path, &header, records); +} + +/// Two mapped-consensus records with CD/CE per-base tags that pass `filter` +/// (`cons1` depth 10, `cons2` depth 5). Shared by the query-grouped writer and +/// the coordinate-sorted guard test. +fn filter_consensus_records() -> Vec { let r1 = { let mut b = RawSamBuilder::new(); b.read_name(b"cons1") @@ -75,7 +87,13 @@ pub(crate) fn write_filter_consensus_bam(path: &Path) { b.add_array_u16(SamTag::CD_BASES, &[5; 8]).add_array_u16(SamTag::CE_BASES, &[0; 8]); b.build() }; - create_consensus_bam(path, vec![r1, r2]); + vec![r1, r2] +} + +/// Write a small mapped-consensus BAM (two reads with CD/CE per-base tags that +/// pass `filter`). Shared by `test_filter_command` and the stdin parity test. +pub(crate) fn write_filter_consensus_bam(path: &Path) { + create_consensus_bam(path, filter_consensus_records()); } /// Test basic filter command with passing reads. @@ -114,6 +132,43 @@ fn test_filter_command_basic() { assert_eq!(count, 2, "Both reads should be present in output"); } +/// FILT3-02: filter must reject coordinate-sorted (non-query-grouped) input, +/// matching fgbio's `Bams.requireQueryGrouped`. On coordinate-sorted input a +/// template's mates scatter, every read degrades to its own "template", and the +/// both-primaries-pass logic silently corrupts the result — so we hard-fail. +#[test] +fn test_filter_rejects_coordinate_sorted_input() { + let temp_dir = TempDir::new().unwrap(); + let input_bam = temp_dir.path().join("input.bam"); + let output_bam = temp_dir.path().join("output.bam"); + let ref_path = create_test_reference(temp_dir.path()); + + // Same consensus records as write_filter_consensus_bam, but a coordinate header. + let header = create_coordinate_sorted_header("chr1", 10000); + create_consensus_bam_with_header(&input_bam, &header, filter_consensus_records()); + + let cmd = Filter::try_parse_from([ + "filter", + "--input", + input_bam.to_str().unwrap(), + "--output", + output_bam.to_str().unwrap(), + "--ref", + ref_path.to_str().unwrap(), + "--min-reads", + "1", + "--max-no-call-fraction", + "1.0", + "--compression-level", + "1", + ]) + .expect("failed to parse filter args"); + + let err = cmd.execute("fgumi filter").expect_err("must reject coordinate-sorted input"); + let msg = err.to_string(); + assert!(msg.contains("queryname sorted or query grouped"), "unexpected error message: {msg}"); +} + /// Test filter command with reads that fail due to low per-base depth. /// /// The filter masks bases where the per-base depth (cd tag) is below --min-reads.