Correct DMR site statistics and segmentation completeness

book/src/advanced_usage.md

@@ -2683,7 +2683,8 @@ Single-site Options:
           Prior distribution for estimating MAP-based p-value. Should be two
           arguments for alpha and beta (e.g. 1.0 1.0). See
           `dmr_scoring_details.md` for additional details on how the metric is
-          calculated
+          calculated. Alpha and beta must each be positive, and their sum must
+          be >= 1.0
 
       --delta <DELTA>
           Consider only effect sizes greater than this when calculating the

book/src/dmr_scoring_details.md

@@ -55,6 +55,11 @@ Where \\(X\\) is the observations (\\(N_{\text{mod}}\\) and \\(N_{\text{canonica
     P(p | X) = \text{Beta}(\alpha_0 + N_{\text{mod}}, \beta_0 + N_{\text{can}})
 \\]
 Where \\(\alpha_0\\) and \\(\beta_0\\) are the parameters for the prior distribution \\(\text{Beta}(\alpha_0, \beta_0)\\).
+
+A user-specified prior must have \\(\alpha_0 > 0\\), \\(\beta_0 > 0\\), and \\(\alpha_0 + \beta_0 \ge 1\\).
+This input contract differs from the finite-domain condition for the closed-form density at zero: the paired posterior parameters must independently satisfy \\(\alpha_1 + \alpha_2 > 1\\) and \\(\beta_1 + \beta_2 > 1\\).
+Therefore, a positive prior whose parameters sum to exactly one remains valid input, while equality at either posterior pair-sum boundary is rejected.
+
 The advantage to this model is that as you collect more coverage, the variance of the posterior gets smaller - you're more confident that the true value of \\(p\\) is near the empirical mean.
 But when you have low coverage, you keep the uncertainty around.
 
@@ -168,4 +173,3 @@ To provide another metric that is more robust to high counts, Modkit DMR will ou
 In addition to the statistic, the high and low bound of the 95% confidence interval are reported.
 A CI value (high or low) of zero indicates that there is little certainty about there being a difference between the two conditions.
 Generally speaking, filtering or sorting on the lower bound is a good test for finding important changes.
-

modkit-core/src/dmr/beta_diff.rs

@@ -33,6 +33,8 @@ impl Counts {
     pub fn new(n_mod: usize, coverage: usize) -> anyhow::Result<Self> {
         if n_mod > coverage {
             bail!("n_mod cannot be > coverage")
+        } else if coverage == 0 {
+            bail!("coverage must be > 0")
         } else {
             let frac_modified = n_mod as f64 / coverage as f64;
             Ok(Self { n_mod, coverage, frac_modified })
@@ -44,6 +46,7 @@ impl Counts {
     }
 
     fn resize(&self, max_coverage: usize) -> Self {
+        assert!(max_coverage > 0, "max coverage must be greater than zero");
         if self.coverage > max_coverage {
             let n_mod = (self.frac_modified * max_coverage as f64).round();
             let frac_modified = n_mod / max_coverage as f64;
@@ -138,12 +141,37 @@ impl PMapEstimator {
         prior: BetaParams,
         rope: f64,
         cap_coverages: bool,
-    ) -> Self {
+    ) -> anyhow::Result<Self> {
         let mut max_coverages = if cap_coverages {
             max_coverages
         } else {
-            [max_coverages[0] * a_num_reps, max_coverages[1] * b_num_reps]
+            [
+                max_coverages[0].checked_mul(a_num_reps).ok_or_else(|| {
+                    anyhow!(
+                        "maximum coverage overflow while scaling control \
+                         coverage {} by {} replicates",
+                        max_coverages[0],
+                        a_num_reps
+                    )
+                })?,
+                max_coverages[1].checked_mul(b_num_reps).ok_or_else(|| {
+                    anyhow!(
+                        "maximum coverage overflow while scaling experiment \
+                         coverage {} by {} replicates",
+                        max_coverages[1],
+                        b_num_reps
+                    )
+                })?,
+            ]
         };
+        if max_coverages.iter().any(|coverage| *coverage == 0) {
+            bail!(
+                "resolved maximum coverage must be greater than zero for both \
+                 conditions, got control {} and experiment {}",
+                max_coverages[0],
+                max_coverages[1]
+            )
+        }
         for x in max_coverages.iter_mut() {
             if *x > MAX_COV_ALLOWED {
                 info!(
@@ -154,7 +182,7 @@ impl PMapEstimator {
             }
         }
 
-        Self { max_coverages, prior, rope }
+        Ok(Self { max_coverages, prior, rope })
     }
 
     fn calc_posterior_params(&self, counts: &Counts) -> BetaParams {
@@ -174,8 +202,8 @@ impl PMapEstimator {
         let ln_A = ln_beta(params1.alpha, params1.beta)
             + ln_beta(params2.alpha, params2.beta);
         if d.abs() < self.rope {
-            if (params1.alpha + params2.alpha < 1f64)
-                || (params1.beta + params2.beta < 1f64)
+            if (params1.alpha + params2.alpha <= 1f64)
+                || (params1.beta + params2.beta <= 1f64)
             {
                 bail!(
                     "alpha1 + alpha2 <= 1 or beta1 + beta2 <= 1, params1 \
@@ -278,10 +306,139 @@ impl PMapEstimator {
 
 #[cfg(test)]
 mod tests {
-    use crate::dmr::beta_diff::{appell_f1_stable, LOWER, UPPER};
+    use crate::dmr::beta_diff::{
+        appell_f1_stable, BetaParams, PMapEstimator, LOWER, UPPER,
+    };
     use assert_approx_eq::assert_approx_eq;
     use rv::misc::gauss_legendre_quadrature;
 
+    fn estimator() -> PMapEstimator {
+        PMapEstimator::new(
+            [10, 10],
+            1,
+            1,
+            BetaParams::new(1.0, 1.0).unwrap(),
+            0.05,
+            true,
+        )
+        .unwrap()
+    }
+
+    fn make_estimator(
+        max_coverages: [usize; 2],
+        a_num_reps: usize,
+        b_num_reps: usize,
+        cap_coverages: bool,
+    ) -> anyhow::Result<PMapEstimator> {
+        PMapEstimator::new(
+            max_coverages,
+            a_num_reps,
+            b_num_reps,
+            BetaParams::new(1.0, 1.0).unwrap(),
+            0.05,
+            cap_coverages,
+        )
+    }
+
+    #[test]
+    fn estimator_rejects_zero_resolved_max_coverages() {
+        for max_coverages in [[0, 0], [0, 10], [10, 0]] {
+            let error = make_estimator(max_coverages, 1, 1, true)
+                .err()
+                .expect("zero maximum coverage must be rejected");
+            assert_eq!(
+                error.to_string(),
+                format!(
+                    "resolved maximum coverage must be greater than zero for \
+                     both conditions, got control {} and experiment {}",
+                    max_coverages[0], max_coverages[1]
+                )
+            );
+        }
+    }
+
+    #[test]
+    fn estimator_accepts_and_preserves_positive_max_coverages() {
+        for max_coverages in [[1, 1], [10, 10]] {
+            let estimator = make_estimator(max_coverages, 1, 1, true).unwrap();
+            assert_eq!(estimator.max_coverages, max_coverages);
+        }
+    }
+
+    #[test]
+    fn estimator_rejects_replicate_scaling_overflow() {
+        for (max_coverages, a_num_reps, b_num_reps, expected) in [
+            (
+                [usize::MAX, 10],
+                2,
+                1,
+                format!(
+                    "maximum coverage overflow while scaling control \
+                     coverage {} by 2 replicates",
+                    usize::MAX
+                ),
+            ),
+            (
+                [10, usize::MAX],
+                1,
+                2,
+                format!(
+                    "maximum coverage overflow while scaling experiment \
+                     coverage {} by 2 replicates",
+                    usize::MAX
+                ),
+            ),
+        ] {
+            let error =
+                make_estimator(max_coverages, a_num_reps, b_num_reps, false)
+                    .err()
+                    .expect("overflowing maximum coverage must be rejected");
+            assert_eq!(error.to_string(), expected);
+        }
+    }
+
+    #[test]
+    #[should_panic(expected = "max coverage must be greater than zero")]
+    fn counts_resize_requires_positive_max_coverage() {
+        super::Counts::new(1, 1).unwrap().resize(0);
+    }
+
+    #[test]
+    fn beta_diff_at_zero_rejects_exact_posterior_pair_sum_boundaries() {
+        let estimator = estimator();
+        let boundary_pairs = [
+            (
+                BetaParams::new(0.4, 2.0).unwrap(),
+                BetaParams::new(0.6, 3.0).unwrap(),
+            ),
+            (
+                BetaParams::new(2.0, 0.4).unwrap(),
+                BetaParams::new(3.0, 0.6).unwrap(),
+            ),
+        ];
+
+        for (params1, params2) in boundary_pairs {
+            let err = estimator
+                .calc_beta_diff(0.0, &params1, &params2)
+                .expect_err("pair-sum boundary must be rejected");
+            assert!(err
+                .to_string()
+                .starts_with("alpha1 + alpha2 <= 1 or beta1 + beta2 <= 1,"));
+        }
+    }
+
+    #[test]
+    fn beta_diff_at_zero_matches_interior_closed_form_value() {
+        let estimator = estimator();
+        let params1 = BetaParams::new(2.0, 3.0).unwrap();
+        let params2 = BetaParams::new(4.0, 5.0).unwrap();
+
+        let actual = estimator.calc_beta_diff(0.0, &params1, &params2).unwrap();
+
+        // B(5, 7) / (B(2, 3) * B(4, 5)) = 16 / 11.
+        assert_approx_eq!(actual, (16f64 / 11f64).ln(), 1e-12);
+    }
+
     #[test]
     fn test_appell_f1_stable() {
         let answers = vec![