diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7aa426b..0720007 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,21 @@ All notable changes to WASP2 will be documented in this file.
 ## [Unreleased]
 
 ### Fixed
+- Compare validation count columns exactly before floating-point conversion,
+  including UInt64 values above 2**53. Distinguish archived narrow CSV input
+  from explicitly widened reference counts; retain strict legacy discrepancies.
+- Document Rust counting/bulk inference and Python/SciPy single-cell inference,
+  and correct the single-cell AnnData axes. Precision regressions cover balanced
+  likelihoods, overflow and conditional unphased symmetry.
+- Single-cell sparse counts and totals use UInt64 so H5AD output preserves values
+  above 65,535. Metadata conversion no longer requires an undeclared PyArrow install.
+- Automatic single-sample selection now applies the same heterozygote filtering
+  and feature reindexing as explicit selection, correcting inherited Python behavior.
+- Single-cell counting normalizes repeated, comma-separated and file-based sample options.
+- Shared-SNP differential analysis applies minimum depth per region, as archived
+  Python did. Unphased Python likelihoods preserve its sequential arithmetic while
+  retaining log-space fallback for underflow. Expanded comparisons explicitly widen
+  the archive's UInt16 accumulators; this also avoids incidental float32 dispatch.
 - Donor feature deduplication uses the parsed first GT allele, as Python does.
   Original GT tokens remain available for the backend's pre-dedup phasing
   decision; fitting weights, raw counts, QC and provenance reflect retained rows.
@@ -15,8 +30,8 @@ All notable changes to WASP2 will be documented in this file.
   checks, parses numeric first alleles consistently, and rejects malformed mixed
   phased input instead of silently changing a region's likelihood. Interior
   dispersion values are no longer clipped to the endpoint regularization range.
-- Full 32-run public ATAC comparisons execute archived Python in pinned Linux
-  CI, including CLI inference and donor routes. Heterogeneous donor checks cover
+- Full 32-run public ATAC comparisons execute archived Python on wide-integer
+  DataFrames in pinned Linux CI; narrow CSV precision is checked separately, including CLI inference and donor routes. Heterogeneous donor checks cover
   both dispersion models, phasing modes and pseudocounts with independent fits.
 - Container checks exercise reordered ATAC headers, region grouping, phased
   likelihoods and raw counts. Images carry package-version and source-revision
diff --git a/README.md b/README.md
index a42e994..f1abf85 100644
--- a/README.md
+++ b/README.md
@@ -56,6 +56,11 @@ WASP2 installs three command-line entry points:
 - `wasp2-count`
 - `wasp2-analyze`
 
+The Python CLI coordinates compiled Rust counting and bulk statistical inference.
+Single-cell counting uses Rust; single-cell imbalance and differential inference
+currently use Python/SciPy. A Python command entry point does not change which
+backend performs the calculations.
+
 ## Quick Start
 
 > **Prerequisite — normalize your VCF first.** WASP2's allelic-imbalance model is biallelic.
diff --git a/docs/source/user_guide/analysis.rst b/docs/source/user_guide/analysis.rst
index 1b11b7e..d7180b6 100644
--- a/docs/source/user_guide/analysis.rst
+++ b/docs/source/user_guide/analysis.rst
@@ -5,7 +5,14 @@ Overview
 --------
 
 ``wasp2-analyze`` runs allelic imbalance analysis on bulk or single-cell count
-outputs.
+outputs. The Python CLI calls Rust for bulk inference. Single-cell imbalance
+and differential inference currently run in Python/SciPy; single-cell counting
+uses Rust.
+
+Numerical compatibility retains the archived statistical procedures with wide
+counts and consistent precision. Untouched legacy file readers and explicitly
+widened references are evaluated separately; the latter do not claim to reproduce
+historical UInt16 overflow or mixed scalar/array precision.
 
 Commands:
 
@@ -171,16 +178,14 @@ sets the reference orientation; the first SNP is not symmetrized.
 - Beta-binomial LRT: :math:`\rho` is held at its **null-model MLE** while
   maximizing the alternative likelihood over :math:`\mu` (profile likelihood,
   df = 1). :math:`\rho` is not jointly re-estimated under :math:`H_1`.
-- Scalar dispersion is fit on :math:`(0,1)` with boundary values clamped for
-  numerical stability. The linear predictor retains the :math:`\pm 10` guard. The
-  linear model uses raw-depth standardization,
-  :math:`\rho(N)=\operatorname{expit}(d_1 + d_2(N-c)/s)`, and persists the fitted
-  center :math:`c` and positive scale :math:`s` with :math:`d_1` and :math:`d_2`.
-- FDR correction uses :func:`scipy.stats.false_discovery_control` with
-  ``method='bh'`` (Benjamini–Hochberg). ``scipy`` raises on NaN p-values,
-  but a hand-written BH loop using ``np.minimum.accumulate`` silently
-  propagates NaN through the cumulative minimum — always drop or impute
-  NaN p-values before BH correction.
+- Scalar dispersion is fit on :math:`(0,1)`; exact endpoints receive numerical
+  regularization. The linear fit uses pseudocounted depth and a :math:`\pm 10`
+  predictor guard during fitting. Stored coefficients preserve that predictor
+  using :math:`\rho(N)=\operatorname{expit}(d_1+d_2(N-c)/s)` with raw depth
+  :math:`N`, center :math:`c=-2\,\mathrm{pseudocount}` and scale :math:`s=1`.
+- Both backends use Benjamini–Hochberg FDR correction. Bulk uses the Rust
+  implementation; single-cell analysis uses
+  :func:`scipy.stats.false_discovery_control` with ``method='bh'``.
 
 See the upstream mapping-filter docs (:doc:`/methods/mapping_filter`) for the
 WASP2 canonical filter contract that applies to both remapping and counting
diff --git a/docs/source/user_guide/single_cell.rst b/docs/source/user_guide/single_cell.rst
index a243cdc..12daada 100644
--- a/docs/source/user_guide/single_cell.rst
+++ b/docs/source/user_guide/single_cell.rst
@@ -4,7 +4,15 @@ Single-Cell Analysis
 Overview
 --------
 
-WASP2 provides specialized tools for allele-specific analysis in single-cell RNA-seq (scRNA-seq) data. This guide covers the barcode file format requirements and single-cell-specific workflows.
+WASP2 supports barcode-based allele counting and group-level imbalance analysis.
+The validated workflow is scATAC-seq. Rust performs allele counting; Python/SciPy
+performs single-cell imbalance and differential inference through the Python CLI.
+The current counter does not implement UMI-aware molecule counting for scRNA-seq.
+
+Counting uses wide integer storage and inference uses wide count sums. Statistical
+compatibility follows the archived likelihood and testing procedures with consistent
+precision; it does not reproduce historical UInt16 overflow or mixed-precision
+rounding. The original and explicitly widened archive references are tested separately.
 
 Barcode File Format
 -------------------
@@ -300,12 +308,14 @@ The single-cell counting module outputs an AnnData (``.h5ad``) file containing:
 
 **Observations (obs):**
 
-* SNP information (chrom, pos, ref, alt)
-* Aggregate counts per SNP
+* Cell barcodes and cell-level annotations
 
 **Variables (var):**
 
-* Cell barcodes
+* SNP information (chrom, pos, ref, alt)
+* Aggregate counts per SNP
+
+Count layers have shape ``cells × SNPs``.
 
 **Unstructured (uns):**
 
diff --git a/scripts/compare_archive_analysis.py b/scripts/compare_archive_analysis.py
index 0bafbbe..080aa2a 100644
--- a/scripts/compare_archive_analysis.py
+++ b/scripts/compare_archive_analysis.py
@@ -9,6 +9,7 @@
 
 from __future__ import annotations
 
+import ast
 import contextlib
 import io
 import subprocess
@@ -27,11 +28,30 @@
 ROOT = Path(__file__).resolve().parents[1]
 ARCHIVE_REF = "41354fc"
 STAT_COLUMNS = ["null_ll", "alt_ll", "mu", "lrt", "pval", "fdr_pval"]
-COUNT_COLUMNS = ["ref_count", "alt_count", "N", "snp_count", "num_snps"]
-
-
-def load_archive(ref: str = ARCHIVE_REF) -> dict[str, types.ModuleType]:
-    """Load unmodified historical modules, restoring their legacy import names."""
+COUNT_COLUMNS = [
+    "ref_count",
+    "alt_count",
+    "N",
+    "snp_count",
+    "num_snps",
+    "num_snps_group1",
+    "num_snps_group2",
+]
+
+
+def load_archive(
+    ref: str = ARCHIVE_REF,
+    *,
+    widen_single_cell_sums: bool = False,
+    widen_bulk_csv_counts: bool = False,
+) -> dict[str, types.ModuleType]:
+    """Load historical modules with separately selectable count-width adapters.
+
+    The explicit dtype repair prevents count overflow and SciPy float32 dispatch.
+    Single-cell changes only sum accumulator keywords; bulk CSV changes only
+    the three count-column reader dtypes. Neither changes statistical algorithms.
+    The default executes the unmodified archive.
+    """
     names = ["as_analysis", "as_analysis_sc", "compare_ai"]
     saved = {name: sys.modules.get(name) for name in names}
     modules = {}
@@ -43,7 +63,60 @@ def load_archive(ref: str = ARCHIVE_REF) -> dict[str, types.ModuleType]:
             module = types.ModuleType(f"archive_{name}")
             modules[name] = module
             sys.modules[name] = module
-            exec(compile(source, f"{ref}:src/analysis/{name}.py", "exec"), module.__dict__)
+            tree = ast.parse(source)
+            if widen_bulk_csv_counts and name == "as_analysis":
+                changed = set()
+                function = next(
+                    node
+                    for node in tree.body
+                    if isinstance(node, ast.FunctionDef) and node.name == "get_imbalance"
+                )
+                for node in ast.walk(function):
+                    if not (
+                        isinstance(node, ast.Call)
+                        and isinstance(node.func, ast.Attribute)
+                        and isinstance(node.func.value, ast.Name)
+                        and node.func.value.id == "pd"
+                        and node.func.attr == "read_csv"
+                    ):
+                        continue
+                    for keyword in node.keywords:
+                        if keyword.arg != "dtype" or not isinstance(keyword.value, ast.Dict):
+                            continue
+                        for key, value in zip(
+                            keyword.value.keys, keyword.value.values, strict=True
+                        ):
+                            if (
+                                isinstance(key, ast.Constant)
+                                and key.value in {"ref_count", "alt_count", "other_count"}
+                                and isinstance(value, ast.Attribute)
+                                and isinstance(value.value, ast.Name)
+                                and value.value.id == "np"
+                                and value.attr == "uint16"
+                            ):
+                                value.attr = "uint64"
+                                changed.add(key.value)
+                if changed != {"ref_count", "alt_count", "other_count"}:
+                    raise ValueError("Cannot locate the three archived bulk CSV count dtypes")
+            if widen_single_cell_sums and name in {"as_analysis_sc", "compare_ai"}:
+                for node in ast.walk(tree):
+                    if not (
+                        isinstance(node, ast.Call)
+                        and isinstance(node.func, ast.Attribute)
+                        and node.func.attr == "sum"
+                    ):
+                        continue
+                    for keyword in node.keywords:
+                        value = keyword.value
+                        if (
+                            keyword.arg == "dtype"
+                            and isinstance(value, ast.Attribute)
+                            and isinstance(value.value, ast.Name)
+                            and value.value.id == "np"
+                            and value.attr == "uint16"
+                        ):
+                            value.attr = "uint64"
+            exec(compile(tree, f"{ref}:src/analysis/{name}.py", "exec"), module.__dict__)
     finally:
         for name, module in saved.items():
             if module is None:
@@ -219,12 +292,24 @@ def compare_frames(expected: pd.DataFrame, observed: pd.DataFrame, key: str = "r
         if column not in observed:
             failures.append(f"missing {column}")
             continue
-        left = expected[column].to_numpy(dtype=float)
-        right = observed[column].to_numpy(dtype=float)
-        deltas[column] = float(np.max(np.abs(left - right))) if len(left) else 0.0
         if column in COUNT_COLUMNS:
-            equal = np.array_equal(left, right)
+            # Python integers retain UInt64 values and signed/unsigned differences.
+            # Casting first to float64 hides one-count changes above 2**53.
+            left_counts = expected[column].tolist()
+            right_counts = observed[column].tolist()
+            valid = all(
+                pd.notna(a) and pd.notna(b) for a, b in zip(left_counts, right_counts, strict=True)
+            )
+            equal = valid and left_counts == right_counts
+            deltas[column] = (
+                max((abs(a - b) for a, b in zip(left_counts, right_counts, strict=True)), default=0)
+                if valid
+                else float("nan")
+            )
         else:
+            left = expected[column].to_numpy(dtype=float)
+            right = observed[column].to_numpy(dtype=float)
+            deltas[column] = float(np.max(np.abs(left - right))) if len(left) else 0.0
             equal = (
                 np.all(np.isfinite(left))
                 and np.all(np.isfinite(right))
@@ -242,6 +327,81 @@ def compare_frames(expected: pd.DataFrame, observed: pd.DataFrame, key: str = "r
     }
 
 
+def classify_reference_differences(
+    legacy: pd.DataFrame, widened: pd.DataFrame, observed: pd.DataFrame, key: str = "region"
+) -> dict:
+    """Retain strict failures and localize changes caused by a count-width adapter.
+
+    A complementary estimate is only a candidate: equivalence requires a separate
+    check of the actual anchored objective. This helper never alters an output.
+    """
+    strict = compare_frames(legacy, observed, key)
+    corrected = compare_frames(widened, observed, key)
+    frames = [
+        frame.set_index(key, drop=False).rename_axis(None) for frame in (legacy, widened, observed)
+    ]
+    if any(not frame.index.is_unique for frame in frames):
+        raise ValueError("Reference classification requires unique result identities")
+    rows = []
+    for identity in sorted(set().union(*(frame.index for frame in frames))):
+        present = [identity in frame.index for frame in frames]
+        if not all(present):
+            rows.append(
+                {
+                    "identity": identity,
+                    "stage": "tested_identity",
+                    "present": present,
+                    "reference_width_explains": present[1:] == [True, True]
+                    and corrected["status"] == "PASS",
+                }
+            )
+            continue
+        parts = [frame.loc[[identity]] for frame in frames]
+        comparison = compare_frames(parts[0], parts[2], key)
+        if comparison["status"] == "PASS":
+            continue
+        fields = comparison.get("failures", [])
+        stages = []
+        if any(field in COUNT_COLUMNS for field in fields):
+            stages.append("count_precision")
+        if any(field in fields for field in ("null_ll", "alt_ll", "lrt")):
+            stages.append("likelihood_precision")
+        effects = [field for field in ("mu", "mu1", "mu2", "combined_mu") if field in fields]
+        if effects:
+            stages.append("effect_estimate_precision")
+        if "pval" in fields:
+            stages.append("p_value_precision")
+        if "fdr_pval" in fields or any(field.startswith("significance at") for field in fields):
+            stages.append("multiple_testing")
+        candidates = [
+            field
+            for field in effects
+            if np.isclose(
+                float(parts[0][field].iloc[0]) + float(parts[2][field].iloc[0]),
+                1.0,
+                rtol=5e-5,
+                atol=5e-5,
+            )
+        ]
+        rows.append(
+            {
+                "identity": identity,
+                "stages": stages,
+                "fields": fields,
+                "complementary_estimate_candidates": candidates,
+                "reference_width_explains": compare_frames(parts[1], parts[2], key)["status"]
+                == "PASS",
+                "max_absolute_differences": comparison["max_absolute_differences"],
+            }
+        )
+    return {
+        "legacy_strict_comparison": strict,
+        "widened_comparison": corrected,
+        "mismatched_rows": rows,
+        "unexplained_rows": sum(not row["reference_width_explains"] for row in rows),
+    }
+
+
 def prepare_bulk_reference_counts(
     counts: pd.DataFrame,
     *,
diff --git a/src/analysis/__main__.py b/src/analysis/__main__.py
index a34b230..d39887e 100644
--- a/src/analysis/__main__.py
+++ b/src/analysis/__main__.py
@@ -306,9 +306,7 @@ def find_imbalance_sc(
     try:
         from .run_analysis_sc import run_ai_analysis_sc
     except (ImportError, ModuleNotFoundError) as error:
-        raise typer.BadParameter(
-            "Single-cell analysis requires: pip install 'wasp2[single-cell]'"
-        ) from error
+        raise typer.BadParameter("Single-cell analysis requires: pip install wasp2") from error
 
     run_ai_analysis_sc(
         count_file=counts,
@@ -404,9 +402,7 @@ def compare_imbalance(
     try:
         from .run_compare_ai import run_ai_comparison
     except (ImportError, ModuleNotFoundError) as error:
-        raise typer.BadParameter(
-            "Single-cell comparison requires: pip install 'wasp2[single-cell]'"
-        ) from error
+        raise typer.BadParameter("Single-cell comparison requires: pip install wasp2") from error
 
     run_ai_comparison(
         count_file=counts,
diff --git a/src/analysis/as_analysis.py b/src/analysis/as_analysis.py
index fcd5a4a..ae7b1bc 100644
--- a/src/analysis/as_analysis.py
+++ b/src/analysis/as_analysis.py
@@ -173,6 +173,19 @@ def opt_unphased_dp(
     # Get likelihood of first pos
     first_ll = opt_prob(prob, first_disp, first_ref[0], first_n[0])
 
+    # Preserve the archived sequential arithmetic while products remain normal.
+    # Equivalent log-space operations can select another symmetric optimum.
+    forward = np.atleast_1d(opt_prob(prob, phase_disp, phase_ref, phase_n, log=False))
+    reverse = np.atleast_1d(opt_prob(1 - prob, phase_disp, phase_ref, phase_n, log=False))
+    previous = 1.0
+    for p1, p2 in zip(forward, reverse, strict=True):
+        phase1, phase2 = previous * p1, previous * p2
+        previous = 0.5 * phase1 + 0.5 * phase2
+        if previous < np.finfo(float).tiny:
+            break
+    else:
+        return float(first_ll - np.log(previous))
+
     # Get LOG likelihood with regard to phasing of first pos.
     # opt_prob(..., log=True) returns the NEGATIVE log-pmf, so negate to get the
     # log-pmf itself for stable log-space accumulation.
diff --git a/src/analysis/as_analysis_sc.py b/src/analysis/as_analysis_sc.py
index e42ede4..fcfaee4 100644
--- a/src/analysis/as_analysis_sc.py
+++ b/src/analysis/as_analysis_sc.py
@@ -3,7 +3,7 @@
 Provides functions for analyzing allelic imbalance in single-cell data
 stored in AnnData format with SNP counts in layers.
 
-Orientation contract (scverse convention, see ISSUE2_ORIENTATION_DECISION §2b):
+Orientation contract (scverse convention):
 
     obs = cells (barcodes); var = SNPs (variants).
 
@@ -126,7 +126,7 @@ def get_imbalance_sc(
     n_counts = ref_counts + alt_counts
 
     # Calculate dispersion across dataset
-    def opt_disp(rho: float, ref_data: NDArray[np.uint16], n_data: NDArray[np.uint16]) -> float:
+    def opt_disp(rho: float, ref_data: NDArray[np.uint64], n_data: NDArray[np.uint64]) -> float:
         rho_safe = float(clamp_rho(rho))  # Prevent division by zero
         return float(
             -np.sum(
diff --git a/src/analysis/compare_ai.py b/src/analysis/compare_ai.py
index abc7b0a..4c67c69 100644
--- a/src/analysis/compare_ai.py
+++ b/src/analysis/compare_ai.py
@@ -114,7 +114,8 @@ def get_compared_imbalance(
         raise ValueError("Please provide 2 or more groups to compare.")
 
     # Process initial minimums for whole data dispersion
-    snp_cutoff: int = min_count + (2 * pseudocount)
+    # SNP eligibility removes zero-count observations; min_count applies per region.
+    snp_cutoff: int = 2 * pseudocount
 
     ref_counts: NDArray[np.uint64] = (
         adata.layers["ref"].sum(axis=0, dtype=np.uint64).T.A1 + pseudocount
diff --git a/src/analysis/run_analysis_sc.py b/src/analysis/run_analysis_sc.py
index 5a24212..733eba4 100644
--- a/src/analysis/run_analysis_sc.py
+++ b/src/analysis/run_analysis_sc.py
@@ -123,7 +123,7 @@ def process_adata_inputs(
         else:
             sample = sample_list[0]
 
-    else:
+    if sample is not None:
         sample_list = adata.uns["samples"]
 
         if sample not in sample_list:
diff --git a/src/counting/count_alleles_sc.py b/src/counting/count_alleles_sc.py
index eb0e0ae..51ff7cb 100644
--- a/src/counting/count_alleles_sc.py
+++ b/src/counting/count_alleles_sc.py
@@ -34,11 +34,11 @@ def _sparse_from_counts(
     shape: tuple[int, int],
 ) -> csr_matrix:
     if not counts:
-        return csr_matrix(shape, dtype=np.uint16)
+        return csr_matrix(shape, dtype=np.uint64)
     return csr_matrix(
         (list(counts.values()), list(zip(*counts.keys()))),
         shape=shape,
-        dtype=np.uint16,
+        dtype=np.uint64,
     )
 
 
@@ -189,7 +189,7 @@ def make_count_matrix(
         layers={"ref": sparse_ref, "alt": sparse_alt, "other": sparse_other},
     )
 
-    adata.obs = snp_df.to_pandas()
+    adata.obs = pd.DataFrame(snp_df.to_dict(as_series=False))
     adata.obs["ref_count"] = adata.layers["ref"].sum(axis=1).T.A1
     adata.obs["alt_count"] = adata.layers["alt"].sum(axis=1).T.A1
 
@@ -212,7 +212,7 @@ def make_count_matrix(
         )
         if "region" not in feature_df.columns:
             feature_df = feature_df.with_columns(pl.col(feature_cols[0]).alias("region"))
-        adata.uns["feature"] = feature_df.to_pandas()
+        adata.uns["feature"] = pd.DataFrame(feature_df.to_dict(as_series=False))
 
         # region_snp_dict = dict(
         #     df.join(snp_df, on=["chrom", "pos", "ref", "alt"], how="left"
diff --git a/src/counting/run_counting_sc.py b/src/counting/run_counting_sc.py
index cb19b2b..6de1ff9 100644
--- a/src/counting/run_counting_sc.py
+++ b/src/counting/run_counting_sc.py
@@ -51,19 +51,19 @@ def __init__(
         self.feature_file = feature_file
         self.use_region_names = use_region_names
 
-        # Make sure samples turned into str list
-        # Ideally single sample for single cell
-        normalized_samples: list[str] | None
-        if isinstance(samples, str):
-            # Check if sample file or comma delim string
-            if Path(samples).is_file():
-                with open(samples) as sample_file:
-                    normalized_samples = [l.strip() for l in sample_file]
-            else:
-                normalized_samples = [s.strip() for s in samples.split(",")]
-        else:
-            normalized_samples = samples
-        self.samples: list[str] | None = normalized_samples
+        # Typer supplies repeated options as a list; normalize each value equally.
+        normalized_samples: list[str] | None = None
+        if samples is not None:
+            normalized_samples = []
+            for value in [samples] if isinstance(samples, str) else samples:
+                entries = Path(value).read_text().splitlines() if Path(value).is_file() else [value]
+                normalized_samples.extend(
+                    sample.strip()
+                    for entry in entries
+                    for sample in entry.split(",")
+                    if sample.strip()
+                )
+        self.samples = normalized_samples
 
         # parse output?
         self.out_file: str = (
diff --git a/tests/analysis/fixtures/public_atac/README.md b/tests/analysis/fixtures/public_atac/README.md
index 0a413be..9feb48a 100644
--- a/tests/analysis/fixtures/public_atac/README.md
+++ b/tests/analysis/fixtures/public_atac/README.md
@@ -12,7 +12,13 @@ SNP rows and 556 peak-overlapping rows; the complete counting CLI reproduced all
 556 peak rows exactly. `metadata.json` records source and table checksums.
 
 Tests execute archived Python commit `41354fc` on these inputs rather than storing
-Rust outputs as expected results. Phased comparisons select truly phased SNPs.
+Rust outputs as expected results. The original 32-run matrix supplies wide-integer
+DataFrames to the archive, bypassing its UInt16 CSV reader. It therefore establishes
+wide-input compatibility, not exact reproduction of legacy file-input rounding.
+The bulk tests additionally execute untouched CSV input and a separately labeled
+CSV adapter changing only the three count-column dtypes to UInt64. Per-row strict
+legacy failures and widened-reference comparisons are retained in test properties;
+a legacy failure is not relabeled as a strict parity pass. Phased comparisons select truly phased SNPs.
 Linear-reference shape repairs and forwarding of its ignored phased flag are
 explicit in the test. Strict statistical tolerances remain `rtol=atol=5e-5`;
 counts must match exactly. Sixteen legacy runs plus sixteen donor feature runs
diff --git a/tests/analysis/test_archive_parity.py b/tests/analysis/test_archive_parity.py
index 0c5bcd5..f9dd250 100644
--- a/tests/analysis/test_archive_parity.py
+++ b/tests/analysis/test_archive_parity.py
@@ -1,6 +1,7 @@
 """Production-code parity with outputs frozen from unmodified Python 41354fc.
 
-The fixture is independent of Rust and does not need archived git objects in CI.
+Frozen fixtures are independent of Rust. Live reference tests also require the
+pinned archived git objects.
 These tests cover statistical compatibility on synthetic counts, not calibration.
 """
 
@@ -10,7 +11,12 @@
 import numpy as np
 import pandas as pd
 import pytest
-from scripts.compare_archive_analysis import bulk_counts, compare_frames, single_cell_counts
+from scripts.compare_archive_analysis import (
+    bulk_counts,
+    compare_frames,
+    load_archive,
+    single_cell_counts,
+)
 from typer.testing import CliRunner
 
 from analysis.as_analysis_sc import get_imbalance_sc
@@ -23,6 +29,105 @@
 )
 
 
+@pytest.mark.parametrize("scenario", ["base", "overlap", "zeros", "unequal_depth"])
+@pytest.mark.parametrize("mode", ["single", "shared", "different"])
+@pytest.mark.parametrize("phased", [False, True])
+@pytest.mark.parametrize("pc", [0, 1])
+@pytest.mark.parametrize("min_count", [0, 10, 20])
+def test_live_single_cell_archive_matrix(scenario, mode, phased, pc, min_count):
+    from scipy.sparse import csr_matrix
+
+    historical = single_cell_counts()
+    if scenario == "overlap":
+        historical.uns["feature"] = pd.concat(
+            [historical.uns["feature"], pd.DataFrame({"index": [0, 1], "region": ["overlap"] * 2})],
+            ignore_index=True,
+        )
+    elif scenario in ["zeros", "unequal_depth"]:
+        for layer in ["ref", "alt"]:
+            values = historical.layers[layer].toarray()
+            if scenario == "zeros":
+                values[0, :] = 0
+            else:
+                values[:, 2:] *= 4
+            historical.layers[layer] = csr_matrix(values)
+    # Explicit reference repair: current sums are uint64 to prevent overflow.
+    # The unchanged archive remains covered by the frozen tests below.
+    archive = load_archive(widen_single_cell_sums=True)
+    reference = {
+        "single": archive["as_analysis_sc"].get_imbalance_sc,
+        "shared": archive["compare_ai"].get_compared_imbalance,
+        "different": archive["compare_ai"].get_compared_imbalance_diff_snps,
+    }[mode]
+    current = {
+        "single": get_imbalance_sc,
+        "shared": get_compared_imbalance,
+        "different": get_compared_imbalance_diff_snps,
+    }[mode]
+    options = {
+        "sample": "sample1",
+        "groups": ["A", "B"],
+        "phased": phased,
+        "pseudocount": pc,
+        "min_count": min_count,
+    }
+    expected = reference(historical.copy(), **options)
+    actual = current(historical.T.copy(), **options)
+    assert expected.keys() == actual.keys()
+    for key in expected:
+        comparison = compare_frames(expected[key], actual[key])
+        assert comparison["status"] == "PASS", comparison
+
+
+@pytest.mark.parametrize("phased", [False, True, None])
+def test_automatic_sample_selection_filters_like_explicit_selection(phased):
+    current = single_cell_counts().T.copy()
+    current.var.loc[current.var.index[0], "sample1"] = "0/0"
+    automatic = process_adata_inputs(current.copy(), phased=phased)
+    explicit = process_adata_inputs(current.copy(), sample="sample1", phased=phased)
+    assert automatic.sample == explicit.sample == "sample1"
+    assert automatic.phased == explicit.phased
+    assert automatic.adata.n_vars == explicit.adata.n_vars == 5
+    pd.testing.assert_frame_equal(automatic.adata.var, explicit.adata.var)
+    pd.testing.assert_frame_equal(automatic.adata.uns["feature"], explicit.adata.uns["feature"])
+    np.testing.assert_array_equal(
+        automatic.adata.layers["ref"].toarray(), explicit.adata.layers["ref"].toarray()
+    )
+
+
+@pytest.mark.parametrize("sample", [None, "sample1"])
+def test_sample_selection_rejects_no_heterozygotes(sample):
+    current = single_cell_counts().T.copy()
+    current.var["sample1"] = "0/0"
+    with pytest.raises(ValueError, match="Heterozygous genotypes not found"):
+        process_adata_inputs(current, sample=sample)
+
+
+@pytest.mark.parametrize("cutoff", [0.5, 1, 3])
+def test_outlier_filtering_matches_archived_feature_reindexing(cutoff):
+    from analysis.as_analysis_sc import adata_count_qc
+
+    historical = single_cell_counts()
+    historical.obs.loc[historical.obs.index[0], "ref_count"] = 10000
+    expected = load_archive()["as_analysis_sc"].adata_count_qc(historical.copy(), z_cutoff=cutoff)
+    actual = adata_count_qc(historical.T.copy(), z_cutoff=cutoff)
+    np.testing.assert_array_equal(
+        actual.layers["ref"].toarray(), expected.layers["ref"].toarray().T
+    )
+    pd.testing.assert_frame_equal(actual.uns["feature"], expected.uns["feature"])
+
+
+@pytest.mark.parametrize("phased", [False, True, None])
+def test_missing_genotypes_force_unphased_and_missing_groups_fail(phased):
+    current = single_cell_counts().T.copy()
+    del current.uns["samples"]
+    prepared = process_adata_inputs(current.copy(), phased=phased)
+    assert prepared.phased is False
+    assert prepared.sample is None
+    with pytest.raises(KeyError, match="not found"):
+        process_adata_inputs(current, groups=["missing"], phased=phased)
+
+
 @pytest.mark.rust
 @pytest.mark.parametrize("unit", ["snv", "atac", "rna"])
 @pytest.mark.parametrize("phased", [False, True])
diff --git a/tests/analysis/test_public_atac_matrix.py b/tests/analysis/test_public_atac_matrix.py
index 62d0e59..eaaa2a9 100644
--- a/tests/analysis/test_public_atac_matrix.py
+++ b/tests/analysis/test_public_atac_matrix.py
@@ -1,7 +1,8 @@
 """Repeat the full 32-run real ATAC matrix against executed archived Python.
 
 Count tables are frozen measurements from public reads, not Rust-generated
-expectations. Linux Python 3.11 CI pins the original numerical environment.
+expectations. DataFrame, untouched CSV and widened CSV references are distinct;
+legacy precision differences are recorded separately from wide-input parity. Linux Python 3.11 CI pins the original numerical environment.
 """
 
 import contextlib
@@ -18,7 +19,12 @@
 import pytest
 import scipy
 from scipy.special import expit
-from scripts.compare_archive_analysis import adapt_linear_reference, compare_frames, load_archive
+from scripts.compare_archive_analysis import (
+    adapt_linear_reference,
+    classify_reference_differences,
+    compare_frames,
+    load_archive,
+)
 
 from analysis.run_analysis_per_donor import run_per_donor_analysis
 
@@ -35,7 +41,7 @@
 @pytest.mark.parametrize("model", ["single", "linear"])
 @pytest.mark.parametrize("phased", [False, True])
 @pytest.mark.parametrize("pc", [0, 1])
-def test_full_public_atac_matrix(tmp_path, unit, model, phased, pc):
+def test_full_public_atac_matrix(tmp_path, unit, model, phased, pc, record_property):
     import wasp2_rust
 
     source = METADATA["counts"][unit]
@@ -77,6 +83,35 @@ def capture_fit(*args, **kwargs):
     output = wasp2_rust.analyze_imbalance_run(
         str(path), method=model, phased=phased, region_col=region, pseudocount=pc
     )
+    # DataFrame input uses inferred wide counts. Also execute the original CSV
+    # reader and a separately labeled adapter changing only its count dtypes.
+    csv_results = {}
+    for wide in (False, True):
+        csv_modules = load_archive(METADATA["archive_commit"], widen_bulk_csv_counts=wide)
+        if model == "linear":
+            adapt_linear_reference(csv_modules)
+            if phased:
+                original_parse_csv = csv_modules["as_analysis"].parse_opt
+                csv_modules["as_analysis"].parse_opt = (
+                    lambda df, disp=None, phased=False, _parse=original_parse_csv: _parse(
+                        df, disp=disp, phased=True
+                    )
+                )
+        with contextlib.redirect_stdout(io.StringIO()):
+            csv_results[wide] = (
+                csv_modules["as_analysis"]
+                .get_imbalance(
+                    str(path), method=model, phased=phased, region_col=region, pseudocount=pc
+                )
+                .rename(columns={"variant": "region"})
+            )
+    classification = classify_reference_differences(
+        csv_results[False], csv_results[True], pd.DataFrame(output["results"])
+    )
+    record_property("csv_reference_comparisons", json.dumps(classification))
+    assert classification["widened_comparison"]["status"] == "PASS", classification
+    assert classification["unexplained_rows"] == 0, classification
+    # A legacy FAIL remains a FAIL in the diagnostic; it is not reported as parity.
     comparison = compare_frames(expected, pd.DataFrame(output["results"]))
     assert comparison["status"] == "PASS", json.dumps(
         {
diff --git a/tests/analysis/test_reference_adapters.py b/tests/analysis/test_reference_adapters.py
index d4465ea..8197344 100644
--- a/tests/analysis/test_reference_adapters.py
+++ b/tests/analysis/test_reference_adapters.py
@@ -38,3 +38,136 @@ def test_vector_adapter_only_slices_dispersion(mu):
         mu, dispersions, first_ref, first_n, rest_ref, rest_n
     )
     assert actual == expected
+
+
+@pytest.mark.parametrize("mu", [0.06, 0.25, 0.38196601125, 0.7, 0.94])
+def test_python_unphased_likelihood_preserves_archive_recurrence(mu):
+    from analysis.as_analysis import opt_unphased_dp
+
+    archive = load_archive()["as_analysis"]
+    ref, depth = np.array([10, 3, 18, 0]), np.array([20, 20, 24, 10])
+    expected = archive.opt_unphased_dp(mu, 0.03, ref[:1], depth[:1], ref[1:], depth[1:])
+    actual = opt_unphased_dp(mu, 0.03, ref[:1], depth[:1], ref[1:], depth[1:])
+    assert actual == float(np.asarray(expected).item())
+
+
+def test_python_unphased_underflow_fallback_remains_finite():
+    from scipy.stats import betabinom
+
+    from analysis.as_analysis import opt_unphased_dp
+
+    mu, rho = 0.2, 0.03
+    ref, depth = np.arange(2000, 2200), np.full(200, 100_000)
+    shape = (1 - rho) / rho
+    forward = betabinom.logpmf(ref, depth, mu * shape, (1 - mu) * shape)
+    reverse = betabinom.logpmf(ref, depth, (1 - mu) * shape, mu * shape)
+    expected = -forward[0] - np.sum(np.logaddexp(forward[1:], reverse[1:]) - np.log(2))
+    actual = opt_unphased_dp(mu, rho, ref[:1], depth[:1], ref[1:], depth[1:])
+    assert np.isfinite(actual)
+    assert np.isclose(actual, expected, rtol=5e-5, atol=5e-5)
+
+
+def test_wide_archive_accumulators_prevent_overflow(monkeypatch):
+    from scipy.sparse import csr_matrix
+    from scripts.compare_archive_analysis import single_cell_counts
+
+    archive = load_archive(widen_single_cell_sums=True)["as_analysis_sc"]
+    historical = single_cell_counts()
+    historical.layers["ref"] = csr_matrix(np.full(historical.shape, 40000, dtype=np.uint16))
+
+    def inspect_fit(function, args, **kwargs):
+        ref, depth = args
+        np.testing.assert_array_equal(ref, np.full(historical.n_obs, 160001))
+        assert ref.dtype == np.uint64
+        assert np.all(depth >= ref)
+        raise RuntimeError("inspected reference inputs")
+
+    monkeypatch.setattr(archive, "minimize_scalar", inspect_fit)
+    with pytest.raises(RuntimeError, match="inspected reference inputs"):
+        archive.get_imbalance_sc(historical, sample="sample1", groups=["A", "B"], pseudocount=1)
+
+
+@pytest.mark.parametrize("column", ["num_snps_group1", "num_snps_group2"])
+def test_differential_snp_membership_counts_require_exact_agreement(column):
+    import pandas as pd
+    from scripts.compare_archive_analysis import compare_frames
+
+    expected = pd.DataFrame({"region": ["feature"], column: [2], "fdr_pval": [0.5]})
+    observed = expected.copy()
+    observed[column] = 3
+    comparison = compare_frames(expected, observed)
+    assert comparison["status"] == "FAIL"
+    assert column in comparison["failures"]
+
+
+@pytest.mark.parametrize("column", ["ref_count", "N", "num_snps_group1", "num_snps_group2"])
+@pytest.mark.parametrize("value", [2**53, 2**63, 2**64 - 2])
+def test_count_verifier_preserves_uint64_precision(column, value):
+    import pandas as pd
+    from scripts.compare_archive_analysis import compare_frames
+
+    expected = pd.DataFrame(
+        {"region": ["feature"], column: pd.Series([value], dtype="uint64"), "fdr_pval": [0.5]}
+    )
+    observed = expected.copy()
+    assert compare_frames(expected, observed)["status"] == "PASS"
+    observed[column] = pd.Series([value + 1], dtype="uint64")
+    result = compare_frames(expected, observed)
+    assert result["status"] == "FAIL"
+    assert result["max_absolute_differences"][column] == 1
+
+
+@pytest.mark.parametrize("wide", [False, True])
+def test_bulk_csv_adapter_only_changes_count_reader_dtypes(tmp_path, monkeypatch, wide):
+    import pandas as pd
+
+    path = tmp_path / "counts.tsv"
+    pd.DataFrame(
+        {
+            "chrom": ["chr1"],
+            "pos": [101],
+            "ref": ["A"],
+            "alt": ["G"],
+            "ref_count": [40000],
+            "alt_count": [40000],
+            "other_count": [0],
+        }
+    ).to_csv(path, sep="\t", index=False)
+    module = load_archive(widen_bulk_csv_counts=wide)["as_analysis"]
+
+    def inspect_fit(frame, *args, **kwargs):
+        expected_type = np.dtype("uint64" if wide else "uint16")
+        assert frame.ref_count.dtype == frame.alt_count.dtype == frame.N.dtype == expected_type
+        assert int(frame.ref_count.iloc[0]) == 40001
+        assert int(frame.N.iloc[0]) == (80002 if wide else 14466)
+        raise RuntimeError("inspected bulk CSV fitting inputs")
+
+    monkeypatch.setattr(module, "single_model", inspect_fit)
+    with pytest.raises(RuntimeError, match="inspected bulk CSV"):
+        module.get_imbalance(str(path), pseudocount=1)
+    # Other modules are untouched by this independently selected adapter.
+    other = load_archive(widen_bulk_csv_counts=wide)["as_analysis_sc"]
+    from scripts.compare_archive_analysis import single_cell_counts
+
+    def inspect_sc(function, args, **kwargs):
+        assert all(value.dtype == np.uint16 for value in args)
+        raise RuntimeError("single-cell sums untouched")
+
+    monkeypatch.setattr(other, "minimize_scalar", inspect_sc)
+    with pytest.raises(RuntimeError, match="single-cell sums untouched"):
+        other.get_imbalance_sc(single_cell_counts(), sample="sample1", groups=["A", "B"])
+
+
+def test_reference_classification_keeps_strict_failures():
+    import pandas as pd
+    from scripts.compare_archive_analysis import classify_reference_differences
+
+    legacy = pd.DataFrame({"region": ["feature"], "mu": [0.2], "pval": [0.8], "fdr_pval": [0.8]})
+    wide = legacy.assign(mu=0.8)
+    result = classify_reference_differences(legacy, wide, wide.copy())
+    assert result["legacy_strict_comparison"]["status"] == "FAIL"
+    assert result["widened_comparison"]["status"] == "PASS"
+    assert result["unexplained_rows"] == 0
+    assert result["mismatched_rows"][0]["complementary_estimate_candidates"] == ["mu"]
+    damaged = wide.assign(pval=0.3)
+    assert classify_reference_differences(legacy, wide, damaged)["unexplained_rows"] == 1
diff --git a/tests/analysis/test_scipy_numerical_kernel.py b/tests/analysis/test_scipy_numerical_kernel.py
index 871612a..e6bc32c 100644
--- a/tests/analysis/test_scipy_numerical_kernel.py
+++ b/tests/analysis/test_scipy_numerical_kernel.py
@@ -135,3 +135,152 @@ def test_many_high_depth_unphased_observations_remain_finite(tmp_path):
     for p1, p2 in zip(np.exp(forward[1:]), np.exp(reverse[1:]), strict=True):
         product = 0.5 * (product * p1) + 0.5 * (product * p2)
     assert product == 0.0
+
+
+def decimal_beta_binomial_logpmf(k, n, alpha, beta):
+    """Independent rising-factorial PMF; no SciPy or gamma/beta function reuse."""
+    from decimal import Decimal, localcontext
+    from math import comb
+
+    with localcontext() as context:
+        context.prec = 70
+        a, b = Decimal(str(alpha)), Decimal(str(beta))
+        value = Decimal(comb(n, k))
+        for i in range(k):
+            value *= a + i
+        for i in range(n - k):
+            value *= b + i
+        for i in range(n):
+            value /= a + b + i
+        return float(value.ln())
+
+
+@pytest.mark.rust
+@pytest.mark.parametrize("ref_count,total", [(0, 0), (0, 14), (4, 8), (7, 14), (14, 14)])
+def test_likelihood_precision_against_independent_decimal(tmp_path, ref_count, total):
+    import wasp2_rust
+    from analysis.as_analysis import opt_prob
+
+    rho, alpha = 0.1, 4.5
+    expected = decimal_beta_binomial_logpmf(ref_count, total, alpha, alpha)
+    vector = opt_prob(
+        0.5, rho, np.array([ref_count], dtype=np.uint64), np.array([total], dtype=np.uint64)
+    )
+    scalar = opt_prob(0.5, rho, np.uint64(ref_count), np.uint64(total))
+    assert float(np.asarray(vector).item()) == scalar
+    assert scalar == pytest.approx(-expected, abs=1e-12, rel=1e-12)
+    path = tmp_path / "counts.tsv"
+    pd.DataFrame(
+        {
+            "chrom": ["chr1"],
+            "pos": [101],
+            "ref": ["A"],
+            "alt": ["G"],
+            "region": ["peak"],
+            "ref_count": [ref_count],
+            "alt_count": [total - ref_count],
+        }
+    ).to_csv(path, sep="\t", index=False)
+    result = wasp2_rust.analyze_imbalance_run(
+        str(path), method="single", min_count=0, pseudocount=0, region_col="region", rho=rho
+    )["results"][0]
+    assert result["null_ll"] == pytest.approx(expected, abs=1e-12, rel=1e-12)
+    if 2 * ref_count == total:
+        assert result["alt_ll"] == pytest.approx(result["null_ll"], abs=1e-12, rel=1e-12)
+        assert result["pval"] == pytest.approx(1.0, abs=1e-12)
+
+
+@pytest.mark.parametrize("first_ref,symmetric", [(3, True), (4, False)])
+def test_unphased_anchor_symmetry_is_conditional(first_ref, symmetric):
+    from analysis.as_analysis import opt_unphased_dp
+
+    rho = 0.1
+    ref = np.array([first_ref, 1, 0, 0, 0], dtype=np.uint64)
+    depth = np.array([6, 20, 0, 0, 1], dtype=np.uint64)
+    values = []
+    for mu in (0.2, 0.8):
+        expected_first = decimal_beta_binomial_logpmf(
+            int(ref[0]), int(depth[0]), mu * 9, (1 - mu) * 9
+        )
+        rest = [
+            np.logaddexp(
+                decimal_beta_binomial_logpmf(int(k), int(n), mu * 9, (1 - mu) * 9),
+                decimal_beta_binomial_logpmf(int(k), int(n), (1 - mu) * 9, mu * 9),
+            )
+            - np.log(2)
+            for k, n in zip(ref[1:], depth[1:], strict=True)
+        ]
+        expected = -expected_first - sum(rest)
+        actual = opt_unphased_dp(mu, rho, ref[:1], depth[:1], ref[1:], depth[1:])
+        assert actual == pytest.approx(expected, abs=1e-12, rel=1e-12)
+        values.append(actual)
+    if symmetric:
+        assert values[0] == pytest.approx(values[1], abs=1e-12, rel=1e-12)
+    else:
+        assert abs(values[0] - values[1]) > 1.0
+
+
+@pytest.mark.skipif(
+    (np.__version__, scipy.__version__) != ("1.26.4", "1.13.1"),
+    reason="Historical dtype dispatch reproduction uses pinned NumPy/SciPy",
+)
+def test_untouched_archive_balanced_precision_defect_is_reported(tmp_path):
+    from scripts.compare_archive_analysis import load_archive
+
+    path = tmp_path / "balanced.tsv"
+    pd.DataFrame(
+        {
+            "chrom": ["chr1"],
+            "pos": [101],
+            "ref": ["A"],
+            "alt": ["G"],
+            "region": ["peak"],
+            "ref_count": [7],
+            "alt_count": [7],
+            "other_count": [0],
+        }
+    ).to_csv(path, sep="\t", index=False)
+    narrow = load_archive()["as_analysis"].get_imbalance(
+        str(path), region_col="region", min_count=0, pseudocount=0
+    )
+    wide = load_archive(widen_bulk_csv_counts=True)["as_analysis"].get_imbalance(
+        str(path), region_col="region", min_count=0, pseudocount=0
+    )
+    assert float(narrow.pval.iloc[0]) < 0.9999
+    assert float(wide.pval.iloc[0]) == 1.0
+    assert compare_frames(narrow, wide)["status"] == "FAIL"
+
+
+@pytest.mark.rust
+@pytest.mark.parametrize("ref_count,alt_count,pc", [(40000, 40000, 1), (65536, 65537, 0)])
+def test_wide_bulk_counts_preserve_raw_outputs_and_fitting_depth(
+    tmp_path, ref_count, alt_count, pc
+):
+    import wasp2_rust
+
+    path = tmp_path / "counts.tsv"
+    pd.DataFrame(
+        {
+            "chrom": ["chr1"],
+            "pos": [101],
+            "ref": ["A"],
+            "alt": ["G"],
+            "region": ["peak"],
+            "ref_count": [ref_count],
+            "alt_count": [alt_count],
+        }
+    ).to_csv(path, sep="\t", index=False)
+    rho = 0.1
+    result = wasp2_rust.analyze_imbalance_run(
+        str(path), method="single", min_count=0, pseudocount=pc, region_col="region", rho=rho
+    )["results"][0]
+    assert (result["ref_count"], result["alt_count"], result["N"]) == (
+        ref_count,
+        alt_count,
+        ref_count + alt_count,
+    )
+    shape = 0.5 * (1 - rho) / rho
+    expected_null = betabinom.logpmf(
+        np.uint64(ref_count + pc), np.uint64(ref_count + alt_count + 2 * pc), shape, shape
+    )
+    assert result["null_ll"] == pytest.approx(expected_null, abs=1e-9, rel=1e-9)
diff --git a/tests/counting/test_single_cell_pipeline.py b/tests/counting/test_single_cell_pipeline.py
new file mode 100644
index 0000000..e072ff5
--- /dev/null
+++ b/tests/counting/test_single_cell_pipeline.py
@@ -0,0 +1,252 @@
+"""Single-cell counting and command execution against independent Python inputs."""
+
+import ast
+import contextlib
+import io
+import subprocess
+from bisect import bisect_left
+from collections import defaultdict
+from pathlib import Path
+
+import anndata as ad
+import numpy as np
+import pandas as pd
+import polars as pl
+import pysam
+import pytest
+from scripts.compare_archive_analysis import compare_frames, load_archive, single_cell_counts
+from typer.testing import CliRunner
+
+from analysis.__main__ import app as analysis_app
+from counting.__main__ import app as count_app
+from counting.count_alleles_sc import CountStatsSC, _sparse_from_counts, count_bc_snp_alleles
+from counting.run_counting_sc import WaspCountSC
+
+ROOT = Path(__file__).resolve().parents[2]
+ARCHIVE = "41354fc087404ba1fa99542803ee9cb9b7e0e6bf"
+
+
+def archive_counter():
+    namespace = {"bisect_left": bisect_left}
+    for path, name in [
+        ("src/counting/count_alleles.py", "find_read_aln_pos"),
+        ("src/counting/count_alleles_sc.py", "count_bc_snp_alleles"),
+    ]:
+        source = subprocess.check_output(["git", "show", f"{ARCHIVE}:{path}"], cwd=ROOT)
+        tree = ast.parse(source)
+        function = next(n for n in tree.body if isinstance(n, ast.FunctionDef) and n.name == name)
+        exec(compile(ast.Module(body=[function], type_ignores=[]), path, "exec"), namespace)
+    return namespace["count_bc_snp_alleles"]
+
+
+@pytest.mark.rust
+@pytest.mark.parametrize(
+    "scenario", ["overlap", "paired", "deletion", "skip", "barcodes", "missing_sequence"]
+)
+@pytest.mark.parametrize("reverse", [False, True])
+@pytest.mark.parametrize("threads", [1, 4])
+def test_three_way_single_cell_counting(tmp_path, scenario, reverse, threads):
+    import wasp2_rust
+
+    bam_path = tmp_path / "reads.bam"
+    header = {
+        "HD": {"VN": "1.6", "SO": "coordinate"},
+        "SQ": [{"SN": c, "LN": 1000} for c in ["chr1", "chr2"]],
+    }
+    with pysam.AlignmentFile(str(bam_path), "wb", header=header) as out:
+        for chrom in [0, 1]:
+            for i in range(6):
+                read = pysam.AlignedSegment(out.header)
+                read.query_name = "pair" if scenario == "paired" and i < 2 else f"read{i}"
+                read.reference_id, read.reference_start = chrom, 90 + i
+                read.mapping_quality, read.flag = 60, 0
+                if scenario == "paired" and i < 2:
+                    read.flag = 65 if i == 0 else 129
+                    read.next_reference_id = chrom
+                    read.next_reference_start = 91 if i == 0 else 90
+                read.cigarstring = {"deletion": "10M1D10M", "skip": "10M1N10M"}.get(scenario, "30M")
+                length = 20 if scenario in ["deletion", "skip"] else 30
+                read.query_sequence = ("A" if i % 2 == 0 else "G") * length
+                if scenario == "missing_sequence" and i == 0:
+                    read.query_sequence = None
+                if scenario != "barcodes" or i != 0:
+                    read.set_tag(
+                        "CB", "unknown" if scenario == "barcodes" and i == 1 else ["A", "B"][i % 2]
+                    )
+                out.write(read)
+    pysam.index(str(bam_path))
+    snps = [(0, 101, "A", "G"), (1, 111, "A", "G")]
+    if reverse:
+        snps.reverse()
+    barcodes = {"A": 0, "B": 1}
+    expected = []
+    for counter in [archive_counter(), count_bc_snp_alleles]:
+        counts = CountStatsSC()
+        with pysam.AlignmentFile(str(bam_path), "rb") as bam:
+            for chrom_index, chrom in enumerate(["chr1", "chr2"]):
+                chromosome_snps = [
+                    (idx + 2 * chrom_index, pos, ref, alt) for idx, pos, ref, alt in snps
+                ]
+                counter(bam, barcodes, chrom, iter(chromosome_snps), counts)
+        expected.append([dict(counts.ref_count), dict(counts.alt_count), dict(counts.other_count)])
+    actual = [{}, {}, {}]
+    regions = [
+        (idx + 2 * chrom_index, chrom, pos, ref, alt)
+        for chrom_index, chrom in enumerate(["chr1", "chr2"])
+        for idx, pos, ref, alt in snps
+    ]
+    for snp, bc, ref, alt, other in wasp2_rust.BamCounterSC(str(bam_path)).count_bc_alleles(
+        regions, barcodes, 0, threads
+    ):
+        for layer, value in zip(actual, [ref, alt, other], strict=True):
+            if value:
+                layer[snp, bc] = value
+    assert expected[0] == expected[1] == actual
+
+
+@pytest.mark.parametrize("value", [0, 65535, 65536, 70000, 2**32])
+def test_sparse_counts_preserve_large_values_through_h5ad(tmp_path, value):
+    counts = defaultdict(int, {(0, 0): value}) if value else defaultdict(int)
+    matrix = _sparse_from_counts(counts, (1, 1))
+    adata = ad.AnnData(X=matrix + matrix, layers={"ref": matrix, "alt": matrix})
+    path = tmp_path / "counts.h5ad"
+    adata.write_h5ad(path)
+    restored = ad.read_h5ad(path)
+    assert restored.layers["ref"][0, 0] == value
+    assert restored.X[0, 0] == 2 * value
+
+
+@pytest.mark.parametrize("form", ["string", "list", "mixed", "file", "file_list"])
+def test_sample_arguments_normalize_consistently(tmp_path, form):
+    path = tmp_path / "samples.txt"
+    path.write_text("sample1\n sample2 \n\n")
+    samples = {
+        "string": "sample1, sample2",
+        "list": ["sample1", "sample2"],
+        "mixed": ["sample1,sample2"],
+        "file": str(path),
+        "file_list": [str(path)],
+    }[form]
+    options = WaspCountSC("reads.bam", "variants.vcf", "barcodes.tsv", None, samples=samples)
+    assert options.samples == ["sample1", "sample2"]
+
+
+@pytest.fixture
+def synthetic_count_bundle(tmp_path):
+    historical = single_cell_counts()
+    positions = [100 + 100 * i for i in range(6)]
+    bam = tmp_path / "reads.bam"
+    with pysam.AlignmentFile(
+        str(bam),
+        "wb",
+        header={"HD": {"VN": "1.6", "SO": "coordinate"}, "SQ": [{"SN": "chr1", "LN": 1000}]},
+    ) as out:
+        for snp, pos in enumerate(positions):
+            for cell in range(4):
+                for allele, layer in [("A", "ref"), ("G", "alt")]:
+                    for k in range(int(historical.layers[layer][snp, cell])):
+                        read = pysam.AlignedSegment(out.header)
+                        read.query_name = f"{snp}_{cell}_{allele}_{k}"
+                        read.query_sequence, read.flag = allele, 0
+                        read.reference_id, read.reference_start = 0, pos - 1
+                        read.mapping_quality, read.cigarstring = 60, "1M"
+                        read.set_tag("CB", f"cell{cell}")
+                        out.write(read)
+    pysam.index(str(bam))
+    vcf = tmp_path / "variants.vcf"
+    vcf.write_text(
+        '##fileformat=VCFv4.2\n##contig=<ID=chr1,length=1000>\n##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tsample1\n'
+        + "".join(
+            f"chr1\t{pos}\t.\tA\tG\t60\tPASS\t.\tGT\t" + ["0|1", "1|0"][i % 2] + "\n"
+            for i, pos in enumerate(positions)
+        )
+    )
+    barcodes, groups, features = [
+        tmp_path / name for name in ["barcodes.tsv", "groups.tsv", "features.bed"]
+    ]
+    barcodes.write_text("".join(f"cell{i}\n" for i in range(4)))
+    groups.write_text("cell0\tA\ncell1\tA\ncell2\tB\ncell3\tB\n")
+    features.write_text("chr1\t90\t210\tp0\nchr1\t290\t310\tp1\nchr1\t390\t610\tp2\n")
+    return bam, vcf, barcodes, groups, features, historical
+
+
+@pytest.mark.rust
+@pytest.mark.parametrize("command", ["find-imbalance-sc", "compare-imbalance"])
+@pytest.mark.parametrize("phased", [False, True])
+@pytest.mark.parametrize("pc", [0, 1])
+def test_bam_to_h5ad_to_analysis_matches_archive(
+    tmp_path, synthetic_count_bundle, monkeypatch, command, phased, pc
+):
+    import shutil
+
+    if not shutil.which("bedtools"):
+        pytest.skip("bedtools required for feature intersection")
+    bam, vcf, barcodes, groups, features, historical = synthetic_count_bundle
+    # Metadata conversion must work without the undeclared optional PyArrow package.
+    monkeypatch.setattr(
+        pl.DataFrame, "to_pandas", lambda *a, **k: pytest.fail("PyArrow-dependent conversion")
+    )
+    counts = tmp_path / "counts.h5ad"
+    result = CliRunner().invoke(
+        count_app,
+        [
+            "count-variants-sc",
+            str(bam),
+            str(vcf),
+            str(barcodes),
+            "--samples",
+            "sample1",
+            "--feature",
+            str(features),
+            "-o",
+            str(counts),
+        ],
+    )
+    assert result.exit_code == 0, (result.output, result.exception)
+    measured = ad.read_h5ad(counts)
+    assert measured.shape == (4, 6)
+    assert measured.obs_names.tolist() == [f"cell{i}" for i in range(4)]
+    for layer in ["ref", "alt"]:
+        np.testing.assert_array_equal(
+            measured.layers[layer].toarray(), historical.layers[layer].toarray().T
+        )
+    output = tmp_path / "results.tsv"
+    result = CliRunner().invoke(
+        analysis_app,
+        [
+            command,
+            str(counts),
+            str(groups),
+            "--sample",
+            "sample1",
+            "--phased" if phased else "--unphased",
+            "--pseudocount",
+            str(pc),
+            "-o",
+            str(output),
+        ],
+    )
+    assert result.exit_code == 0, (result.output, result.exception)
+    archive = load_archive()
+    function = (
+        archive["as_analysis_sc"].get_imbalance_sc
+        if command == "find-imbalance-sc"
+        else archive["compare_ai"].get_compared_imbalance
+    )
+    reference = measured.T.copy()
+    reference.var["group"] = ["A", "A", "B", "B"]
+    with contextlib.redirect_stdout(io.StringIO()):
+        expected = function(
+            reference,
+            sample="sample1",
+            groups=["A", "B"],
+            phased=phased,
+            pseudocount=pc,
+            min_count=10,
+        )
+    assert expected
+    for key, frame in expected.items():
+        suffix = key if isinstance(key, str) else "_".join(key)
+        observed = pd.read_csv(tmp_path / f"results_{suffix}.tsv", sep="\t")
+        comparison = compare_frames(frame, observed)
+        assert comparison["status"] == "PASS", comparison