refactor(commands): deduplicate simulate utilities and clip header logic

src/commands/clip.rs

@@ -287,46 +287,10 @@ impl Clip {
 
         // Open input BAM with MT BGZF decompression
         let reader_threads = self.threading.num_threads();
-        let (mut reader, mut header) = create_bam_reader(&self.io.input, reader_threads)?;
+        let (mut reader, header) = create_bam_reader(&self.io.input, reader_threads)?;
 
         // Update header sort order if specified
-        if let Some(ref sort_order) = self.sort_order {
-            use bstr::BString;
-            use noodles::sam::header::record::value::Map;
-            use noodles::sam::header::record::value::map::header::tag::SORT_ORDER;
-
-            // Get or create the header map
-            let mut header_map = if let Some(hd) = header.header() {
-                hd.clone()
-            } else {
-                Map::<noodles::sam::header::record::value::map::Header>::default()
-            };
-
-            // Update sort order
-            *header_map.other_fields_mut().entry(SORT_ORDER).or_insert(BString::from("")) =
-                BString::from(sort_order.as_str());
-
-            // Rebuild header with new header map
-            let mut builder = noodles::sam::Header::builder();
-
-            // Copy existing components
-            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());
-            }
-
-            // Set the modified header
-            builder = builder.set_header(header_map);
-            header = builder.build();
-        }
+        let header = self.update_header_sort_order(header)?;
 
         // Add @PG record with PP chaining to input's last program
         let header = crate::commands::common::add_pg_record(header, command_line)?;

src/commands/simulate/common.rs

@@ -1,6 +1,8 @@
-//! Shared CLI arguments for simulation commands.
+//! Shared CLI arguments and utilities for simulation commands.
 
+use super::sort::TemplateCoordKey;
 use clap::Args;
+use rand::{Rng, RngExt};
 use std::path::PathBuf;
 
 /// Common simulation options shared across all simulate subcommands.
@@ -205,6 +207,67 @@ pub struct PositionDistArgs {
     pub umis_per_position: usize,
 }
 
+/// Generate a random DNA sequence of the given length.
+pub(super) fn generate_random_sequence(len: usize, rng: &mut impl Rng) -> Vec<u8> {
+    const BASES: &[u8] = b"ACGT";
+    let mut seq = Vec::with_capacity(len);
+    for _ in 0..len {
+        seq.push(BASES[rng.random_range(0..4)]);
+    }
+    seq
+}
+
+/// Pad a sequence to the target length with random bases, or truncate if too long.
+pub(super) fn pad_sequence(mut seq: Vec<u8>, target_len: usize, rng: &mut impl Rng) -> Vec<u8> {
+    while seq.len() < target_len {
+        seq.push(b"ACGT"[rng.random_range(0..4)]);
+    }
+    seq.truncate(target_len);
+    seq
+}
+
+/// Compute the genomic position for a molecule based on its ID.
+#[inline]
+pub(super) fn compute_position(mol_id: usize, num_positions: usize, ref_length: usize) -> usize {
+    if num_positions == 0 {
+        return 100;
+    }
+    let usable_span = ref_length.saturating_sub(1000);
+    if usable_span == 0 {
+        return 100;
+    }
+    let position_idx = mol_id % num_positions;
+    ((position_idx as f64 / num_positions as f64) * usable_span as f64) as usize + 100
+}
+
+/// Lightweight molecule info for position-first sorting.
+#[derive(Debug)]
+pub(super) struct MoleculeInfo {
+    pub mol_id: usize,
+    pub seed: u64,
+    pub sort_key: TemplateCoordKey,
+}
+
+impl Ord for MoleculeInfo {
+    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
+        self.sort_key.cmp(&other.sort_key)
+    }
+}
+
+impl PartialOrd for MoleculeInfo {
+    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
+        Some(self.cmp(other))
+    }
+}
+
+impl PartialEq for MoleculeInfo {
+    fn eq(&self, other: &Self) -> bool {
+        self.sort_key == other.sort_key
+    }
+}
+
+impl Eq for MoleculeInfo {}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -585,4 +648,67 @@ mod tests {
         // Should not panic
         let _ = args_lower.to_family_size_distribution().unwrap();
     }
+
+    #[test]
+    fn test_compute_position_zero_num_positions() {
+        // Should not panic on num_positions == 0
+        let pos = compute_position(5, 0, 250_000_000);
+        assert_eq!(pos, 100);
+    }
+
+    #[test]
+    fn test_compute_position_small_ref_length() {
+        // Should not underflow when ref_length < 1000
+        let pos = compute_position(0, 10, 500);
+        assert_eq!(pos, 100);
+    }
+
+    #[test]
+    fn test_compute_position_normal() {
+        // First position should be at offset 100
+        assert_eq!(compute_position(0, 10, 10_000), 100);
+        // Last position should be near ref_length - 900
+        let pos = compute_position(9, 10, 10_000);
+        assert!(pos > 100);
+        assert!(pos < 10_000);
+    }
+
+    fn make_molecule_info(mol_id: usize, tid1: i32, pos1: i64) -> MoleculeInfo {
+        MoleculeInfo {
+            mol_id,
+            seed: 0,
+            sort_key: TemplateCoordKey {
+                tid1,
+                tid2: 0,
+                pos1,
+                pos2: 0,
+                neg1: false,
+                neg2: false,
+                mid: String::new(),
+                name: String::new(),
+                is_upper_of_pair: false,
+            },
+        }
+    }
+
+    #[test]
+    fn test_molecule_info_ordering() {
+        let a = make_molecule_info(0, 1, 100);
+        let b = make_molecule_info(1, 1, 200);
+        let c = make_molecule_info(2, 2, 50);
+
+        // Same tid, earlier position comes first
+        assert!(a < b);
+        assert!(b > a);
+        // Higher tid comes later regardless of position
+        assert!(b < c);
+
+        // PartialOrd consistent with Ord
+        assert_eq!(a.partial_cmp(&b), Some(std::cmp::Ordering::Less));
+
+        // Equal sort keys are equal
+        let a2 = make_molecule_info(99, 1, 100);
+        assert_eq!(a, a2);
+        assert_eq!(a.cmp(&a2), std::cmp::Ordering::Equal);
+    }
 }

src/commands/simulate/consensus_reads.rs

@@ -2,11 +2,12 @@
 
 use crate::commands::command::Command;
 use crate::commands::common::CompressionOptions;
-use crate::commands::simulate::common::StrandBiasArgs;
+use crate::commands::simulate::common::{StrandBiasArgs, generate_random_sequence};
 use anyhow::{Context, Result};
 use clap::Parser;
 use crossbeam_channel::bounded;
 use fgumi_lib::bam_io::create_bam_writer;
+use fgumi_lib::dna::reverse_complement;
 use fgumi_lib::progress::ProgressTracker;
 use fgumi_lib::sam::builder::{ConsensusTagsBuilder, RecordBuilder};
 use fgumi_lib::simulate::{StrandBiasModel, create_rng};
@@ -395,29 +396,6 @@ fn build_consensus_record(
         .build()
 }
 
-fn generate_random_sequence(len: usize, rng: &mut impl Rng) -> Vec<u8> {
-    const BASES: &[u8] = b"ACGT";
-    let mut seq = Vec::with_capacity(len);
-    for _ in 0..len {
-        seq.push(BASES[rng.random_range(0..4)]);
-    }
-    seq
-}
-
-fn reverse_complement(seq: &[u8]) -> Vec<u8> {
-    let mut result = Vec::with_capacity(seq.len());
-    for &b in seq.iter().rev() {
-        result.push(match b {
-            b'A' => b'T',
-            b'T' => b'A',
-            b'C' => b'G',
-            b'G' => b'C',
-            _ => b'N',
-        });
-    }
-    result
-}
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -708,7 +686,7 @@ mod tests {
     #[test]
     fn test_reverse_complement_unknown_bases() {
         assert_eq!(reverse_complement(b"N"), b"N");
-        assert_eq!(reverse_complement(b"X"), b"N");
+        assert_eq!(reverse_complement(b"X"), b"X");
         assert_eq!(reverse_complement(b"ANCG"), b"CGNT");
     }
 

src/commands/simulate/correct_reads.rs

@@ -2,6 +2,7 @@
 
 use crate::commands::command::Command;
 use crate::commands::common::CompressionOptions;
+use crate::commands::simulate::common::generate_random_sequence;
 use anyhow::{Context, Result};
 use clap::Parser;
 use crossbeam_channel::bounded;
@@ -397,15 +398,6 @@ fn generate_correct_read_pair(
     }
 }
 
-fn generate_random_sequence(len: usize, rng: &mut impl Rng) -> Vec<u8> {
-    const BASES: &[u8] = b"ACGT";
-    let mut seq = Vec::with_capacity(len);
-    for _ in 0..len {
-        seq.push(BASES[rng.random_range(0..4)]);
-    }
-    seq
-}
-
 fn introduce_n_errors(umi: &str, n: usize, rng: &mut impl Rng) -> String {
     const BASES: &[u8] = b"ACGT";
     let mut result: Vec<u8> = umi.as_bytes().to_vec();

src/commands/simulate/fastq_reads.rs

@@ -2,11 +2,12 @@
 
 use crate::commands::command::Command;
 use crate::commands::simulate::common::{
-    FamilySizeArgs, InsertSizeArgs, QualityArgs, SimulationCommon,
+    FamilySizeArgs, InsertSizeArgs, QualityArgs, SimulationCommon, generate_random_sequence,
 };
 use anyhow::{Context, Result, bail};
 use clap::Parser;
 use crossbeam_channel::bounded;
+use fgumi_dna::dna::complement_base;
 use fgumi_lib::progress::ProgressTracker;
 use fgumi_lib::simulate::{FastqWriter, create_rng};
 use log::info;
@@ -532,27 +533,11 @@ fn generate_molecule_reads(
     records
 }
 
-/// Generate a random DNA sequence.
-fn generate_random_sequence(len: usize, rng: &mut impl Rng) -> Vec<u8> {
-    const BASES: &[u8] = b"ACGT";
-    let mut seq = Vec::with_capacity(len);
-    for _ in 0..len {
-        seq.push(BASES[rng.random_range(0..4)]);
-    }
-    seq
-}
-
 /// Reverse complement a DNA sequence into an existing buffer.
 fn reverse_complement_into(seq: &[u8], output: &mut Vec<u8>) {
     output.reserve(seq.len());
     for &b in seq.iter().rev() {
-        output.push(match b {
-            b'A' => b'T',
-            b'T' => b'A',
-            b'C' => b'G',
-            b'G' => b'C',
-            _ => b'N',
-        });
+        output.push(complement_base(b));
     }
 }
 
@@ -698,9 +683,9 @@ mod tests {
 
     #[test]
     fn test_reverse_complement_unknown_bases() {
-        // Unknown bases should become N
+        // Unknown bases are returned unchanged by the library complement_base
         assert_eq!(reverse_complement(b"N"), b"N");
-        assert_eq!(reverse_complement(b"X"), b"N");
+        assert_eq!(reverse_complement(b"X"), b"X");
         assert_eq!(reverse_complement(b"ANT"), b"ANT"); // A->T, N->N, T->A, reversed = ANT
     }