diff --git a/book/src/advanced_usage.md b/book/src/advanced_usage.md
index 6745d3a..40781da 100644
--- a/book/src/advanced_usage.md
+++ b/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
diff --git a/book/src/dmr_scoring_details.md b/book/src/dmr_scoring_details.md
index 7822c67..31ab48d 100644
--- a/book/src/dmr_scoring_details.md
+++ b/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.
-
diff --git a/modkit-core/src/dmr/beta_diff.rs b/modkit-core/src/dmr/beta_diff.rs
index 1c85355..5d26503 100644
--- a/modkit-core/src/dmr/beta_diff.rs
+++ b/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![
diff --git a/modkit-core/src/dmr/single_site.rs b/modkit-core/src/dmr/single_site.rs
index 6e0f6e6..1d7474d 100644
--- a/modkit-core/src/dmr/single_site.rs
+++ b/modkit-core/src/dmr/single_site.rs
@@ -17,7 +17,7 @@ use rustc_hash::FxHashMap;
 use crate::dmr::beta_diff::{BetaParams, PMapEstimator};
 use crate::dmr::llr_model::{llk_ratio, AggregatedCounts};
 use crate::dmr::tabix::{
-    MultiSampleIndex, SampleToChromBMLines, SingleSiteSampleIndex,
+    MultiSampleIndex, SampleCount, SampleToChromBMLines, SingleSiteSampleIndex,
 };
 use crate::dmr::util::{cohen_h, DmrBatchOfPositions};
 use crate::errs::{MkError, MkResult};
@@ -75,7 +75,7 @@ impl SingleSiteDmrAnalysis {
         }
         let prior = if let Some(raw_prior_params) = prior {
             if raw_prior_params[0] + raw_prior_params[1] < 1.0 {
-                bail!("alpha + beta must be > 1.0 for numerical stability")
+                bail!("alpha + beta must be >= 1.0 for numerical stability")
             }
             let prior =
                 BetaParams::new(raw_prior_params[0], raw_prior_params[1])?;
@@ -121,7 +121,7 @@ impl SingleSiteDmrAnalysis {
             prior,
             rope,
             cap_coverages,
-        ));
+        )?);
 
         Ok(Self {
             sample_index,
@@ -570,25 +570,25 @@ impl SingleSiteDmrScore {
     }
 
     fn new_multi(
-        counts_a: &[AggregatedCounts],
-        counts_b: &[AggregatedCounts],
+        counts_a: &[SampleCount],
+        counts_b: &[SampleCount],
         sample_index: &SingleSiteSampleIndex,
         position: u64,
         strand: Strand,
         estimator: &PMapEstimator,
     ) -> MkResult<Self> {
         let (replicate_epmap, replicate_effect_sizes) = if sample_index
-            .matched_replicate_samples()
-            && counts_a.len() == counts_b.len()
+            .has_complete_positive_matched_counts(counts_a, counts_b)
         {
             let n_samples = counts_a.len();
             let mut replicate_epmap = Vec::with_capacity(n_samples);
             let mut replicate_effect_size = Vec::with_capacity(n_samples);
             for (a, b) in counts_a.iter().zip(counts_b) {
-                let epmap = estimator.predict(a, b).map_err(|e| {
-                    debug!("failed to calculate MAP-based p-value, {e}");
-                    MkError::BetaDiffCalcError
-                })?;
+                let epmap =
+                    estimator.predict(&a.counts, &b.counts).map_err(|e| {
+                        debug!("failed to calculate MAP-based p-value, {e}");
+                        MkError::BetaDiffCalcError
+                    })?;
                 replicate_epmap.push(epmap.e_pmap);
                 replicate_effect_size.push(epmap.effect_size);
             }
@@ -596,14 +596,18 @@ impl SingleSiteDmrScore {
         } else {
             (Vec::new(), Vec::new())
         };
-        let pct_a_samples = ((counts_a.len() as f32
-            / sample_index.num_a_samples() as f32)
-            * 100f32)
-            .floor() as usize;
-        let pct_b_samples = ((counts_b.len() as f32
-            / sample_index.num_b_samples() as f32)
-            * 100f32)
-            .floor() as usize;
+        let positive_a =
+            counts_a.iter().filter(|sample| sample.counts.total > 0).count();
+        let positive_b =
+            counts_b.iter().filter(|sample| sample.counts.total > 0).count();
+        let pct_a_samples = represented_sample_percentage(
+            positive_a,
+            sample_index.num_a_samples(),
+        );
+        let pct_b_samples = represented_sample_percentage(
+            positive_b,
+            sample_index.num_b_samples(),
+        );
         let balanced_counts_a = collapse_counts(counts_a, true);
         let balanced_counts_b = collapse_counts(counts_b, true);
         let epmap_balanced = estimator
@@ -781,33 +785,196 @@ impl SingleSiteDmrScore {
     }
 }
 
-fn collapse_counts(
-    counts: &[AggregatedCounts],
-    balance: bool,
-) -> AggregatedCounts {
-    if counts.len() == 1 {
-        counts[0].clone()
+fn represented_sample_percentage(
+    positive_count: usize,
+    configured_count: usize,
+) -> usize {
+    debug_assert!(configured_count > 0);
+    debug_assert!(positive_count <= configured_count);
+    let percentage = (positive_count as u128 * 100) / configured_count as u128;
+    usize::try_from(percentage)
+        .expect("represented sample percentage must fit in usize")
+}
+
+fn collapse_counts(counts: &[SampleCount], balance: bool) -> AggregatedCounts {
+    let positive_count =
+        counts.iter().filter(|sample| sample.counts.total > 0).count();
+    if positive_count == 0 {
+        AggregatedCounts::zero()
+    } else if positive_count == 1 {
+        counts
+            .iter()
+            .find(|sample| sample.counts.total > 0)
+            .map(|sample| sample.counts.clone())
+            .unwrap()
     } else if balance {
-        let total_cov = counts.iter().map(|ac| ac.total).sum::<usize>();
-        let n = counts.len();
-        let target_cov = total_cov as f32 / n as f32;
-        counts.iter().fold(AggregatedCounts::zero(), |agg, next| {
-            let counts = next
-                .iter_mod_fractions()
-                .map(|(code, frac)| {
-                    (code, (frac * target_cov).floor() as usize)
-                })
-                .collect::<HashMap<ModCodeRepr, usize>>();
-            let total = target_cov as usize;
-            let ac = AggregatedCounts::try_new(counts, total).unwrap();
-            agg.op(&ac)
-        })
+        let total_cov = counts
+            .iter()
+            .filter(|sample| sample.counts.total > 0)
+            .map(|sample_count| sample_count.counts.total)
+            .sum::<usize>();
+        let target_cov = total_cov as f32 / positive_count as f32;
+        counts.iter().filter(|sample| sample.counts.total > 0).fold(
+            AggregatedCounts::zero(),
+            |agg, next| {
+                let counts = next
+                    .counts
+                    .iter_mod_fractions()
+                    .map(|(code, frac)| {
+                        (code, (frac * target_cov).floor() as usize)
+                    })
+                    .collect::<HashMap<ModCodeRepr, usize>>();
+                let total = target_cov as usize;
+                let ac = AggregatedCounts::try_new(counts, total).unwrap();
+                agg.op(&ac)
+            },
+        )
     } else {
-        counts.iter().fold(AggregatedCounts::zero(), |agg, next| agg.op(next))
+        counts
+            .iter()
+            .filter(|sample| sample.counts.total > 0)
+            .fold(AggregatedCounts::zero(), |agg, next| agg.op(&next.counts))
+    }
+}
+
+#[cfg(test)]
+mod positive_coverage_tests {
+    use std::collections::HashMap;
+
+    use super::SingleSiteDmrScore;
+    use crate::dmr::beta_diff::{BetaParams, PMapEstimator};
+    use crate::dmr::llr_model::AggregatedCounts;
+    use crate::dmr::tabix::{
+        MultiSampleIndex, SampleCount, SingleSiteSampleIndex,
+    };
+    use crate::errs::MkError;
+    use crate::mod_base_code::ModCodeRepr;
+    use crate::util::Strand;
+    use rustc_hash::FxHashMap;
+
+    fn counts(code: char, n_mod: usize, total: usize) -> AggregatedCounts {
+        AggregatedCounts::try_new(
+            HashMap::from([(ModCodeRepr::Code(code), n_mod)]),
+            total,
+        )
+        .unwrap()
+    }
+
+    fn score(
+        counts_a: &[(char, usize, usize)],
+        counts_b: &[(char, usize, usize)],
+        num_a: usize,
+        num_b: usize,
+    ) -> Result<SingleSiteDmrScore, MkError> {
+        let sample_index = SingleSiteSampleIndex::new(
+            MultiSampleIndex::new(Vec::new(), FxHashMap::default(), 0, 1),
+            num_a,
+            num_b,
+            None,
+        )
+        .unwrap();
+        let counts_a = counts_a
+            .iter()
+            .enumerate()
+            .map(|(sample_id, &(code, n_mod, total))| SampleCount {
+                sample_id,
+                counts: counts(code, n_mod, total),
+            })
+            .collect::<Vec<_>>();
+        let counts_b = counts_b
+            .iter()
+            .enumerate()
+            .map(|(sample_id, &(code, n_mod, total))| SampleCount {
+                sample_id: sample_id + num_a,
+                counts: counts(code, n_mod, total),
+            })
+            .collect::<Vec<_>>();
+        let estimator = PMapEstimator::new(
+            [10, 10],
+            num_a,
+            num_b,
+            BetaParams::new(0.55, 0.55).unwrap(),
+            0.05,
+            true,
+        )
+        .unwrap();
+        SingleSiteDmrScore::new_multi(
+            &counts_a,
+            &counts_b,
+            &sample_index,
+            0,
+            Strand::Positive,
+            &estimator,
+        )
+    }
+
+    #[test]
+    fn balanced_score_ignores_zero_coverage_replicates() {
+        let score = score(&[('m', 10, 10), ('h', 0, 0)], &[('m', 0, 10)], 2, 1)
+            .unwrap();
+
+        assert_eq!(score.counts_a.total, 10);
+        assert_eq!(score.counts_a.modified_counts(), 10);
+        assert_eq!(score.counts_a.string_counts(), "m:10");
+        assert_eq!(score.counts_b.total, 10);
+        assert_eq!(score.effect_size, 1.0);
+        assert_eq!(score.balanced_effect_size, score.effect_size);
+        assert_eq!(score.pct_a_samples, 50);
+        assert_eq!(score.pct_b_samples, 100);
+    }
+
+    #[test]
+    fn reverse_balanced_score_ignores_zero_coverage_replicates() {
+        let score = score(&[('m', 0, 10)], &[('m', 10, 10), ('h', 0, 0)], 1, 2)
+            .unwrap();
+
+        assert_eq!(score.effect_size, -1.0);
+        assert_eq!(score.balanced_effect_size, score.effect_size);
+        assert_eq!(score.pct_a_samples, 100);
+        assert_eq!(score.pct_b_samples, 50);
+    }
+
+    #[test]
+    fn positive_canonical_and_absent_replicates_keep_their_meaning() {
+        let score =
+            score(&[('m', 0, 10), ('h', 0, 0)], &[('m', 0, 10)], 3, 2).unwrap();
+
+        assert_eq!(score.counts_a.total, 10);
+        assert_eq!(score.counts_a.modified_counts(), 0);
+        assert_eq!(score.effect_size, 0.0);
+        assert_eq!(score.balanced_effect_size, score.effect_size);
+        assert_eq!(score.pct_a_samples, 33);
+        assert_eq!(score.pct_b_samples, 50);
+    }
+
+    #[test]
+    fn condition_without_positive_coverage_is_a_recoverable_site_failure() {
+        let result = score(&[('m', 0, 0), ('h', 0, 0)], &[('m', 0, 10)], 2, 1);
+
+        assert!(matches!(result, Err(MkError::BetaDiffCalcError)));
+    }
+
+    #[test]
+    fn represented_sample_percentages_use_exact_integer_floor() {
+        for represented in [53, 59] {
+            let counts_a = vec![('m', 10, 10); represented];
+            let score = score(&counts_a, &[('m', 0, 10)], 100, 1).unwrap();
+
+            assert_eq!(score.pct_a_samples, represented);
+            assert_eq!(score.pct_b_samples, 100);
+        }
     }
 }
 
 type ChromToSingleScores = (String, Vec<MkResult<SingleSiteDmrScore>>);
+
+fn sort_chrom_to_site_scores(chrom_to_site_scores: &mut [ChromToSingleScores]) {
+    // The batch maps use hash iteration order, while the writer and stateful
+    // segmenter require the lexical contig order used by SingleSiteBatches.
+    chrom_to_site_scores
+        .sort_unstable_by(|(a_chrom, _), (b_chrom, _)| a_chrom.cmp(b_chrom));
+}
+
 fn process_batch_of_positions(
     batch: DmrBatchOfPositions,
     sample_index: Arc<SingleSiteSampleIndex>,
@@ -816,7 +983,7 @@ fn process_batch_of_positions(
     let (a_lines, b_lines) =
         sample_index.read_bedmethyl_lines_organized_by_position(batch)?;
 
-    let chrom_to_site_scores = a_lines
+    let mut chrom_to_site_scores = a_lines
         .into_iter()
         // intersect a_lines and b_lines on contig/chrom, there should be a
         // filter upstream of this to make sure that this is not ever a miss
@@ -852,9 +1019,32 @@ fn process_batch_of_positions(
         })
         .collect::<Vec<ChromToSingleScores>>();
 
+    sort_chrom_to_site_scores(&mut chrom_to_site_scores);
+
     Ok(chrom_to_site_scores)
 }
 
+#[cfg(test)]
+mod chrom_score_order_tests {
+    use super::{sort_chrom_to_site_scores, ChromToSingleScores};
+
+    #[test]
+    fn chrom_scores_are_sorted_lexically() {
+        let mut chrom_to_site_scores = ["cc", "aa", "bb"]
+            .into_iter()
+            .map(|chrom| (chrom.to_string(), Vec::new()))
+            .collect::<Vec<ChromToSingleScores>>();
+
+        sort_chrom_to_site_scores(&mut chrom_to_site_scores);
+
+        let chroms = chrom_to_site_scores
+            .iter()
+            .map(|(chrom, _)| chrom.as_str())
+            .collect::<Vec<_>>();
+        assert_eq!(chroms, vec!["aa", "bb", "cc"]);
+    }
+}
+
 struct Coverages {
     a_coverages: Vec<u64>,
     b_coverages: Vec<u64>,
@@ -1319,7 +1509,7 @@ fn path_to_region_labels(
     path: &[States],
     positions: &[u64],
 ) -> Vec<(u64, u64, States)> {
-    assert_eq!(path.len(), positions.len() - 1);
+    assert_eq!(path.len(), positions.len());
     if path.is_empty() {
         return Vec::new();
     } else {
@@ -1345,3 +1535,75 @@ fn path_to_region_labels(
         agg
     }
 }
+
+#[cfg(test)]
+mod segmentation_path_tests {
+    use super::{path_to_region_labels, States};
+
+    fn assert_every_position_covered_once(
+        path: &[States],
+        positions: &[u64],
+        regions: &[(u64, u64, States)],
+    ) {
+        for (&position, &state) in positions.iter().zip(path) {
+            let covering_regions = regions
+                .iter()
+                .filter(|(start, end, _)| *start <= position && position < *end)
+                .collect::<Vec<_>>();
+            assert_eq!(
+                covering_regions.len(),
+                1,
+                "position {position} in state {state:?} was not covered exactly once"
+            );
+            assert_eq!(covering_regions[0].2, state);
+        }
+    }
+
+    fn check_case(
+        path: &[States],
+        positions: &[u64],
+        expected: &[(u64, u64, States)],
+    ) {
+        let regions = path_to_region_labels(path, positions);
+        assert_eq!(regions, expected);
+        assert_every_position_covered_once(path, positions, &regions);
+    }
+
+    #[test]
+    fn singleton_position_becomes_a_single_base_region() {
+        check_case(&[States::Same], &[10], &[(10, 11, States::Same)]);
+    }
+
+    #[test]
+    fn two_positions_are_both_integrated() {
+        check_case(
+            &[States::Same, States::Same],
+            &[10, 20],
+            &[(10, 21, States::Same)],
+        );
+        check_case(
+            &[States::Same, States::Different],
+            &[10, 20],
+            &[(10, 11, States::Same), (20, 21, States::Different)],
+        );
+    }
+
+    #[test]
+    fn multi_position_transitions_cover_each_position_once() {
+        check_case(
+            &[
+                States::Same,
+                States::Same,
+                States::Different,
+                States::Different,
+                States::Same,
+            ],
+            &[10, 20, 30, 40, 50],
+            &[
+                (10, 21, States::Same),
+                (30, 41, States::Different),
+                (50, 51, States::Same),
+            ],
+        );
+    }
+}
diff --git a/modkit-core/src/dmr/subcommands.rs b/modkit-core/src/dmr/subcommands.rs
index 66658e1..0c63d82 100644
--- a/modkit-core/src/dmr/subcommands.rs
+++ b/modkit-core/src/dmr/subcommands.rs
@@ -275,7 +275,8 @@ pub struct PairwiseDmr {
     /// 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.
+    /// is calculated. Alpha and beta must each be positive, and their sum
+    /// must be >= 1.0.
 
     #[clap(help_heading = "Single-site Options")]
     #[arg(
diff --git a/modkit-core/src/dmr/tabix.rs b/modkit-core/src/dmr/tabix.rs
index 6f67056..7155c90 100644
--- a/modkit-core/src/dmr/tabix.rs
+++ b/modkit-core/src/dmr/tabix.rs
@@ -9,7 +9,7 @@ use rustc_hash::{FxHashMap, FxHashSet};
 use crate::dmr::bedmethyl::{aggregate_counts2, BedMethylLine};
 use crate::dmr::llr_model::AggregatedCounts;
 use crate::dmr::util::{n_choose_2, DmrBatch, DmrBatchOfPositions};
-use crate::errs::MkResult;
+use crate::errs::{MkError, MkResult};
 use crate::genome_positions::StrandedPosition;
 use crate::mod_base_code::{DnaBase, ModCodeRepr};
 use crate::monoid::Moniod;
@@ -21,11 +21,14 @@ pub(super) type ChromToSampleBMLines =
 /// Sample id -> {Chrom -> <BedMethyl_records>}
 pub(super) type SampleToChromBMLines =
     FxHashMap<usize, FxHashMap<String, Vec<BedMethylLine>>>;
+#[derive(Debug)]
+pub(super) struct SampleCount {
+    pub(super) sample_id: usize,
+    pub(super) counts: AggregatedCounts,
+}
 /// Chrom -> {StrandedPosition -> <X_i>} for all i in samples
-pub(super) type ChromToPosAggregatedCounts = FxHashMap<
-    String,
-    BTreeMap<StrandedPosition<DnaBase>, Vec<AggregatedCounts>>,
->;
+pub(super) type ChromToPosAggregatedCounts =
+    FxHashMap<String, BTreeMap<StrandedPosition<DnaBase>, Vec<SampleCount>>>;
 /// Usually (control, experiment)
 pub(super) type BedMethylLinesResult<T> = MkResult<(T, T)>;
 
@@ -206,13 +209,21 @@ impl SingleSiteSampleIndex {
 
     fn organize_bedmethy_lines(
         &self,
-        sample: SampleToChromBMLines,
+        mut sample: SampleToChromBMLines,
+        configured_sample_ids: &[usize],
         code_lookup: &FxHashMap<ModCodeRepr, DnaBase>,
     ) -> MkResult<ChromToPosAggregatedCounts> {
         let mut agg = FxHashMap::default();
 
         // samples should be length ~1-5
-        for chrom_to_filtered_bm_records in sample.into_values() {
+        for sample_id in configured_sample_ids {
+            let chrom_to_filtered_bm_records =
+                sample.remove(sample_id).ok_or_else(|| {
+                    MkError::InvalidBedMethyl(format!(
+                        "missing configured sample ID {sample_id} while \
+                         organizing bedMethyl records"
+                    ))
+                })?;
             let chrom_to_counts = chrom_to_filtered_bm_records
                 .into_iter()
                 .map(|(chrom, lines)| {
@@ -264,12 +275,23 @@ impl SingleSiteSampleIndex {
                         Ok(aggregated_counts) => chrom_agg
                             .entry(position)
                             .or_insert(Vec::new())
-                            .push(aggregated_counts),
+                            .push(SampleCount {
+                                sample_id: *sample_id,
+                                counts: aggregated_counts,
+                            }),
                         Err(e) => return Err(e),
                     }
                 }
             }
         }
+        if !sample.is_empty() {
+            let mut unexpected_ids = sample.keys().copied().collect::<Vec<_>>();
+            unexpected_ids.sort_unstable();
+            return Err(MkError::InvalidBedMethyl(format!(
+                "found unexpected sample IDs while organizing bedMethyl \
+                 records: {unexpected_ids:?}"
+            )));
+        }
         Ok(agg)
     }
 
@@ -341,10 +363,12 @@ impl SingleSiteSampleIndex {
         // group by chrom, this can fail if the records are deemed invalid
         let counts_a = self.organize_bedmethy_lines(
             bedmethyl_lines_a,
+            &self.control_idxs,
             &self.multi_sample_index.code_lookup,
         )?;
         let counts_b = self.organize_bedmethy_lines(
             bedmethyl_lines_b,
+            &self.exp_idxs,
             &self.multi_sample_index.code_lookup,
         )?;
         Ok((counts_a, counts_b))
@@ -367,8 +391,93 @@ impl SingleSiteSampleIndex {
         self.num_a_samples() == self.num_b_samples() && self.num_a_samples() > 1
     }
 
+    pub(super) fn has_complete_positive_matched_counts(
+        &self,
+        counts_a: &[SampleCount],
+        counts_b: &[SampleCount],
+    ) -> bool {
+        let aligned_and_positive =
+            |expected_ids: &[usize], counts: &[SampleCount]| {
+                expected_ids.len() == counts.len()
+                    && expected_ids.iter().zip(counts).all(
+                        |(expected_id, sample_count)| {
+                            expected_id == &sample_count.sample_id
+                                && sample_count.counts.total > 0
+                        },
+                    )
+            };
+        self.matched_replicate_samples()
+            && aligned_and_positive(&self.control_idxs, counts_a)
+            && aligned_and_positive(&self.exp_idxs, counts_b)
+    }
+
     #[inline]
     pub(super) fn multiple_samples(&self) -> bool {
         self.num_a_samples() > 1 || self.num_b_samples() > 1
     }
 }
+
+#[cfg(test)]
+mod sample_identity_tests {
+    use rustc_hash::FxHashMap;
+
+    use super::{
+        MultiSampleIndex, SampleToChromBMLines, SingleSiteSampleIndex,
+    };
+    use crate::dmr::bedmethyl::BedMethylLine;
+    use crate::mod_base_code::{DnaBase, ModCodeRepr};
+
+    fn line(modified: usize) -> BedMethylLine {
+        BedMethylLine::parse(&format!(
+            "chr1\t0\t1\tm\t10\t+\t0\t1\t255,0,0\t10\t0.00\t{}\t{}\t0\t0\t0\t0\t0",
+            modified,
+            10 - modified
+        ))
+        .unwrap()
+    }
+
+    #[test]
+    fn organizer_preserves_configured_sample_order() {
+        let mut code_lookup = FxHashMap::default();
+        code_lookup.insert(ModCodeRepr::Code('m'), DnaBase::C);
+        let sample_index = SingleSiteSampleIndex::new(
+            MultiSampleIndex::new(Vec::new(), code_lookup.clone(), 0, 1),
+            3,
+            3,
+            None,
+        )
+        .unwrap();
+        let mut samples = SampleToChromBMLines::default();
+        for (sample_id, modified) in [(5, 9), (3, 1), (4, 5)] {
+            samples.insert(
+                sample_id,
+                FxHashMap::from_iter([(
+                    "chr1".to_string(),
+                    vec![line(modified)],
+                )]),
+            );
+        }
+
+        let organized = sample_index
+            .organize_bedmethy_lines(
+                samples,
+                &sample_index.exp_idxs,
+                &code_lookup,
+            )
+            .unwrap();
+        let counts = organized.get("chr1").unwrap().values().next().unwrap();
+
+        assert_eq!(
+            counts
+                .iter()
+                .map(|sample_count| {
+                    (
+                        sample_count.sample_id,
+                        sample_count.counts.modified_counts(),
+                    )
+                })
+                .collect::<Vec<_>>(),
+            vec![(3, 1), (4, 5), (5, 9)]
+        );
+    }
+}
diff --git a/modkit-core/src/hmm.rs b/modkit-core/src/hmm.rs
index 5a86669..9d91d78 100644
--- a/modkit-core/src/hmm.rs
+++ b/modkit-core/src/hmm.rs
@@ -187,7 +187,7 @@ impl HmmModel {
         assert_eq!(probs.len(), transitions.len());
         let (dp_matrix, pointers) = self.viterbi_forward(&probs, &transitions);
         let path = self.viterbi_decode(&dp_matrix, &pointers);
-        assert_eq!(path.len(), scores.len() - 1);
+        assert_eq!(path.len(), scores.len());
         path
     }
 
@@ -197,21 +197,16 @@ impl HmmModel {
         pointers: &[PointerCell],
     ) -> Vec<States> {
         let final_state = dp_matrix.last().unwrap().argmax();
-        // dbg!(final_state);
         let mut path = vec![final_state];
-        let mut curr_pointer =
-            pointers.last().unwrap().get_value(final_state).unwrap();
-        for pointers in pointers.iter().rev().skip(1) {
-            let pointer = pointers.get_value(curr_pointer);
-            if let Some(pointer) = pointer {
-                path.push(pointer);
-                curr_pointer = pointer;
-            } else {
-                break;
-            }
+        let mut current_state = final_state;
+        // The first pointer cell is the empty start cell, and the second
+        // points back to the un-emitted start state. Decode only the emitted
+        // states, from the final score back through the second score.
+        for pointers in pointers.iter().skip(2).rev() {
+            current_state = pointers.get_value(current_state).unwrap();
+            path.push(current_state);
         }
 
-        path.pop();
         path.reverse();
         path
     }
@@ -413,7 +408,106 @@ impl Projection {
 
 #[cfg(test)]
 mod hmm_tests {
-    use crate::hmm::HmmModel;
+    use crate::hmm::{HmmModel, States};
+
+    fn test_model() -> HmmModel {
+        HmmModel::new(0.1, 0.9, 0.3, -0.1, 0.01, 500, true).unwrap()
+    }
+
+    fn emission_score(model: &HmmModel, p: f64, state: States) -> f64 {
+        let p = if p == 0f64 { 1e-5 } else { p };
+        let (factor, log_probability) = match state {
+            States::Same => (model.same_state_factor, p.ln()),
+            States::Different => {
+                (model.diff_state_factor, (1f64 - p + 1e-5).ln())
+            }
+        };
+        factor * (log_probability - model.significance_factor)
+    }
+
+    fn transition_score(
+        model: &HmmModel,
+        previous: States,
+        current: States,
+        diff_stay: f64,
+    ) -> f64 {
+        match (previous, current) {
+            (States::Same, States::Same) => model.same_to_same,
+            (States::Same, States::Different) => model.same_to_diff,
+            (States::Different, States::Same) => (1f64 - diff_stay).ln(),
+            (States::Different, States::Different) => diff_stay.ln(),
+        }
+    }
+
+    /// Exhaustively score the hidden start state and every emitted state.
+    /// This deliberately does not use the dynamic-programming matrix or its
+    /// back-pointers, so it independently checks both state order and length.
+    fn brute_force_path(
+        model: &HmmModel,
+        scores: &[f64],
+        positions: &[u64],
+    ) -> Vec<States> {
+        assert!(!scores.is_empty());
+        assert_eq!(scores.len(), positions.len());
+
+        let probabilities = scores
+            .iter()
+            .map(|&score| (-score.max(0f64)).exp())
+            .collect::<Vec<_>>();
+        let diff_stays = positions.windows(2).fold(
+            vec![model.dmr_prior],
+            |mut transitions, window| {
+                let gap = (window[1] - window[0]) as f64;
+                transitions.push(if model.linear_proj {
+                    model.projection.linear_project_prob(gap)
+                } else {
+                    model.projection.ln_project_prob(gap)
+                });
+                transitions
+            },
+        );
+
+        let state_count = scores.len() + 1;
+        let mut best: Option<(f64, Vec<States>)> = None;
+        for encoded_path in 0..(1usize << state_count) {
+            let states = (0..state_count)
+                .map(|i| {
+                    if encoded_path & (1usize << i) == 0 {
+                        States::Same
+                    } else {
+                        States::Different
+                    }
+                })
+                .collect::<Vec<_>>();
+            let initial_score = match states[0] {
+                States::Same => model.same_to_same,
+                States::Different => model.same_to_diff,
+            };
+            let total_score = probabilities
+                .iter()
+                .zip(diff_stays.iter())
+                .enumerate()
+                .fold(initial_score, |total, (i, (&p, &diff_stay))| {
+                    total
+                        + transition_score(
+                            model,
+                            states[i],
+                            states[i + 1],
+                            diff_stay,
+                        )
+                        + emission_score(model, p, states[i + 1])
+                });
+
+            if best
+                .as_ref()
+                .map(|(best_score, _)| total_score > *best_score)
+                .unwrap_or(true)
+            {
+                best = Some((total_score, states[1..].to_vec()));
+            }
+        }
+        best.unwrap().1
+    }
 
     #[test]
     fn test_prob_to_factor() {
@@ -421,4 +515,34 @@ mod hmm_tests {
         let fact = HmmModel::prob_to_factor(sig_fact).unwrap();
         dbg!(fact);
     }
+
+    #[test]
+    fn viterbi_path_matches_independent_oracle_for_tiny_sequences() {
+        let model = test_model();
+        let cases = [
+            (vec![0.0], vec![10]),
+            (vec![12.0], vec![10]),
+            (vec![0.0, 12.0], vec![10, 20]),
+            (vec![12.0, 0.0], vec![10, 20]),
+            (vec![0.0, 0.0, 12.0, 12.0, 0.0], vec![10, 20, 30, 40, 50]),
+        ];
+
+        for (scores, positions) in cases {
+            let expected = brute_force_path(&model, &scores, &positions);
+            let actual = model.viterbi_path(&scores, &positions);
+            assert_eq!(actual.len(), scores.len());
+            assert_eq!(actual, expected, "scores: {scores:?}");
+        }
+    }
+
+    #[test]
+    fn viterbi_path_preserves_state_transitions_in_order() {
+        let model = test_model();
+        let scores = vec![0.0, 0.0, 12.0, 12.0, 0.0];
+        let positions = vec![10, 20, 30, 40, 50];
+        let expected = brute_force_path(&model, &scores, &positions);
+        assert!(expected.windows(2).any(|states| states[0] != states[1]));
+
+        assert_eq!(model.viterbi_path(&scores, &positions), expected);
+    }
 }
diff --git a/modkit/tests/test_dmr.rs b/modkit/tests/test_dmr.rs
index 5a2991d..9e3fae8 100644
--- a/modkit/tests/test_dmr.rs
+++ b/modkit/tests/test_dmr.rs
@@ -1,9 +1,29 @@
 use crate::common::{
     check_against_expected_text_file, check_legal_csv, run_modkit,
 };
+use std::fs;
+use std::path::Path;
+use std::process::{Command, Output};
 
 mod common;
 
+const ZERO_COVERAGE_BED: &str =
+    "../tests/resources/dmr_zero_coverage_chr1.bed.gz";
+const CANONICAL_ONLY_BED: &str =
+    "../tests/resources/dmr_canonical_only_chr1.bed.gz";
+const MODIFIED_10X_BED: &str =
+    "../tests/resources/dmr_10x_modified_chr1.bed.gz";
+const HALF_MODIFIED_10X_BED: &str =
+    "../tests/resources/dmr_10x_half_modified_chr1.bed.gz";
+const CANONICAL_10X_BED: &str =
+    "../tests/resources/dmr_10x_canonical_chr1.bed.gz";
+const EMPTY_BED: &str = "../tests/resources/dmr_empty_chr1.bed.gz";
+const ZERO_COVERAGE_REF: &str = "../tests/resources/dmr_zero_coverage_chr1.fa";
+const ZERO_PERCENTILE_BED: &str =
+    "../tests/resources/dmr_zero_percentile_21_sites_chr1.bed.gz";
+const ZERO_PERCENTILE_REF: &str =
+    "../tests/resources/dmr_zero_percentile_21_sites_chr1.fa";
+
 #[test]
 fn test_dmr_helps() {
     let _ = run_modkit(&["dmr", "pair", "--help"])
@@ -75,6 +95,839 @@ fn test_dmr_regression() {
     );
 }
 
+fn run_dmr_with_prior(output: &Path, alpha: &str, beta: &str) -> Output {
+    Command::new(env!("CARGO_BIN_EXE_modkit"))
+        .args([
+            "dmr",
+            "pair",
+            "-a",
+            "../tests/resources/lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz",
+            "-b",
+            "../tests/resources/lung_00733-m_primary-tumour_5mc-5hmc_chr20_cpg_pileup.bed.gz",
+            "-o",
+            output.to_str().unwrap(),
+            "--ref",
+            "../tests/resources/GRCh38_chr20.fa",
+            "--base",
+            "C",
+            "--prior",
+            alpha,
+            beta,
+            "--delta",
+            "1",
+            "--max-coverages",
+            "100",
+            "100",
+            "--threads",
+            "1",
+            "--io-threads",
+            "1",
+            "--suppress-progress",
+            "--force",
+        ])
+        .output()
+        .unwrap()
+}
+
+#[test]
+fn dmr_prior_cli_accepts_boundary_and_interior_but_rejects_invalid_inputs() {
+    let temp_dir = tempfile::tempdir().unwrap();
+
+    for (label, alpha, beta) in
+        [("boundary", "0.5", "0.5"), ("interior", "0.55", "0.55")]
+    {
+        let output_path = temp_dir.path().join(format!("{label}.bed"));
+        let output = run_dmr_with_prior(&output_path, alpha, beta);
+        assert!(
+            output.status.success(),
+            "{label} prior ({alpha}, {beta}) was rejected: {}",
+            String::from_utf8_lossy(&output.stderr)
+        );
+    }
+
+    let below_boundary = run_dmr_with_prior(
+        &temp_dir.path().join("below-boundary.bed"),
+        "0.4",
+        "0.5",
+    );
+    assert!(!below_boundary.status.success());
+    assert!(String::from_utf8_lossy(&below_boundary.stderr)
+        .contains("alpha + beta must be >= 1.0 for numerical stability"));
+
+    let non_positive =
+        run_dmr_with_prior(&temp_dir.path().join("non-positive.bed"), "0", "1");
+    assert!(!non_positive.status.success());
+    assert!(String::from_utf8_lossy(&non_positive.stderr)
+        .contains("invalid beta parameters 0, 1"));
+}
+
+fn run_dmr_with_max_coverages(
+    output: &Path,
+    max_coverages: [usize; 2],
+    prior: Option<(&str, &str)>,
+) -> Output {
+    let mut command = Command::new(env!("CARGO_BIN_EXE_modkit"));
+    command.args([
+        "dmr",
+        "pair",
+        "-a",
+        MODIFIED_10X_BED,
+        "-b",
+        CANONICAL_10X_BED,
+        "-o",
+        output.to_str().unwrap(),
+        "--ref",
+        ZERO_COVERAGE_REF,
+        "--base",
+        "C",
+        "--header",
+        "--threads",
+        "1",
+        "--io-threads",
+        "1",
+        "--suppress-progress",
+        "--force",
+    ]);
+    command
+        .arg("--max-coverages")
+        .arg(max_coverages[0].to_string())
+        .arg(max_coverages[1].to_string());
+    if let Some((alpha, beta)) = prior {
+        command.args(["--prior", alpha, beta]);
+    }
+    command.output().unwrap()
+}
+
+fn assert_no_scientific_dmr_rows(output_path: &Path) {
+    if output_path.exists() {
+        let output = fs::read_to_string(output_path).unwrap();
+        assert!(
+            output.lines().all(|line| line.is_empty() || line.starts_with('#')),
+            "unexpected scientific DMR output:\n{output}"
+        );
+    }
+}
+
+#[test]
+fn explicit_zero_max_coverage_fails_clearly_without_scientific_rows() {
+    let temp_dir = tempfile::tempdir().unwrap();
+
+    for (prior_label, prior) in
+        [("default", None), ("boundary", Some(("0.5", "0.5")))]
+    {
+        for max_coverages in [[0, 0], [0, 10], [10, 0]] {
+            let output_path = temp_dir.path().join(format!(
+                "{prior_label}-{}-{}.bed",
+                max_coverages[0], max_coverages[1]
+            ));
+            let output =
+                run_dmr_with_max_coverages(&output_path, max_coverages, prior);
+            let stderr = String::from_utf8_lossy(&output.stderr);
+            assert!(!output.status.success(), "{stderr}");
+            assert!(
+                stderr.contains(&format!(
+                    "resolved maximum coverage must be greater than zero for \
+                     both conditions, got control {} and experiment {}",
+                    max_coverages[0], max_coverages[1]
+                )),
+                "{stderr}"
+            );
+            assert!(!stderr.contains("NaN"), "{stderr}");
+            assert_no_scientific_dmr_rows(&output_path);
+        }
+    }
+}
+
+#[test]
+fn explicit_positive_max_coverage_still_produces_exact_finite_output() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let output_path = temp_dir.path().join("positive-10-10.bed");
+
+    let output = run_dmr_with_max_coverages(&output_path, [10, 10], None);
+    let stderr = String::from_utf8_lossy(&output.stderr);
+    assert!(output.status.success(), "{stderr}");
+    assert!(
+        stderr.contains("processed 1 sites successfully, 0 failed"),
+        "{stderr}"
+    );
+
+    let output = fs::read_to_string(output_path).unwrap();
+    assert!(!output.contains("NaN"), "{output}");
+    let lines = output.lines().collect::<Vec<_>>();
+    assert_eq!(lines.len(), 2, "{output}");
+    let header = lines[0].split('\t').collect::<Vec<_>>();
+    let row = lines[1].split('\t').collect::<Vec<_>>();
+    let field = |name| header.iter().position(|field| *field == name).unwrap();
+    assert_eq!(row[field("a_counts")], "m:10");
+    assert_eq!(row[field("a_total")], "10");
+    assert_eq!(row[field("b_counts")], "m:0");
+    assert_eq!(row[field("b_total")], "10");
+    assert_eq!(row[field("a_pct_modified")], "1");
+    assert_eq!(row[field("b_pct_modified")], "0");
+    assert_eq!(row[field("effect_size")], "1");
+    assert_eq!(row[field("map_pvalue")], "0.0000006230948043897833");
+}
+
+fn run_automatic_zero_percentile_dmr(
+    output: &Path,
+    threads: &str,
+    io_threads: &str,
+    max_coverages: Option<[&str; 2]>,
+) -> Output {
+    let mut command = Command::new(env!("CARGO_BIN_EXE_modkit"));
+    command.args([
+        "dmr",
+        "pair",
+        "-a",
+        ZERO_PERCENTILE_BED,
+        "-b",
+        ZERO_PERCENTILE_BED,
+        "-o",
+        output.to_str().unwrap(),
+        "--ref",
+        ZERO_PERCENTILE_REF,
+        "--base",
+        "C",
+        "--header",
+        "--threads",
+        threads,
+        "--io-threads",
+        io_threads,
+        "--suppress-progress",
+        "--force",
+    ]);
+    if let Some([control, experiment]) = max_coverages {
+        command.args(["--max-coverages", control, experiment]);
+    } else {
+        command.args(["-N", "21", "--interval-size", "4"]);
+    }
+    command.output().unwrap()
+}
+
+#[test]
+fn automatic_zero_percentile_cap_fails_and_explicit_positive_cap_recovers() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let mut error_lines = Vec::new();
+
+    for (label, threads, io_threads) in
+        [("one", "1", "1"), ("four-a", "4", "2"), ("four-b", "4", "2")]
+    {
+        let output_path = temp_dir.path().join(format!("{label}.bed"));
+        let output = run_automatic_zero_percentile_dmr(
+            &output_path,
+            threads,
+            io_threads,
+            None,
+        );
+        let stderr = String::from_utf8_lossy(&output.stderr);
+        assert!(!output.status.success(), "{stderr}");
+        assert!(!stderr.contains("NaN"), "{stderr}");
+        let error_line = stderr
+            .lines()
+            .find(|line| {
+                line.contains(
+                    "resolved maximum coverage must be greater than zero",
+                )
+            })
+            .unwrap_or_else(|| {
+                panic!("missing maximum-coverage error: {stderr}")
+            });
+        error_lines.push(error_line.to_string());
+        assert_no_scientific_dmr_rows(&output_path);
+    }
+
+    assert!(error_lines.iter().all(|line| line == &error_lines[0]));
+    assert!(error_lines[0].contains("control 0 and experiment 0"));
+
+    let output_path = temp_dir.path().join("explicit-positive.bed");
+    let output = run_automatic_zero_percentile_dmr(
+        &output_path,
+        "1",
+        "1",
+        Some(["10", "10"]),
+    );
+    let stderr = String::from_utf8_lossy(&output.stderr);
+    assert!(output.status.success(), "{stderr}");
+    let beta_error_line = stderr
+        .lines()
+        .find(|line| line.contains("beta-diff-calc-error"))
+        .unwrap_or_else(|| panic!("missing beta-diff error count: {stderr}"));
+    assert!(beta_error_line.contains("20"), "{stderr}");
+    assert!(
+        stderr.contains("processed 1 sites successfully, 20 failed"),
+        "{stderr}"
+    );
+    assert!(!stderr.contains("NaN"), "{stderr}");
+
+    let output = fs::read_to_string(output_path).unwrap();
+    assert!(!output.contains("NaN"), "{output}");
+    let lines = output.lines().collect::<Vec<_>>();
+    assert_eq!(lines.len(), 2, "{output}");
+    let header = lines[0].split('\t').collect::<Vec<_>>();
+    let row = lines[1].split('\t').collect::<Vec<_>>();
+    assert_eq!(header.len(), row.len());
+    let field = |name| header.iter().position(|field| *field == name).unwrap();
+    assert_eq!(row[field("#chrom")], "chr1");
+    assert_eq!(row[field("start")], "20");
+    assert_eq!(row[field("end")], "21");
+    assert_eq!(row[field("a_counts")], "C:0");
+    assert_eq!(row[field("a_total")], "1");
+    assert_eq!(row[field("b_counts")], "C:0");
+    assert_eq!(row[field("b_total")], "1");
+    assert_eq!(row[field("map_pvalue")], "1");
+    assert_eq!(row[field("effect_size")], "0");
+}
+
+fn run_zero_coverage_dmr(
+    output: &Path,
+    prior: Option<(&str, &str)>,
+    delta: Option<&str>,
+) -> Output {
+    let mut command = Command::new(env!("CARGO_BIN_EXE_modkit"));
+    command.args([
+        "dmr",
+        "pair",
+        "-a",
+        ZERO_COVERAGE_BED,
+        "-b",
+        CANONICAL_ONLY_BED,
+        "-o",
+        output.to_str().unwrap(),
+        "--ref",
+        ZERO_COVERAGE_REF,
+        "--base",
+        "C",
+        "--max-coverages",
+        "1",
+        "1",
+        "--threads",
+        "1",
+        "--io-threads",
+        "1",
+        "--suppress-progress",
+        "--force",
+    ]);
+    if let Some((alpha, beta)) = prior {
+        command.args(["--prior", alpha, beta]);
+    }
+    if let Some(delta) = delta {
+        command.args(["--delta", delta]);
+    }
+    command.output().unwrap()
+}
+
+fn assert_zero_coverage_is_failed(output: Output, output_path: &Path) {
+    let stderr = String::from_utf8_lossy(&output.stderr);
+    assert!(
+        output.status.success(),
+        "zero-coverage DMR site should be a recoverable site failure: {stderr}"
+    );
+    assert!(stderr.contains("beta-diff-calc-error"), "{stderr}");
+    assert!(
+        stderr.contains("processed 0 sites successfully, 1 failed"),
+        "{stderr}"
+    );
+
+    let output_bytes = fs::read(output_path).unwrap();
+    assert!(
+        output_bytes.is_empty(),
+        "zero-coverage site produced output: {}",
+        String::from_utf8_lossy(&output_bytes)
+    );
+}
+
+#[test]
+fn zero_coverage_site_is_failed_without_nan_under_default_prior() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let output_path = temp_dir.path().join("default-prior.bed");
+
+    let output = run_zero_coverage_dmr(&output_path, None, None);
+
+    assert_zero_coverage_is_failed(output, &output_path);
+}
+
+#[test]
+fn zero_coverage_site_does_not_abort_at_boundary_prior() {
+    let temp_dir = tempfile::tempdir().unwrap();
+
+    for (label, delta) in [("default-delta", None), ("delta-one", Some("1"))] {
+        let output_path = temp_dir.path().join(format!("{label}.bed"));
+        let output =
+            run_zero_coverage_dmr(&output_path, Some(("0.5", "0.5")), delta);
+
+        assert_zero_coverage_is_failed(output, &output_path);
+    }
+}
+
+fn run_replicate_dmr(
+    output_path: &Path,
+    samples_a: &[&str],
+    samples_b: &[&str],
+    min_coverage: &str,
+    threads: &str,
+    io_threads: &str,
+) -> Output {
+    let mut command = Command::new(env!("CARGO_BIN_EXE_modkit"));
+    command.args(["dmr", "pair"]);
+    for sample in samples_a {
+        command.args(["-a", sample]);
+    }
+    for sample in samples_b {
+        command.args(["-b", sample]);
+    }
+    command.args([
+        "-o",
+        output_path.to_str().unwrap(),
+        "--ref",
+        ZERO_COVERAGE_REF,
+        "--base",
+        "C",
+        "--header",
+        "--max-coverages",
+        "10",
+        "10",
+        "--min-valid-coverage",
+        min_coverage,
+        "--threads",
+        threads,
+        "--io-threads",
+        io_threads,
+        "--suppress-progress",
+        "--force",
+    ]);
+    command.output().unwrap()
+}
+
+fn read_single_successful_dmr_row(
+    output: Output,
+    output_path: &Path,
+) -> String {
+    let stderr = String::from_utf8_lossy(&output.stderr);
+    assert!(output.status.success(), "{stderr}");
+    assert!(
+        stderr.contains("processed 1 sites successfully, 0 failed"),
+        "{stderr}"
+    );
+    let output = fs::read_to_string(output_path).unwrap();
+    assert_eq!(output.lines().count(), 2, "{output}");
+    output
+}
+
+fn dmr_field<'a>(output: &'a str, field_name: &str) -> &'a str {
+    let mut lines = output.lines();
+    let header = lines.next().unwrap().split('\t').collect::<Vec<_>>();
+    let row = lines.next().unwrap().split('\t').collect::<Vec<_>>();
+    assert_eq!(header.len(), row.len(), "{output}");
+    let index = header.iter().position(|field| *field == field_name).unwrap();
+    row[index]
+}
+
+#[test]
+fn complete_matched_replicates_follow_cli_order_and_independent_oracles() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let samples_a =
+        [MODIFIED_10X_BED, CANONICAL_10X_BED, HALF_MODIFIED_10X_BED];
+    let samples_b =
+        [CANONICAL_10X_BED, HALF_MODIFIED_10X_BED, MODIFIED_10X_BED];
+    let expected_maps = [
+        "0.0000006230948043897833",
+        "0.020888505198542493",
+        "0.01995939262165838",
+    ];
+    let expected_effects = ["1", "-0.5", "-0.5"];
+
+    for (pair_index, ((sample_a, sample_b), (map, effect))) in samples_a
+        .iter()
+        .zip(samples_b.iter())
+        .zip(expected_maps.iter().zip(expected_effects.iter()))
+        .enumerate()
+    {
+        let output_path =
+            temp_dir.path().join(format!("independent-{pair_index}.bed"));
+        let output = run_replicate_dmr(
+            &output_path,
+            &[*sample_a],
+            &[*sample_b],
+            "0",
+            "1",
+            "1",
+        );
+        let output = read_single_successful_dmr_row(output, &output_path);
+        assert_eq!(dmr_field(&output, "map_pvalue"), *map);
+        assert_eq!(dmr_field(&output, "effect_size"), *effect);
+    }
+
+    let expected_maps = expected_maps.join(",");
+    let expected_effects = expected_effects.join(",");
+    let mut outputs = Vec::new();
+    for (threads, io_threads) in [("1", "1"), ("4", "2")] {
+        let output_path = temp_dir
+            .path()
+            .join(format!("complete-{threads}-{io_threads}.bed"));
+        let output = run_replicate_dmr(
+            &output_path,
+            &samples_a,
+            &samples_b,
+            "0",
+            threads,
+            io_threads,
+        );
+        let output = read_single_successful_dmr_row(output, &output_path);
+        assert_eq!(dmr_field(&output, "replicate_map_pvalues"), expected_maps);
+        assert_eq!(
+            dmr_field(&output, "replicate_effect_sizes"),
+            expected_effects
+        );
+        for (field, expected) in [
+            ("a_counts", "m:15"),
+            ("a_total", "30"),
+            ("b_counts", "m:15"),
+            ("b_total", "30"),
+            ("a_pct_modified", "0.5"),
+            ("b_pct_modified", "0.5"),
+            ("map_pvalue", "1"),
+            ("effect_size", "0"),
+            ("balanced_map_pvalue", "1"),
+            ("balanced_effect_size", "0"),
+            ("pct_a_samples", "100"),
+            ("pct_b_samples", "100"),
+        ] {
+            assert_eq!(dmr_field(&output, field), expected);
+        }
+        outputs.push(output);
+    }
+    assert_eq!(outputs[0].as_bytes(), outputs[1].as_bytes());
+}
+
+#[test]
+fn incomplete_matched_replicates_suppress_both_fields_and_keep_group_row() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let complete_b =
+        [CANONICAL_10X_BED, HALF_MODIFIED_10X_BED, MODIFIED_10X_BED];
+    let scenarios = [
+        (
+            "aligned-missing",
+            [MODIFIED_10X_BED, EMPTY_BED, HALF_MODIFIED_10X_BED],
+            [CANONICAL_10X_BED, EMPTY_BED, MODIFIED_10X_BED],
+            "0",
+        ),
+        (
+            "crossed-missing",
+            [EMPTY_BED, CANONICAL_10X_BED, HALF_MODIFIED_10X_BED],
+            [CANONICAL_10X_BED, EMPTY_BED, MODIFIED_10X_BED],
+            "0",
+        ),
+        (
+            "zero-member",
+            [MODIFIED_10X_BED, ZERO_COVERAGE_BED, HALF_MODIFIED_10X_BED],
+            complete_b,
+            "0",
+        ),
+        (
+            "threshold-filtered",
+            [MODIFIED_10X_BED, CANONICAL_ONLY_BED, HALF_MODIFIED_10X_BED],
+            complete_b,
+            "2",
+        ),
+    ];
+
+    for (label, samples_a, samples_b, min_coverage) in scenarios {
+        let mut outputs = Vec::new();
+        for (threads, io_threads) in [("1", "1"), ("4", "2")] {
+            let output_path = temp_dir
+                .path()
+                .join(format!("{label}-{threads}-{io_threads}.bed"));
+            let output = run_replicate_dmr(
+                &output_path,
+                &samples_a,
+                &samples_b,
+                min_coverage,
+                threads,
+                io_threads,
+            );
+            let output = read_single_successful_dmr_row(output, &output_path);
+            assert_eq!(
+                dmr_field(&output, "replicate_map_pvalues"),
+                "-",
+                "{label}: {output}"
+            );
+            assert_eq!(
+                dmr_field(&output, "replicate_effect_sizes"),
+                "-",
+                "{label}: {output}"
+            );
+            outputs.push(output);
+        }
+        assert_eq!(outputs[0].as_bytes(), outputs[1].as_bytes(), "{label}");
+    }
+}
+
+#[test]
+fn positive_coverage_oracle_is_filter_and_thread_invariant() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let scenarios = [
+        ("zero", [ZERO_COVERAGE_BED, MODIFIED_10X_BED], "0"),
+        ("min-one", [ZERO_COVERAGE_BED, MODIFIED_10X_BED], "1"),
+        ("absent", [EMPTY_BED, MODIFIED_10X_BED], "0"),
+    ];
+    let mut outputs = Vec::new();
+
+    for (label, samples_a, min_coverage) in scenarios {
+        for (threads, io_threads) in [("1", "1"), ("4", "2")] {
+            let output_path = temp_dir
+                .path()
+                .join(format!("{label}-{threads}-{io_threads}.bed"));
+            let output = run_replicate_dmr(
+                &output_path,
+                &samples_a,
+                &[CANONICAL_10X_BED],
+                min_coverage,
+                threads,
+                io_threads,
+            );
+            let output = read_single_successful_dmr_row(output, &output_path);
+            for (field, expected) in [
+                ("a_counts", "m:10"),
+                ("a_total", "10"),
+                ("b_counts", "m:0"),
+                ("b_total", "10"),
+                ("map_pvalue", "0.0000006230948043897833"),
+                ("effect_size", "1"),
+                ("balanced_map_pvalue", "0.0000006230948043897833"),
+                ("balanced_effect_size", "1"),
+                ("pct_a_samples", "50"),
+                ("pct_b_samples", "100"),
+            ] {
+                assert_eq!(dmr_field(&output, field), expected);
+            }
+            assert!(!output.contains("h:0"), "{label}: {output}");
+            outputs.push(output);
+        }
+    }
+
+    for output in &outputs[1..] {
+        assert_eq!(&outputs[0], output);
+    }
+}
+
+fn run_segmented_dmr(
+    output: &Path,
+    segments: &Path,
+    threads: &str,
+    io_threads: &str,
+) {
+    run_modkit(&[
+        "dmr",
+        "pair",
+        "-a",
+        "../tests/resources/\
+         lung_00733-m_adjacent-normal_5mc-5hmc_chr20_cpg_pileup.bed.gz",
+        "-b",
+        "../tests/resources/\
+         lung_00733-m_primary-tumour_5mc-5hmc_chr20_cpg_pileup.bed.gz",
+        "-o",
+        output.to_str().unwrap(),
+        "--segment",
+        segments.to_str().unwrap(),
+        "--ref",
+        "../tests/resources/GRCh38_chr20.fa",
+        "--header",
+        "--base",
+        "C",
+        "--max-coverages",
+        "100",
+        "100",
+        "--threads",
+        threads,
+        "--io-threads",
+        io_threads,
+        "--suppress-progress",
+        "--force",
+    ])
+    .expect("segmented DMR run should succeed");
+}
+
+#[test]
+fn segmentation_includes_last_site_and_is_thread_deterministic() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let sites_one = temp_dir.path().join("sites-threads-1.bed");
+    let segments_one = temp_dir.path().join("segments-threads-1.bed");
+    let sites_four = temp_dir.path().join("sites-threads-4.bed");
+    let segments_four = temp_dir.path().join("segments-threads-4.bed");
+
+    run_segmented_dmr(&sites_one, &segments_one, "1", "1");
+    run_segmented_dmr(&sites_four, &segments_four, "4", "2");
+
+    let sites_one = fs::read_to_string(sites_one).unwrap();
+    let sites_four = fs::read_to_string(sites_four).unwrap();
+    let segments_one = fs::read_to_string(segments_one).unwrap();
+    let segments_four = fs::read_to_string(segments_four).unwrap();
+    assert_eq!(sites_one.as_bytes(), sites_four.as_bytes());
+    assert_eq!(segments_one.as_bytes(), segments_four.as_bytes());
+
+    let site_rows = sites_one.lines().skip(1).collect::<Vec<_>>();
+    let segment_rows = segments_one.lines().skip(1).collect::<Vec<_>>();
+    assert_eq!(site_rows.len(), 17_271);
+
+    let segment_site_total = segment_rows
+        .iter()
+        .map(|row| row.split('\t').nth(5).unwrap().parse::<usize>().unwrap())
+        .sum::<usize>();
+    assert_eq!(segment_site_total, site_rows.len());
+
+    let last_site_end = site_rows.last().unwrap().split('\t').nth(2).unwrap();
+    let last_segment_end =
+        segment_rows.last().unwrap().split('\t').nth(2).unwrap();
+    assert_eq!(last_site_end, "10804378");
+    assert_eq!(last_segment_end, last_site_end);
+}
+
+struct DmrOutput {
+    sites: String,
+    segments: String,
+}
+
+fn run_multi_contig_dmr(
+    output_dir: &Path,
+    label: &str,
+    interval_size: &str,
+    threads: &str,
+    io_threads: &str,
+) -> DmrOutput {
+    let sites = output_dir.join(format!("sites-{label}.bed"));
+    let segments = output_dir.join(format!("segments-{label}.bed"));
+
+    run_modkit(&[
+        "dmr",
+        "pair",
+        "-a",
+        "../tests/resources/dmr_contig_order_a.bed.gz",
+        "-b",
+        "../tests/resources/dmr_contig_order_b.bed.gz",
+        "-o",
+        sites.to_str().unwrap(),
+        "--segment",
+        segments.to_str().unwrap(),
+        "--ref",
+        "../tests/resources/dmr_contig_order.fa",
+        "--header",
+        "--base",
+        "C",
+        "--max-coverages",
+        "100",
+        "100",
+        "--interval-size",
+        interval_size,
+        "--batch-size",
+        "1",
+        "--threads",
+        threads,
+        "--io-threads",
+        io_threads,
+        "--suppress-progress",
+        "--force",
+    ])
+    .expect("multi-contig segmented DMR run should succeed");
+
+    DmrOutput {
+        sites: fs::read_to_string(sites).unwrap(),
+        segments: fs::read_to_string(segments).unwrap(),
+    }
+}
+
+fn parse_site_keys(output: &str) -> Vec<(String, u64, u64)> {
+    output
+        .lines()
+        .skip(1)
+        .map(|row| {
+            let fields = row.split('\t').collect::<Vec<_>>();
+            (
+                fields[0].to_string(),
+                fields[1].parse::<u64>().unwrap(),
+                fields[2].parse::<u64>().unwrap(),
+            )
+        })
+        .collect()
+}
+
+fn parse_segments(output: &str) -> Vec<(String, u64, u64, String, usize)> {
+    output
+        .lines()
+        .skip(1)
+        .map(|row| {
+            let fields = row.split('\t').collect::<Vec<_>>();
+            (
+                fields[0].to_string(),
+                fields[1].parse::<u64>().unwrap(),
+                fields[2].parse::<u64>().unwrap(),
+                fields[3].to_string(),
+                fields[5].parse::<usize>().unwrap(),
+            )
+        })
+        .collect()
+}
+
+#[test]
+fn multi_contig_segmentation_is_invariant_to_batch_geometry_and_threads() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let interval_ten_one =
+        run_multi_contig_dmr(temp_dir.path(), "i10-t1", "10", "1", "1");
+    let interval_ten_four =
+        run_multi_contig_dmr(temp_dir.path(), "i10-t4", "10", "4", "2");
+    let interval_three_one =
+        run_multi_contig_dmr(temp_dir.path(), "i3-t1", "3", "1", "1");
+    let interval_three_four =
+        run_multi_contig_dmr(temp_dir.path(), "i3-t4", "3", "4", "2");
+
+    for (label, other) in [
+        ("interval 10, four threads", &interval_ten_four),
+        ("interval 3, one thread", &interval_three_one),
+        ("interval 3, four threads", &interval_three_four),
+    ] {
+        assert!(
+            interval_ten_one.sites.as_bytes() == other.sites.as_bytes(),
+            "site output changed for {label}"
+        );
+        assert!(
+            interval_ten_one.segments.as_bytes() == other.segments.as_bytes(),
+            "segment output changed for {label}"
+        );
+    }
+
+    let site_keys = parse_site_keys(&interval_ten_one.sites);
+    let expected_site_keys = [("aa", 6u64), ("bb", 15u64), ("cc", 6u64)]
+        .into_iter()
+        .flat_map(|(chrom, size)| {
+            (0..size).map(move |position| {
+                (chrom.to_string(), position, position + 1)
+            })
+        })
+        .collect::<Vec<_>>();
+    assert_eq!(site_keys, expected_site_keys);
+
+    let segments = parse_segments(&interval_ten_one.segments);
+    assert_eq!(
+        segments,
+        vec![
+            ("aa".to_string(), 0, 6, "different".to_string(), 6),
+            ("bb".to_string(), 0, 15, "different".to_string(), 15),
+            ("cc".to_string(), 0, 6, "different".to_string(), 6),
+        ]
+    );
+
+    for (chrom, start, end) in site_keys {
+        let covering_segments = segments
+            .iter()
+            .filter(|(segment_chrom, segment_start, segment_end, _, _)| {
+                segment_chrom == &chrom
+                    && *segment_start <= start
+                    && end <= *segment_end
+            })
+            .count();
+        assert_eq!(covering_segments, 1, "site {chrom}:{start}-{end}");
+    }
+}
+
 // todo
 //  test pair with explicit index
 //  test multi
diff --git a/tests/resources/dmr_10x_canonical_chr1.bed b/tests/resources/dmr_10x_canonical_chr1.bed
new file mode 100644
index 0000000..b78a983
--- /dev/null
+++ b/tests/resources/dmr_10x_canonical_chr1.bed
@@ -0,0 +1 @@
+chr1	0	1	m	10	+	0	1	255,0,0	10	0.00	0	10	0	0	0	0	0
diff --git a/tests/resources/dmr_10x_canonical_chr1.bed.gz b/tests/resources/dmr_10x_canonical_chr1.bed.gz
new file mode 100644
index 0000000..f3b998a
Binary files /dev/null and b/tests/resources/dmr_10x_canonical_chr1.bed.gz differ
diff --git a/tests/resources/dmr_10x_canonical_chr1.bed.gz.tbi b/tests/resources/dmr_10x_canonical_chr1.bed.gz.tbi
new file mode 100644
index 0000000..0e8b2e6
Binary files /dev/null and b/tests/resources/dmr_10x_canonical_chr1.bed.gz.tbi differ
diff --git a/tests/resources/dmr_10x_half_modified_chr1.bed b/tests/resources/dmr_10x_half_modified_chr1.bed
new file mode 100644
index 0000000..c7f0366
--- /dev/null
+++ b/tests/resources/dmr_10x_half_modified_chr1.bed
@@ -0,0 +1 @@
+chr1	0	1	m	10	+	0	1	255,0,0	10	50.00	5	5	0	0	0	0	0
diff --git a/tests/resources/dmr_10x_half_modified_chr1.bed.gz b/tests/resources/dmr_10x_half_modified_chr1.bed.gz
new file mode 100644
index 0000000..6e52060
Binary files /dev/null and b/tests/resources/dmr_10x_half_modified_chr1.bed.gz differ
diff --git a/tests/resources/dmr_10x_half_modified_chr1.bed.gz.tbi b/tests/resources/dmr_10x_half_modified_chr1.bed.gz.tbi
new file mode 100644
index 0000000..1e36f47
Binary files /dev/null and b/tests/resources/dmr_10x_half_modified_chr1.bed.gz.tbi differ
diff --git a/tests/resources/dmr_10x_modified_chr1.bed b/tests/resources/dmr_10x_modified_chr1.bed
new file mode 100644
index 0000000..6769184
--- /dev/null
+++ b/tests/resources/dmr_10x_modified_chr1.bed
@@ -0,0 +1 @@
+chr1	0	1	m	10	+	0	1	255,0,0	10	100.00	10	0	0	0	0	0	0
diff --git a/tests/resources/dmr_10x_modified_chr1.bed.gz b/tests/resources/dmr_10x_modified_chr1.bed.gz
new file mode 100644
index 0000000..07f2b13
Binary files /dev/null and b/tests/resources/dmr_10x_modified_chr1.bed.gz differ
diff --git a/tests/resources/dmr_10x_modified_chr1.bed.gz.tbi b/tests/resources/dmr_10x_modified_chr1.bed.gz.tbi
new file mode 100644
index 0000000..9aed4bf
Binary files /dev/null and b/tests/resources/dmr_10x_modified_chr1.bed.gz.tbi differ
diff --git a/tests/resources/dmr_canonical_only_chr1.bed.gz b/tests/resources/dmr_canonical_only_chr1.bed.gz
new file mode 100644
index 0000000..7153aa2
Binary files /dev/null and b/tests/resources/dmr_canonical_only_chr1.bed.gz differ
diff --git a/tests/resources/dmr_canonical_only_chr1.bed.gz.tbi b/tests/resources/dmr_canonical_only_chr1.bed.gz.tbi
new file mode 100644
index 0000000..e8fdab9
Binary files /dev/null and b/tests/resources/dmr_canonical_only_chr1.bed.gz.tbi differ
diff --git a/tests/resources/dmr_contig_order.fa b/tests/resources/dmr_contig_order.fa
new file mode 100644
index 0000000..ac83f38
--- /dev/null
+++ b/tests/resources/dmr_contig_order.fa
@@ -0,0 +1,6 @@
+>aa
+CCCCCC
+>bb
+CCCCCCCCCCCCCCC
+>cc
+CCCCCC
diff --git a/tests/resources/dmr_contig_order.fa.fai b/tests/resources/dmr_contig_order.fa.fai
new file mode 100644
index 0000000..84bdbae
--- /dev/null
+++ b/tests/resources/dmr_contig_order.fa.fai
@@ -0,0 +1,3 @@
+aa	6	4	6	7
+bb	15	15	15	16
+cc	6	35	6	7
diff --git a/tests/resources/dmr_contig_order_a.bed b/tests/resources/dmr_contig_order_a.bed
new file mode 100644
index 0000000..653af2d
--- /dev/null
+++ b/tests/resources/dmr_contig_order_a.bed
@@ -0,0 +1,27 @@
+aa	0	1	C	10	+	0	1	255,0,0	10	100.00	10	0	0	0	0	0	0
+aa	1	2	C	10	+	1	2	255,0,0	10	100.00	10	0	0	0	0	0	0
+aa	2	3	C	10	+	2	3	255,0,0	10	100.00	10	0	0	0	0	0	0
+aa	3	4	C	10	+	3	4	255,0,0	10	100.00	10	0	0	0	0	0	0
+aa	4	5	C	10	+	4	5	255,0,0	10	100.00	10	0	0	0	0	0	0
+aa	5	6	C	10	+	5	6	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	0	1	C	10	+	0	1	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	1	2	C	10	+	1	2	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	2	3	C	10	+	2	3	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	3	4	C	10	+	3	4	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	4	5	C	10	+	4	5	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	5	6	C	10	+	5	6	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	6	7	C	10	+	6	7	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	7	8	C	10	+	7	8	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	8	9	C	10	+	8	9	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	9	10	C	10	+	9	10	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	10	11	C	10	+	10	11	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	11	12	C	10	+	11	12	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	12	13	C	10	+	12	13	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	13	14	C	10	+	13	14	255,0,0	10	100.00	10	0	0	0	0	0	0
+bb	14	15	C	10	+	14	15	255,0,0	10	100.00	10	0	0	0	0	0	0
+cc	0	1	C	10	+	0	1	255,0,0	10	100.00	10	0	0	0	0	0	0
+cc	1	2	C	10	+	1	2	255,0,0	10	100.00	10	0	0	0	0	0	0
+cc	2	3	C	10	+	2	3	255,0,0	10	100.00	10	0	0	0	0	0	0
+cc	3	4	C	10	+	3	4	255,0,0	10	100.00	10	0	0	0	0	0	0
+cc	4	5	C	10	+	4	5	255,0,0	10	100.00	10	0	0	0	0	0	0
+cc	5	6	C	10	+	5	6	255,0,0	10	100.00	10	0	0	0	0	0	0
diff --git a/tests/resources/dmr_contig_order_a.bed.gz b/tests/resources/dmr_contig_order_a.bed.gz
new file mode 100644
index 0000000..4e4c92a
Binary files /dev/null and b/tests/resources/dmr_contig_order_a.bed.gz differ
diff --git a/tests/resources/dmr_contig_order_a.bed.gz.tbi b/tests/resources/dmr_contig_order_a.bed.gz.tbi
new file mode 100644
index 0000000..241c105
Binary files /dev/null and b/tests/resources/dmr_contig_order_a.bed.gz.tbi differ
diff --git a/tests/resources/dmr_contig_order_b.bed b/tests/resources/dmr_contig_order_b.bed
new file mode 100644
index 0000000..521a2bd
--- /dev/null
+++ b/tests/resources/dmr_contig_order_b.bed
@@ -0,0 +1,27 @@
+aa	0	1	C	10	+	0	1	255,0,0	10	0.00	0	10	0	0	0	0	0
+aa	1	2	C	10	+	1	2	255,0,0	10	0.00	0	10	0	0	0	0	0
+aa	2	3	C	10	+	2	3	255,0,0	10	0.00	0	10	0	0	0	0	0
+aa	3	4	C	10	+	3	4	255,0,0	10	0.00	0	10	0	0	0	0	0
+aa	4	5	C	10	+	4	5	255,0,0	10	0.00	0	10	0	0	0	0	0
+aa	5	6	C	10	+	5	6	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	0	1	C	10	+	0	1	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	1	2	C	10	+	1	2	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	2	3	C	10	+	2	3	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	3	4	C	10	+	3	4	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	4	5	C	10	+	4	5	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	5	6	C	10	+	5	6	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	6	7	C	10	+	6	7	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	7	8	C	10	+	7	8	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	8	9	C	10	+	8	9	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	9	10	C	10	+	9	10	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	10	11	C	10	+	10	11	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	11	12	C	10	+	11	12	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	12	13	C	10	+	12	13	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	13	14	C	10	+	13	14	255,0,0	10	0.00	0	10	0	0	0	0	0
+bb	14	15	C	10	+	14	15	255,0,0	10	0.00	0	10	0	0	0	0	0
+cc	0	1	C	10	+	0	1	255,0,0	10	0.00	0	10	0	0	0	0	0
+cc	1	2	C	10	+	1	2	255,0,0	10	0.00	0	10	0	0	0	0	0
+cc	2	3	C	10	+	2	3	255,0,0	10	0.00	0	10	0	0	0	0	0
+cc	3	4	C	10	+	3	4	255,0,0	10	0.00	0	10	0	0	0	0	0
+cc	4	5	C	10	+	4	5	255,0,0	10	0.00	0	10	0	0	0	0	0
+cc	5	6	C	10	+	5	6	255,0,0	10	0.00	0	10	0	0	0	0	0
diff --git a/tests/resources/dmr_contig_order_b.bed.gz b/tests/resources/dmr_contig_order_b.bed.gz
new file mode 100644
index 0000000..36a355e
Binary files /dev/null and b/tests/resources/dmr_contig_order_b.bed.gz differ
diff --git a/tests/resources/dmr_contig_order_b.bed.gz.tbi b/tests/resources/dmr_contig_order_b.bed.gz.tbi
new file mode 100644
index 0000000..64ed8bc
Binary files /dev/null and b/tests/resources/dmr_contig_order_b.bed.gz.tbi differ
diff --git a/tests/resources/dmr_empty_chr1.bed b/tests/resources/dmr_empty_chr1.bed
new file mode 100644
index 0000000..61d1e44
--- /dev/null
+++ b/tests/resources/dmr_empty_chr1.bed
@@ -0,0 +1 @@
+# Intentionally contains no bedMethyl records.
diff --git a/tests/resources/dmr_empty_chr1.bed.gz b/tests/resources/dmr_empty_chr1.bed.gz
new file mode 100644
index 0000000..7c49c82
Binary files /dev/null and b/tests/resources/dmr_empty_chr1.bed.gz differ
diff --git a/tests/resources/dmr_empty_chr1.bed.gz.tbi b/tests/resources/dmr_empty_chr1.bed.gz.tbi
new file mode 100644
index 0000000..ca981d8
Binary files /dev/null and b/tests/resources/dmr_empty_chr1.bed.gz.tbi differ
diff --git a/tests/resources/dmr_zero_coverage_chr1.bed b/tests/resources/dmr_zero_coverage_chr1.bed
new file mode 100644
index 0000000..9435428
--- /dev/null
+++ b/tests/resources/dmr_zero_coverage_chr1.bed
@@ -0,0 +1 @@
+chr1	0	1	h	0	+	0	1	255,0,0	0	0.00	0	0	0	0	0	0	0
diff --git a/tests/resources/dmr_zero_coverage_chr1.bed.gz b/tests/resources/dmr_zero_coverage_chr1.bed.gz
new file mode 100644
index 0000000..a5af8c3
Binary files /dev/null and b/tests/resources/dmr_zero_coverage_chr1.bed.gz differ
diff --git a/tests/resources/dmr_zero_coverage_chr1.bed.gz.tbi b/tests/resources/dmr_zero_coverage_chr1.bed.gz.tbi
new file mode 100644
index 0000000..b9fdda7
Binary files /dev/null and b/tests/resources/dmr_zero_coverage_chr1.bed.gz.tbi differ
diff --git a/tests/resources/dmr_zero_coverage_chr1.fa b/tests/resources/dmr_zero_coverage_chr1.fa
new file mode 100644
index 0000000..cebd535
--- /dev/null
+++ b/tests/resources/dmr_zero_coverage_chr1.fa
@@ -0,0 +1,2 @@
+>chr1
+C
diff --git a/tests/resources/dmr_zero_coverage_chr1.fa.fai b/tests/resources/dmr_zero_coverage_chr1.fa.fai
new file mode 100644
index 0000000..26d3676
--- /dev/null
+++ b/tests/resources/dmr_zero_coverage_chr1.fa.fai
@@ -0,0 +1 @@
+chr1	1	6	1	2
diff --git a/tests/resources/dmr_zero_percentile_21_sites_chr1.bed b/tests/resources/dmr_zero_percentile_21_sites_chr1.bed
new file mode 100644
index 0000000..94724a2
--- /dev/null
+++ b/tests/resources/dmr_zero_percentile_21_sites_chr1.bed
@@ -0,0 +1,21 @@
+chr1	0	1	C	0	+	0	1	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	1	2	C	0	+	1	2	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	2	3	C	0	+	2	3	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	3	4	C	0	+	3	4	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	4	5	C	0	+	4	5	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	5	6	C	0	+	5	6	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	6	7	C	0	+	6	7	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	7	8	C	0	+	7	8	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	8	9	C	0	+	8	9	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	9	10	C	0	+	9	10	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	10	11	C	0	+	10	11	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	11	12	C	0	+	11	12	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	12	13	C	0	+	12	13	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	13	14	C	0	+	13	14	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	14	15	C	0	+	14	15	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	15	16	C	0	+	15	16	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	16	17	C	0	+	16	17	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	17	18	C	0	+	17	18	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	18	19	C	0	+	18	19	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	19	20	C	0	+	19	20	255,0,0	0	0.00	0	0	0	0	0	0	0
+chr1	20	21	C	1	+	20	21	255,0,0	1	0.00	0	1	0	0	0	0	0
diff --git a/tests/resources/dmr_zero_percentile_21_sites_chr1.bed.gz b/tests/resources/dmr_zero_percentile_21_sites_chr1.bed.gz
new file mode 100644
index 0000000..730b035
Binary files /dev/null and b/tests/resources/dmr_zero_percentile_21_sites_chr1.bed.gz differ
diff --git a/tests/resources/dmr_zero_percentile_21_sites_chr1.bed.gz.tbi b/tests/resources/dmr_zero_percentile_21_sites_chr1.bed.gz.tbi
new file mode 100644
index 0000000..e9f1a14
Binary files /dev/null and b/tests/resources/dmr_zero_percentile_21_sites_chr1.bed.gz.tbi differ
diff --git a/tests/resources/dmr_zero_percentile_21_sites_chr1.fa b/tests/resources/dmr_zero_percentile_21_sites_chr1.fa
new file mode 100644
index 0000000..555d429
--- /dev/null
+++ b/tests/resources/dmr_zero_percentile_21_sites_chr1.fa
@@ -0,0 +1,2 @@
+>chr1
+CCCCCCCCCCCCCCCCCCCCC
diff --git a/tests/resources/dmr_zero_percentile_21_sites_chr1.fa.fai b/tests/resources/dmr_zero_percentile_21_sites_chr1.fa.fai
new file mode 100644
index 0000000..35ad3d6
--- /dev/null
+++ b/tests/resources/dmr_zero_percentile_21_sites_chr1.fa.fai
@@ -0,0 +1 @@
+chr1	21	6	21	22