fix(simulate): write SS sub-sort tag via SortOrder accessors (R2-HDR-01)

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");
     }
 
     // ========================================================================

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 `<sort-order>:<sub-sort>`, matching fgbio
+    /// (`SamOrder.applyTo`: `sortOrder.name() + ":" + ss`) and samtools. The
+    /// `<sort-order>` 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 `<sort-order>:<sub-sort>` 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]

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 `<sort-order>:<sub-sort>`; pin that fgumi
+        // does too, so fgbio's `SamOrder.apply` re-recognizes the header.
+        let empty = Header::builder().build();
+        let ss_of = |hd: &Map<noodles::sam::header::record::value::map::Header>| {
+            <_ 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");
+    }
 }

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 = <sort-order>:<sub-sort>` 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<HeaderRecord::Header> {
+    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::<HeaderRecord::Header>::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 {

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::<HeaderRecord::Header>::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() {

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::<HeaderRecord::Header>::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

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]