Skip to content

fix(filter,clip): require query-grouped input (FILT3-02, CLIP3-05) #517

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nh13 merged 1 commit into from
Jul 13, 2026
Merged

fix(filter,clip): require query-grouped input (FILT3-02, CLIP3-05) #517

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions crates/fgumi-sam/src/builder.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
163 changes: 135 additions & 28 deletions crates/fgumi-sam/src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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`
Expand All @@ -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::<noodles::sam::header::record::value::map::Header>::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.
Expand All @@ -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).
Expand Down Expand Up @@ -531,6 +593,7 @@ pub fn buf_value_to_smallest_signed_int(value: &BufValue) -> Option<BufValue> {
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() {
Expand Down Expand Up @@ -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()
// =========================================================================
Expand Down
26 changes: 26 additions & 0 deletions src/lib/commands/clip.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)?);

Expand Down Expand Up @@ -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)?;

Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down
Loading
Loading