subsample: Fix skip-ahead bias, reproducibility, bgz, and error handling

- Skip-ahead no longer always starts at record 0; the first jump is
  drawn before selecting any record
- Tile processing order is deterministic: tiles sorted by bin_key,
  per-tile RNG derived from base_seed + bin_key (no shared mutex)
- Remove .bgz from is_gzipped() to match fastq::fs which only
  recognizes .gz
- Temp file creation errors propagate as SubsampleError instead of
  panicking via unwrap()

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: jdidion

src/commands/subsample.rs

@@ -4,7 +4,7 @@ use std::{
     io::{self, BufRead, BufReader, BufWriter, Write},
     ops::{Bound, RangeBounds},
     path::{Path, PathBuf},
-    sync::{Mutex, mpsc},
+    sync::mpsc,
 };
 
 use bitvec::vec::BitVec;
@@ -85,12 +85,7 @@ pub fn subsample(args: SubsampleArgs) -> Result<(), SubsampleError> {
 }
 
 fn is_gzipped(path: &Path) -> bool {
-    path
-        .extension()
-        .and_then(|e| e.to_str())
-        .map(|ext| ["gz", "bgz"]
-        .contains(&ext))
-        .unwrap_or(false)
+    path.extension().and_then(|e| e.to_str()) == Some("gz")
 }
 
 fn subsample_approximate<Rng>(
@@ -563,6 +558,7 @@ impl Write for OutputWriter {
 }
 
 struct TileInfo {
+    bin_key: u64,
     record_count: usize,
     r1_path: PathBuf,
     r2_path: Option<PathBuf>,
@@ -609,24 +605,31 @@ fn write_tile_temp_files(
         }
 
         if let Some(bin_key) = parse_tile_bin(r1_rec.name()) {
-            let tw = tile_writers.entry(bin_key).or_insert_with(|| {
+            if !tile_writers.contains_key(&bin_key) {
                 let r1_path = temp_dir.join(format!("{bin_key}.r1.fq"));
-                let r1_file = BufWriter::new(File::create(&r1_path).unwrap());
+                let r1_file = BufWriter::new(
+                    File::create(&r1_path)
+                        .map_err(|e| SubsampleError::CreateFile(e, r1_path.clone()))?,
+                );
                 let (r2_writer, r2_path) = if paired {
                     let p = temp_dir.join(format!("{bin_key}.r2.fq"));
-                    let f = BufWriter::new(File::create(&p).unwrap());
+                    let f = BufWriter::new(
+                        File::create(&p)
+                            .map_err(|e| SubsampleError::CreateFile(e, p.clone()))?,
+                    );
                     (Some(fastq::io::Writer::new(f)), Some(p))
                 } else {
                     (None, None)
                 };
-                TileWriter {
+                tile_writers.insert(bin_key, TileWriter {
                     r1: fastq::io::Writer::new(r1_file),
                     r1_path,
                     r2: r2_writer,
                     r2_path,
                     record_count: 0,
-                }
-            });
+                });
+            }
+            let tw = tile_writers.get_mut(&bin_key).unwrap();
 
             tw.r1.write_record(&r1_rec)?;
             if let Some(r2w) = tw.r2.as_mut() {
@@ -639,22 +642,25 @@ fn write_tile_temp_files(
     }
 
     // Collect TileInfo (writers are dropped here, flushing buffers)
-    let tiles: Vec<TileInfo> = tile_writers
-        .into_values()
-        .map(|tw| TileInfo {
+    let mut tiles: Vec<TileInfo> = tile_writers
+        .into_iter()
+        .map(|(key, tw)| TileInfo {
+            bin_key: key,
             record_count: tw.record_count,
             r1_path: tw.r1_path,
             r2_path: tw.r2_path,
         })
         .collect();
+    // Sort by bin key for deterministic processing order
+    tiles.sort_unstable_by_key(|t| t.bin_key);
 
     Ok((tiles, parse_failures))
 }
 
 fn subsample_by_tile<Rng>(
     (r1_src, r1_dst): (&Path, &Path),
     (r2_src, r2_dst): (Option<&Path>, Option<&Path>),
-    rng: Rng,
+    mut rng: Rng,
     mode: TileCountMode,
     fast: bool,
     sampling_threads: usize,
@@ -775,8 +781,12 @@ where
         retained_count
     );
 
+    // Derive a base seed from the caller's RNG once, then derive per-tile
+    // seeds deterministically from base_seed + bin_key. This avoids a shared
+    // mutex and makes output reproducible regardless of thread scheduling.
+    let base_seed: u64 = rng.random();
+
     let (tx, rx) = mpsc::sync_channel::<TileResult>(threads * 2);
-    let rng = Mutex::new(rng);
 
     let selected_records = std::thread::scope(|s| -> Result<usize, SubsampleError> {
         // Writer thread: receives tile results and writes to output files
@@ -813,14 +823,11 @@ where
             .into_iter()
             .map(|chunk| {
                 let tx = tx.clone();
-                let rng = &rng;
 
                 s.spawn(move || -> Result<(), SubsampleError> {
                     for tile in chunk {
-                        let mut tile_rng = {
-                            let mut rng = rng.lock().unwrap();
-                            SmallRng::seed_from_u64(rng.random())
-                        };
+                        let mut tile_rng =
+                            SmallRng::seed_from_u64(base_seed.wrapping_add(tile.bin_key));
 
                         // Skip-ahead only works for single-end; paired-end
                         // always uses exact indexed sampling.

src/fastq/io/indexed_reader.rs

@@ -117,9 +117,12 @@ impl<R: Read + Seek> IndexedReader<R> {
 
         let mean_jump = total as f64 / target as f64;
         let mut record = Record::default();
-        let mut current: usize = 0;
         let mut count = 0;
 
+        // Draw the first jump to avoid always starting at record 0
+        let first = (-mean_jump * rng.random::<f64>().ln()).ceil() as usize;
+        let mut current: usize = first.min(total.saturating_sub(1));
+
         while count < target && current < total {
             self.read_record_at(current, &mut record)?;
             callback(&record)?;