From 6004b3c211a10d021c9b0e186fdf80d5df55ef0a Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Wed, 8 Jul 2026 23:23:26 -0700 Subject: [PATCH 1/2] fix(sort): write SS sub-sort tag as : for fgbio/samtools parity MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit fgbio (`SamOrder.applyTo`: `sortOrder.name() + ":" + ss`) and samtools write the `SS` header tag as `:` — e.g. `queryname:natural`, `unsorted:template-coordinate`. fgumi wrote the bare sub-sort (`natural`, `template-coordinate`), which diverges from both parity targets: fgbio's `SamOrder.apply` strips everything up to the first colon (`s.substring(s.indexOf(':') + 1)`), so a bare `SS:natural` is not re-recognized as a queryname order. (Template-coordinate happened to round-trip by luck.) `SortOrder::header_ss_tag` now returns the `:` form, and `create_output_header` uses it for both the queryname and template-coordinate branches. The `` prefix is the value written to the `SO` tag (`queryname` for queryname sub-sorts, `unsorted` for template-coordinate). No reader regression: `is_template_coordinate_sorted` already splits on the first colon and accepts both the prefixed and bare forms (mirroring fgbio), so fgumi-, fgbio-, and samtools-written headers are all accepted. There is no reader that inspects the queryname sub-sort value. R2-HDR-01. --- crates/fgumi-sort/src/external.rs | 2 +- crates/fgumi-sort/src/keys.rs | 26 ++++++++++++++++------ crates/fgumi-sort/src/lib.rs | 37 +++++++++++++++++++++++++++++-- src/lib/commands/sort.rs | 6 ++--- 4 files changed, 58 insertions(+), 13 deletions(-) diff --git a/crates/fgumi-sort/src/external.rs b/crates/fgumi-sort/src/external.rs index 2392fc262..0d0f32a91 100644 --- a/crates/fgumi-sort/src/external.rs +++ b/crates/fgumi-sort/src/external.rs @@ -3873,7 +3873,7 @@ mod tests { assert_eq!(<_ as AsRef<[u8]>>::as_ref(go), b"query"); let ss = fields.get(b"SS").expect("should have SS tag"); - assert_eq!(<_ as AsRef<[u8]>>::as_ref(ss), b"template-coordinate"); + assert_eq!(<_ as AsRef<[u8]>>::as_ref(ss), b"unsorted:template-coordinate"); } // ======================================================================== diff --git a/crates/fgumi-sort/src/keys.rs b/crates/fgumi-sort/src/keys.rs index e2294af02..92facbc03 100644 --- a/crates/fgumi-sort/src/keys.rs +++ b/crates/fgumi-sort/src/keys.rs @@ -196,14 +196,22 @@ impl SortOrder { } } - /// Get the SAM header sub-sort tag value. + /// Get the SAM header sub-sort (`SS`) tag value. + /// + /// The `SS` tag is written as `:`, matching fgbio + /// (`SamOrder.applyTo`: `sortOrder.name() + ":" + ss`) and samtools. The + /// `` prefix is the value written to the `SO` tag for this + /// order (`queryname` for queryname sub-sorts; `unsorted` for + /// template-coordinate). Without the prefix, fgbio's `SamOrder.apply` + /// (which strips everything up to the first `:`) fails to re-recognize the + /// header. #[must_use] pub fn header_ss_tag(&self) -> Option<&'static str> { match self { Self::Coordinate => None, - Self::Queryname(QuerynameComparator::Lexicographic) => Some("lexicographic"), - Self::Queryname(QuerynameComparator::Natural) => Some("natural"), - Self::TemplateCoordinate => Some("template-coordinate"), + Self::Queryname(QuerynameComparator::Lexicographic) => Some("queryname:lexicographic"), + Self::Queryname(QuerynameComparator::Natural) => Some("queryname:natural"), + Self::TemplateCoordinate => Some("unsorted:template-coordinate"), } } } @@ -1489,16 +1497,20 @@ mod tests { #[test] fn test_sort_order_header_ss_tag() { + // The SS tag carries the `:` form (fgbio/samtools). assert_eq!(SortOrder::Coordinate.header_ss_tag(), None); assert_eq!( SortOrder::Queryname(QuerynameComparator::Lexicographic).header_ss_tag(), - Some("lexicographic") + Some("queryname:lexicographic") ); assert_eq!( SortOrder::Queryname(QuerynameComparator::Natural).header_ss_tag(), - Some("natural") + Some("queryname:natural") + ); + assert_eq!( + SortOrder::TemplateCoordinate.header_ss_tag(), + Some("unsorted:template-coordinate") ); - assert_eq!(SortOrder::TemplateCoordinate.header_ss_tag(), Some("template-coordinate")); } #[test] diff --git a/crates/fgumi-sort/src/lib.rs b/crates/fgumi-sort/src/lib.rs index 3d392f0ab..e36a44d95 100644 --- a/crates/fgumi-sort/src/lib.rs +++ b/crates/fgumi-sort/src/lib.rs @@ -137,8 +137,9 @@ pub(crate) fn create_output_header(sort_order: keys::SortOrder, header: &Header) keys::SortOrder::TemplateCoordinate => { hd.other_fields_mut().insert(header_tag::SORT_ORDER, BString::from("unsorted")); hd.other_fields_mut().insert(header_tag::GROUP_ORDER, BString::from("query")); - hd.other_fields_mut() - .insert(header_tag::SUBSORT_ORDER, BString::from("template-coordinate")); + if let Some(ss) = sort_order.header_ss_tag() { + hd.other_fields_mut().insert(header_tag::SUBSORT_ORDER, BString::from(ss)); + } } } @@ -247,4 +248,36 @@ mod tests { let so = hd.other_fields().get(b"SO").expect("should have SO"); assert_eq!(<_ as AsRef<[u8]>>::as_ref(so), b"coordinate"); } + + #[test] + fn test_create_output_header_ss_has_sort_order_prefix() { + // fgbio/samtools write SS as `:`; pin that fgumi + // does too, so fgbio's `SamOrder.apply` re-recognizes the header. + let empty = Header::builder().build(); + let ss_of = |hd: &Map| { + <_ as AsRef<[u8]>>::as_ref(hd.other_fields().get(b"SS").expect("SS")).to_vec() + }; + + // queryname (natural) -> SO:queryname, SS:queryname:natural + let out = create_output_header( + keys::SortOrder::Queryname(keys::QuerynameComparator::Natural), + &empty, + ); + let hd = out.header().expect("@HD"); + assert_eq!( + <_ as AsRef<[u8]>>::as_ref(hd.other_fields().get(b"SO").expect("SO")), + b"queryname" + ); + assert_eq!(ss_of(hd), b"queryname:natural"); + + // template-coordinate -> SO:unsorted, GO:query, SS:unsorted:template-coordinate + let out = create_output_header(keys::SortOrder::TemplateCoordinate, &empty); + let hd = out.header().expect("@HD"); + assert_eq!( + <_ as AsRef<[u8]>>::as_ref(hd.other_fields().get(b"SO").expect("SO")), + b"unsorted" + ); + assert_eq!(<_ as AsRef<[u8]>>::as_ref(hd.other_fields().get(b"GO").expect("GO")), b"query"); + assert_eq!(ss_of(hd), b"unsorted:template-coordinate"); + } } diff --git a/src/lib/commands/sort.rs b/src/lib/commands/sort.rs index 07ecf251d..99d1a718a 100644 --- a/src/lib/commands/sort.rs +++ b/src/lib/commands/sort.rs @@ -1004,14 +1004,14 @@ mod tests { fn test_queryname_lex_header_has_subsort() { let order = SortOrder::from(SortOrderArg::Queryname); assert_eq!(order.header_so_tag(), "queryname"); - assert_eq!(order.header_ss_tag(), Some("lexicographic")); + assert_eq!(order.header_ss_tag(), Some("queryname:lexicographic")); } #[test] fn test_queryname_natural_header_has_subsort() { let order = SortOrder::from(SortOrderArg::QuerynameNatural); assert_eq!(order.header_so_tag(), "queryname"); - assert_eq!(order.header_ss_tag(), Some("natural")); + assert_eq!(order.header_ss_tag(), Some("queryname:natural")); } #[test] @@ -1025,7 +1025,7 @@ mod tests { fn test_template_coordinate_header_subsort() { let order = SortOrder::from(SortOrderArg::TemplateCoordinate); assert_eq!(order.header_so_tag(), "unsorted"); - assert_eq!(order.header_ss_tag(), Some("template-coordinate")); + assert_eq!(order.header_ss_tag(), Some("unsorted:template-coordinate")); } #[test] From 9676dc1e5091c74f78a80da5338a64329f290a56 Mon Sep 17 00:00:00 2001 From: Nils Homer Date: Thu, 9 Jul 2026 21:18:10 -0700 Subject: [PATCH 2/2] fix(simulate): write SS sub-sort tag via SortOrder accessors (R2-HDR-01) simulate grouped-reads / mapped-reads hand-formatted the `SS` sub-sort tag as a bare literal. Reuse `SortOrder::TemplateCoordinate.header_{so,go,ss}_tag()` so the emitted `@HD SS` matches the : form this PR establishes for `sort`, keeping every SS-writing site on one canonical mapping. Stacked on #514 (depends on its keys.rs prefix). --- src/lib/commands/simulate/common.rs | 25 ++++++++++++++++++++++ src/lib/commands/simulate/grouped_reads.rs | 20 ++++++----------- src/lib/commands/simulate/mapped_reads.rs | 20 ++++++----------- 3 files changed, 37 insertions(+), 28 deletions(-) diff --git a/src/lib/commands/simulate/common.rs b/src/lib/commands/simulate/common.rs index f9f4553d9..b06263132 100644 --- a/src/lib/commands/simulate/common.rs +++ b/src/lib/commands/simulate/common.rs @@ -6,13 +6,38 @@ use anyhow::{Context, Result, bail}; use clap::Args; use fgumi_consensus::MethylationMode; use fgumi_consensus::methylation::is_cpg_context; +use fgumi_sort::SortOrder; use log::info; use noodles::fasta; +use noodles::sam::header::record::value::Map; +use noodles::sam::header::record::value::map::header::{self as HeaderRecord, Tag as HeaderTag}; use rand::{Rng, RngExt}; use std::fs::File; use std::io::BufReader; use std::path::PathBuf; +/// Build the `@HD` header map carrying the `SO`/`GO`/`SS` sort-order tags for `order`. +/// +/// Shared by the `simulate` subcommands so they emit the same `SO`/`GO`/`SS` values — +/// including the `SS = :` prefix — that `fgumi sort` writes (see +/// [`SortOrder::header_ss_tag`]). `GO` and `SS` are inserted only when the order defines +/// them, so this is safe for any [`SortOrder`], not just template-coordinate. +#[must_use] +pub fn build_sort_order_header_map(order: SortOrder) -> Map { + let HeaderTag::Other(so_tag) = HeaderTag::from([b'S', b'O']) else { unreachable!() }; + let HeaderTag::Other(go_tag) = HeaderTag::from([b'G', b'O']) else { unreachable!() }; + let HeaderTag::Other(ss_tag) = HeaderTag::from([b'S', b'S']) else { unreachable!() }; + + let mut builder = Map::::builder().insert(so_tag, order.header_so_tag()); + if let Some(go) = order.header_go_tag() { + builder = builder.insert(go_tag, go); + } + if let Some(ss) = order.header_ss_tag() { + builder = builder.insert(ss_tag, ss); + } + builder.build().expect("header map with valid SO/GO/SS tags") +} + /// Common simulation options shared across all simulate subcommands. #[derive(Args, Debug, Clone)] pub struct SimulationCommon { diff --git a/src/lib/commands/simulate/grouped_reads.rs b/src/lib/commands/simulate/grouped_reads.rs index 921599abb..a26a998e8 100644 --- a/src/lib/commands/simulate/grouped_reads.rs +++ b/src/lib/commands/simulate/grouped_reads.rs @@ -10,7 +10,8 @@ use crate::commands::common::{CompressionOptions, parse_bool}; use crate::commands::simulate::common::{ FamilySizeArgs, InsertSizeArgs, MethylationArgs, MethylationConfig, MoleculeInfo, PositionDistArgs, QualityArgs, ReferenceArgs, ReferenceGenome, SimulationCommon, - StrandBiasArgs, apply_methylation_conversion, generate_random_sequence, pad_sequence, + StrandBiasArgs, apply_methylation_conversion, build_sort_order_header_map, + generate_random_sequence, pad_sequence, }; use crate::commands::simulate::region_to_bin; use crate::dna::reverse_complement; @@ -24,9 +25,9 @@ use clap::Parser; use fgumi_bam_io::ProgressTracker; use fgumi_bam_io::create_raw_bam_writer; use fgumi_raw_bam::{RawRecord, SamBuilder, flags as raw_flags}; +use fgumi_sort::SortOrder; use log::info; use noodles::sam::header::Header; -use noodles::sam::header::record::value::map::header::{self as HeaderRecord, Tag as HeaderTag}; use rand::{Rng, RngExt}; use std::fs::File; use std::io::{BufWriter, Write}; @@ -146,18 +147,9 @@ impl Command for GroupedReads { let ref_header = ref_genome.build_bam_header(); let mut header_builder = Header::builder(); - // Add sort order tags: SO:unsorted, GO:query, SS:template-coordinate - let HeaderTag::Other(so_tag) = HeaderTag::from([b'S', b'O']) else { unreachable!() }; - let HeaderTag::Other(go_tag) = HeaderTag::from([b'G', b'O']) else { unreachable!() }; - let HeaderTag::Other(ss_tag) = HeaderTag::from([b'S', b'S']) else { unreachable!() }; - - let header_map = - noodles::sam::header::record::value::Map::::builder() - .insert(so_tag, "unsorted") - .insert(go_tag, "query") - .insert(ss_tag, "template-coordinate") - .build() - .expect("header map with valid SO/GO/SS tags"); + // Add sort order tags (SO/GO/SS) matching template-coordinate order; see + // `SortOrder::header_ss_tag()` for why the `SS` value is prefixed with the `SO` value. + let header_map = build_sort_order_header_map(SortOrder::TemplateCoordinate); header_builder = header_builder.set_header(header_map); for (name, map) in ref_header.reference_sequences() { diff --git a/src/lib/commands/simulate/mapped_reads.rs b/src/lib/commands/simulate/mapped_reads.rs index 2bd98f36f..0fb552eec 100644 --- a/src/lib/commands/simulate/mapped_reads.rs +++ b/src/lib/commands/simulate/mapped_reads.rs @@ -10,7 +10,8 @@ use crate::commands::common::CompressionOptions; use crate::commands::simulate::common::{ FamilySizeArgs, InsertSizeArgs, MethylationArgs, MethylationConfig, MoleculeInfo, PositionDistArgs, QualityArgs, ReferenceArgs, ReferenceGenome, SimulationCommon, - apply_methylation_conversion, generate_random_sequence, pad_sequence, validate_rate, + apply_methylation_conversion, build_sort_order_header_map, generate_random_sequence, + pad_sequence, validate_rate, }; use crate::commands::simulate::region_to_bin; use crate::dna::reverse_complement; @@ -23,9 +24,9 @@ use clap::Parser; use fgumi_bam_io::ProgressTracker; use fgumi_bam_io::create_raw_bam_writer; use fgumi_raw_bam::{RawRecord, RawRecordView, SamBuilder, flags as raw_flags}; +use fgumi_sort::SortOrder; use log::info; use noodles::sam::header::Header; -use noodles::sam::header::record::value::map::header::{self as HeaderRecord, Tag as HeaderTag}; use rand::{Rng, RngExt}; use std::fs::File; use std::io::{BufWriter, Write}; @@ -148,18 +149,9 @@ impl Command for MappedReads { let ref_header = ref_genome.build_bam_header(); let mut header_builder = Header::builder(); - // Add sort order tags: SO:unsorted, GO:query, SS:template-coordinate - let HeaderTag::Other(so_tag) = HeaderTag::from([b'S', b'O']) else { unreachable!() }; - let HeaderTag::Other(go_tag) = HeaderTag::from([b'G', b'O']) else { unreachable!() }; - let HeaderTag::Other(ss_tag) = HeaderTag::from([b'S', b'S']) else { unreachable!() }; - - let header_map = - noodles::sam::header::record::value::Map::::builder() - .insert(so_tag, "unsorted") - .insert(go_tag, "query") - .insert(ss_tag, "template-coordinate") - .build() - .expect("header map with valid SO/GO/SS tags"); + // Add sort order tags (SO/GO/SS) matching template-coordinate order; see + // `SortOrder::header_ss_tag()` for why the `SS` value is prefixed with the `SO` value. + let header_map = build_sort_order_header_map(SortOrder::TemplateCoordinate); header_builder = header_builder.set_header(header_map); // Re-add reference sequences from ref_genome header