Make indexed sampling and threshold caps deterministic

modkit-core/src/command_utils.rs

@@ -16,6 +16,23 @@ use crate::threshold_mod_caller::MultipleThresholdModCaller;
 use crate::thresholds::calc_threshold_from_bam;
 use crate::util::{create_out_directory, Region};
 
+pub(crate) fn parse_sampling_fraction(raw: &str) -> Result<f64, String> {
+    let fraction = raw.parse::<f64>().map_err(|_| {
+        format!(
+            "sampling fraction must be a finite number in the inclusive \
+             range [0, 1]; got '{raw}'"
+        )
+    })?;
+    if fraction.is_finite() && (0.0..=1.0).contains(&fraction) {
+        Ok(fraction)
+    } else {
+        Err(format!(
+            "sampling fraction must be a finite number in the inclusive \
+             range [0, 1]; got '{raw}'"
+        ))
+    }
+}
+
 pub fn parse_per_mod_thresholds(
     raw_per_mod_thresholds: &[String],
 ) -> anyhow::Result<HashMap<ModCodeRepr, f32>> {
@@ -454,3 +471,41 @@ pub(crate) fn parse_raw_motifs(
         e @ _ => e,
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::parse_sampling_fraction;
+
+    #[test]
+    fn parse_sampling_fraction_matrix() {
+        for (raw, expected) in [
+            ("0", 0.0),
+            ("-0.0", -0.0),
+            ("0.25", 0.25),
+            ("1e-3", 0.001),
+            ("1", 1.0),
+            ("1e0", 1.0),
+        ] {
+            assert_eq!(parse_sampling_fraction(raw).unwrap(), expected);
+        }
+        assert!(parse_sampling_fraction("-0.0").unwrap().is_sign_negative());
+
+        for raw in [
+            "-0.0000001",
+            "1.0000001",
+            "NaN",
+            "nan",
+            "inf",
+            "+inf",
+            "-inf",
+            "1e309",
+            "not-a-number",
+            "",
+        ] {
+            assert!(
+                parse_sampling_fraction(raw).is_err(),
+                "expected {raw:?} to be rejected"
+            );
+        }
+    }
+}

modkit-core/src/extract/subcommand.rs

@@ -17,7 +17,7 @@ use modkit_logging::{init_logging, init_logging_smart};
 use crate::command_utils::{
     get_serial_reader, parse_edge_filter_input, parse_per_base_thresholds,
     parse_per_mod_thresholds, parse_raw_thresholds_string_with_default,
-    parse_thresholds, using_stream,
+    parse_sampling_fraction, parse_thresholds, using_stream,
 };
 use crate::extract::args::InputArgs;
 use crate::extract::util::{
@@ -450,12 +450,14 @@ pub struct EntryExtractCalls {
     /// In practice, 10-100 thousand reads is sufficient to estimate the model
     /// output distribution and determine the filtering threshold. See
     /// filtering.md for details on filtering.
+    /// Must be a finite value in the inclusive range [0, 1].
     #[clap(help_heading = "Sampling Options")]
     #[arg(
         group = "sampling_options",
         short = 'f',
         long,
-        hide_short_help = true
+        hide_short_help = true,
+        value_parser = parse_sampling_fraction
     )]
     sampling_frac: Option<f64>,
     /// Sample this many reads when estimating the filtering threshold. If a

modkit-core/src/modbam_util/subcommands.rs

@@ -22,7 +22,8 @@ use crate::command_utils::{
     get_bam_writer, get_motif_lookup_from_parts, get_serial_reader,
     get_threshold_from_options, parse_edge_filter_input, parse_forward_motifs,
     parse_per_mod_thresholds, parse_raw_motifs,
-    parse_raw_thresholds_string_with_default, parse_thresholds, using_stream,
+    parse_raw_thresholds_string_with_default, parse_sampling_fraction,
+    parse_thresholds, using_stream,
 };
 use crate::errs::{MkError, MkResult};
 use crate::interval_chunks::{
@@ -39,6 +40,7 @@ use crate::modbam_util::check_tags::ModTagViews;
 use crate::monoid::Moniod;
 use crate::motifs::motif_bed::RegexMotif;
 use crate::position_filter::StrandedPositionFilter;
+use crate::reads_sampler::deterministic_sampler::resolve_master_seed;
 use crate::reads_sampler::record_sampler::RecordSampler;
 use crate::reads_sampler::sample_reads_from_interval;
 use crate::reads_sampler::sampling_schedule::SamplingSchedule;
@@ -1012,21 +1014,23 @@ pub struct SampleModBaseProbs {
     num_reads: Option<usize>,
     /// Instead of using a defined number of reads, specify a fraction of reads
     /// to sample, for example 0.1 will sample 1/10th of the reads.
+    /// Must be a finite value in the inclusive range [0, 1].
     #[clap(help_heading = "Sampling Options")]
     #[arg(
         short = 'f',
         long,
         alias = "sample-frac",
         group = "sampling_options",
-        conflicts_with = "no_sampling"
+        conflicts_with = "no_sampling",
+        value_parser = parse_sampling_fraction
     )]
     sampling_frac: Option<f64>,
     /// No sampling, use all of the reads to calculate the filter thresholds.
     #[clap(help_heading = "Sampling Options")]
     #[arg(long, alias = "no_filtering", default_value_t = false)]
     no_sampling: bool,
-    /// Random seed for deterministic running, the default is
-    /// non-deterministic, only used when no BAM index is provided.
+    /// Provide a seed to make fractional read sampling decisions repeatable
+    /// for indexed and unindexed inputs.
     #[clap(help_heading = "Sampling Options")]
     #[arg(short, conflicts_with = "no_sampling", long)]
     seed: Option<u64>,
@@ -1247,6 +1251,11 @@ impl SampleModBaseProbs {
                     0f64,
                 )
             };
+        let master_seed = if rng_sample {
+            resolve_master_seed(self.seed)
+        } else {
+            self.seed.unwrap_or_default()
+        };
 
         let mut workers: Vec<Box<dyn ExtractProbsWorker>> =
             Vec::with_capacity(self.threads);
@@ -1273,7 +1282,7 @@ impl SampleModBaseProbs {
                     "collecting base and modification histograms at aligned \
                      positions"
                 );
-                for i in 0..self.threads {
+                for _ in 0..self.threads {
                     let w = RegionMleProbs::<
                         AlignedBaseAndModArgmaxProbs,
                         ProbsExtractor,
@@ -1284,7 +1293,7 @@ impl SampleModBaseProbs {
                         motif_bases,
                         edge_filter.as_ref(),
                         &thread_pool,
-                        self.seed.unwrap_or(i as u64),
+                        master_seed,
                         rng_sample,
                         sample_frac,
                         chrom_to_counts.clone(),
@@ -1293,7 +1302,7 @@ impl SampleModBaseProbs {
                 }
             } else {
                 info!("collecting base-level histograms at aligned positions");
-                for i in 0..self.threads {
+                for _ in 0..self.threads {
                     let w = RegionMleProbs::<
                         AlignedBaseArgmaxProbs,
                         ProbsExtractor,
@@ -1304,7 +1313,7 @@ impl SampleModBaseProbs {
                         motif_bases,
                         edge_filter.as_ref(),
                         &thread_pool,
-                        self.seed.unwrap_or(i as u64),
+                        master_seed,
                         rng_sample,
                         sample_frac,
                         chrom_to_counts.clone(),
@@ -1318,7 +1327,7 @@ impl SampleModBaseProbs {
                     "collecting base and modification histograms, using all \
                      read positions"
                 );
-                for i in 0..self.threads {
+                for _ in 0..self.threads {
                     let w = RegionMleProbs::<
                         BaseAndModArgmaxProbs,
                         ProbsExtractor,
@@ -1329,7 +1338,7 @@ impl SampleModBaseProbs {
                         motif_bases,
                         edge_filter.as_ref(),
                         &thread_pool,
-                        self.seed.unwrap_or(i as u64),
+                        master_seed,
                         rng_sample,
                         sample_frac,
                         chrom_to_counts.clone(),
@@ -1341,7 +1350,7 @@ impl SampleModBaseProbs {
                     "collecting base-level histograms, using all read \
                      positions"
                 );
-                for i in 0..self.threads {
+                for _ in 0..self.threads {
                     let w =
                         RegionMleProbs::<BaseArgmaxProbs, ProbsExtractor>::new(
                             &bam_fp,
@@ -1350,7 +1359,7 @@ impl SampleModBaseProbs {
                             motif_bases,
                             edge_filter.as_ref(),
                             &thread_pool,
-                            self.seed.unwrap_or(i as u64),
+                            master_seed,
                             rng_sample,
                             sample_frac,
                             chrom_to_counts.clone(),
@@ -1638,8 +1647,14 @@ pub struct ModSummarize {
     /// Instead of using a defined number of reads, specify a fraction of reads
     /// to sample when estimating the filter threshold. For example 0.1 will
     /// sample 1/10th of the reads.
+    /// Must be a finite value in the inclusive range [0, 1].
     #[clap(help_heading = "Sampling Options")]
-    #[arg(group = "sampling_options", short = 'f', long)]
+    #[arg(
+        group = "sampling_options",
+        short = 'f',
+        long,
+        value_parser = parse_sampling_fraction
+    )]
     sampling_frac: Option<f64>,
     /// Sets a random seed for deterministic running (when using
     /// --sample-frac).
@@ -2257,15 +2272,17 @@ pub struct CallMods {
     /// filter-percentile. In practice, 50-100 thousand reads is sufficient
     /// to estimate the model output distribution and determine the
     /// filtering threshold. See filtering.md for details on filtering.
+    /// Must be a finite value in the inclusive range [0, 1].
     #[arg(
         group = "sampling_options",
         short = 'f',
         long,
-        hide_short_help = true
+        hide_short_help = true,
+        value_parser = parse_sampling_fraction
     )]
     sampling_frac: Option<f64>,
-    /// Set a random seed for deterministic running, the default is
-    /// non-deterministic, only used when no BAM index is provided.
+    /// Provide a seed to make fractional read sampling decisions repeatable
+    /// for indexed and unindexed inputs.
     #[arg(
         long,
         conflicts_with = "num_reads",

modkit-core/src/pileup/subcommand.rs

@@ -19,8 +19,8 @@ use crate::command_utils::{
     get_motif_lookup_from_parts, get_threshold_from_options,
     parse_edge_filter_input, parse_per_base_thresholds,
     parse_per_mod_thresholds, parse_raw_motifs,
-    parse_raw_thresholds_string_with_default, parse_thresholds,
-    parse_thresholds_values,
+    parse_raw_thresholds_string_with_default, parse_sampling_fraction,
+    parse_thresholds, parse_thresholds_values,
 };
 use crate::fasta::MotifLocationsLookup;
 use crate::interval_chunks::{
@@ -41,6 +41,7 @@ use crate::pileup::pileup_processor::{
 };
 use crate::pileup::{ModBasePileup2, PileupNumericOptions};
 use crate::position_filter::StrandedPositionFilter;
+use crate::reads_sampler::deterministic_sampler::resolve_master_seed;
 use crate::reads_sampler::sampling_schedule::IdxStats;
 use crate::sample_probs::{
     calculate_reads_per_contig, run_extract_probs_workers,
@@ -172,12 +173,14 @@ pub struct ModBamPileup {
     /// In practice, 10-100 thousand reads is sufficient to estimate the model
     /// output distribution and determine the filtering threshold. See
     /// filtering.md for details on filtering.
+    /// Must be a finite value in the inclusive range [0, 1].
     #[clap(help_heading = "Sampling Options")]
     #[arg(
         group = "sampling_options",
         short = 'f',
         long,
-        hide_short_help = true
+        hide_short_help = true,
+        value_parser = parse_sampling_fraction
     )]
     sampling_frac: Option<f64>,
     /// Set a random seed for deterministic running, the default is
@@ -776,6 +779,11 @@ impl ModBamPileup {
                 self.sampling_interval_size,
                 sampling_region,
             )?;
+        let master_seed = if rng_sample {
+            resolve_master_seed(self.seed)
+        } else {
+            self.seed.unwrap_or_default()
+        };
         let chrom_to_counts = chrom_to_counts.map(|x| Arc::new(x));
 
         if let Some(preset) = preset {
@@ -806,7 +814,7 @@ impl ModBamPileup {
                     *motif_bases
                 }
             };
-            for i in 0..n_workers {
+            for _ in 0..n_workers {
                 let worker = RegionMleProbs::<
                     AlignedBaseArgmaxProbs,
                     ProbsExtractor,
@@ -817,7 +825,7 @@ impl ModBamPileup {
                     motif_bases,
                     edge_filter,
                     thread_pool,
-                    self.seed.unwrap_or(i as u64),
+                    master_seed,
                     rng_sample,
                     sample_frac,
                     chrom_to_counts.clone(),
@@ -832,7 +840,7 @@ impl ModBamPileup {
                     motif_bases.iter().unique().join(",")
                 )
             });
-            for i in 0..n_workers {
+            for _ in 0..n_workers {
                 let worker = RegionMleProbs::<
                     AlignedBaseArgmaxProbs,
                     ProbsExtractor,
@@ -843,7 +851,7 @@ impl ModBamPileup {
                     motif_bases,
                     edge_filter,
                     thread_pool,
-                    self.seed.unwrap_or(i as u64),
+                    master_seed,
                     rng_sample,
                     sample_frac,
                     chrom_to_counts.clone(),
@@ -856,7 +864,7 @@ impl ModBamPileup {
                     "calculating threshold value with probabilites from reads",
                 )
             });
-            for i in 0..n_workers {
+            for _ in 0..n_workers {
                 let worker =
                     RegionMleProbs::<BaseArgmaxProbs, ProbsExtractor>::new(
                         &self.in_bam,
@@ -865,7 +873,7 @@ impl ModBamPileup {
                         [DnaBase::A; 4],
                         edge_filter,
                         thread_pool,
-                        self.seed.unwrap_or(i as u64),
+                        master_seed,
                         rng_sample,
                         sample_frac,
                         chrom_to_counts.clone(),
@@ -1787,12 +1795,14 @@ pub struct DuplexModBamPileup {
     /// filter-percentile. In practice, 50-100 thousand reads is sufficient
     /// to estimate the model output distribution and determine the
     /// filtering threshold. See filtering.md for details on filtering.
+    /// Must be a finite value in the inclusive range [0, 1].
     #[clap(help_heading = "Sampling Options")]
     #[arg(
         group = "sampling_options",
         short = 'f',
         long,
-        hide_short_help = true
+        hide_short_help = true,
+        value_parser = parse_sampling_fraction
     )]
     sampling_frac: Option<f64>,
     /// Set a random seed for deterministic running, the default is