Add donor-local dispersion while preserving Python imbalance tests

.github/workflows/ci.yml

@@ -66,6 +66,8 @@ jobs:
       - uses: actions/checkout@de0fac2e4500dabe0009e67214ff5f5447ce83dd  # v6.0.2
         with:
           persist-credentials: false
+          # Runtime reference tests read the pinned pre-Rust archive commit.
+          fetch-depth: 0
 
       - name: Install system dependencies
         run: |
@@ -114,6 +116,14 @@ jobs:
           source .venv/bin/activate
           python -c "import wasp2_rust; print('wasp2_rust loaded:', wasp2_rust.__file__)"
 
+      - name: Pin historical numerical reference for the full real ATAC matrix
+        if: matrix.python-version == '3.11'
+        run: uv pip install --python .venv/bin/python numpy==1.26.4 scipy==1.13.1
+
+      - name: Diagnose platform numerical arithmetic
+        if: matrix.python-version == '3.11'
+        run: .venv/bin/python scripts/diagnose_bulk_numerics.py
+
       - name: Run tests
         run: |
           source .venv/bin/activate

.github/workflows/docker.yml

@@ -19,32 +19,49 @@ permissions:
 jobs:
   build:
     name: Container ${{ matrix.platform }}
-    runs-on: ubuntu-latest
+    runs-on: ${{ matrix.runner }}
     timeout-minutes: 60
     strategy:
       fail-fast: false
       matrix:
-        platform: [linux/amd64, linux/arm64]
+        include:
+          - {runner: ubuntu-24.04, platform: linux/amd64}
+          - {runner: ubuntu-24.04-arm, platform: linux/arm64}
     steps:
       - uses: actions/checkout@de0fac2e4500dabe0009e67214ff5f5447ce83dd  # v6.0.2
         with:
           persist-credentials: false
       - uses: docker/setup-qemu-action@29109295f81e9208d7d86ff1c6c12d2833863392  # v3
       - uses: docker/setup-buildx-action@8d2750c68a42422c14e847fe6c8ac0403b4cbd6f  # v3
+      - name: Read package version
+        id: version
+        run: |
+          python3 - <<'PY' >> "$GITHUB_OUTPUT"
+          import tomllib
+          with open("rust/Cargo.toml", "rb") as handle:
+              print("value=" + tomllib.load(handle)["package"]["version"])
+          PY
       - name: Build local test image
         uses: docker/build-push-action@10e90e3645eae34f1e60eeb005ba3a3d33f178e8  # v6
         with:
           context: .
           load: true
           platforms: ${{ matrix.platform }}
           tags: wasp2:ci
-          build-args: VERSION=0+ci.${{ github.sha }}
+          build-args: |
+            VERSION=${{ steps.version.outputs.value }}
+            REVISION=${{ github.sha }}
           cache-from: type=gha
           cache-to: type=gha,mode=max
       - name: Smoke-test image
         env:
           PLATFORM: ${{ matrix.platform }}
-        run: docker run --rm --platform "$PLATFORM" wasp2:ci /opt/wasp2/scripts/container_smoke_test.sh
+          EXPECTED_VERSION: ${{ steps.version.outputs.value }}
+          EXPECTED_REVISION: ${{ github.sha }}
+        run: |
+          test "$(docker inspect --format '{{ index .Config.Labels "org.opencontainers.image.version" }}' wasp2:ci)" = "$EXPECTED_VERSION"
+          test "$(docker inspect --format '{{ index .Config.Labels "org.opencontainers.image.revision" }}' wasp2:ci)" = "$EXPECTED_REVISION"
+          docker run --rm --platform "$PLATFORM" -e EXPECTED_WASP2_VERSION="$EXPECTED_VERSION" wasp2:ci /opt/wasp2/scripts/container_smoke_test.sh
 
   required:
     name: Required / Container

.github/workflows/release.yml

@@ -281,14 +281,21 @@ jobs:
           load: true
           platforms: ${{ matrix.platform }}
           tags: wasp2:release-smoke-${{ matrix.tag }}
-          build-args: VERSION=${{ needs.validate.outputs.version }}
+          build-args: |
+            VERSION=${{ needs.validate.outputs.version }}
+            REVISION=${{ github.sha }}
           cache-from: type=gha,scope=release-${{ matrix.tag }}
           cache-to: type=gha,mode=max,scope=release-${{ matrix.tag }}
       - name: Smoke-test local release image
         env:
           IMAGE: wasp2:release-smoke-${{ matrix.tag }}
           PLATFORM: ${{ matrix.platform }}
-        run: docker run --rm --platform "$PLATFORM" "$IMAGE" /opt/wasp2/scripts/container_smoke_test.sh
+          EXPECTED_VERSION: ${{ needs.validate.outputs.version }}
+          EXPECTED_REVISION: ${{ github.sha }}
+        run: |
+          test "$(docker inspect --format '{{ index .Config.Labels "org.opencontainers.image.version" }}' "$IMAGE")" = "$EXPECTED_VERSION"
+          test "$(docker inspect --format '{{ index .Config.Labels "org.opencontainers.image.revision" }}' "$IMAGE")" = "$EXPECTED_REVISION"
+          docker run --rm --platform "$PLATFORM" -e EXPECTED_WASP2_VERSION="$EXPECTED_VERSION" "$IMAGE" /opt/wasp2/scripts/container_smoke_test.sh
 
   singularity-smoke:
     name: Singularity smoke
@@ -311,7 +318,9 @@ jobs:
           context: .
           platforms: linux/amd64
           tags: wasp2:release-singularity
-          build-args: VERSION=${{ needs.validate.outputs.version }}
+          build-args: |
+            VERSION=${{ needs.validate.outputs.version }}
+            REVISION=${{ github.sha }}
           outputs: type=docker,dest=${{ runner.temp }}/wasp2-docker.tar
           cache-from: type=gha,scope=release-amd64
           cache-to: type=gha,mode=max,scope=release-amd64
@@ -373,7 +382,9 @@ jobs:
           tags: |
             ${{ env.IMAGE_NAME }}:${{ needs.validate.outputs.version }}
             ${{ env.IMAGE_NAME }}:latest
-          build-args: VERSION=${{ needs.validate.outputs.version }}
+          build-args: |
+            VERSION=${{ needs.validate.outputs.version }}
+            REVISION=${{ github.sha }}
           provenance: mode=max
           sbom: true
       - name: Smoke-test published architectures

CHANGELOG.md

@@ -5,37 +5,81 @@ All notable changes to WASP2 will be documented in this file.
 ## [Unreleased]
 
 ### Fixed
+- 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.
+- Cephes arithmetic follows the measured reference compiler policy: separate
+  multiplication/addition for Linux x86_64 and fused operations for the macOS
+  arm64 reference. Platform probes distinguish this from statistical-model changes.
+- Bulk phasing retains the original genotype tokens for archived table-wide
+  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
+  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
+  labels; smoke tests require a matching Python/Rust distribution version.
+- CLI help now reports the actual donor-local pseudocount default of 1.
+- Bulk beta-binomial calculations use a SciPy 1.13.1-compatible Cephes kernel,
+  reference FMA arithmetic, and NumPy float64 reductions. Unphased inference
+  follows archived probability accumulation when safe and uses log-space fallback
+  when products lose range. Linear fits preserve adjusted-depth coefficients in
+  the existing four-field representation to avoid conversion-rounding changes.
+  Strict real ATAC comparisons include raw allelic proportions.
+- Bulk archive validation separately verifies symmetric unphased optima using
+  retained anchor counts and independent Python objective evaluations, while
+  preserving strict comparisons, raw proportions, and numerical tolerances.
+- Allelic-proportion optimization uses archived Python bounded Brent defaults
+  for single-SNP, phased, and unphased bulk tests, including empty-allele inputs
+  with pseudocount 0. Balanced first-SNP unphased objectives can still select
+  mirrored proportions across numerical libraries; no tie-breaking rule is added.
+- Linear dispersion follows archived Python single-start Nelder-Mead fitting
+  on pseudocounted depth, retains coefficients with a pseudocount depth-origin
+  shift for exact reusable prediction. Equal depths use a unit scale; testing no longer applies
+  the fitting-only predictor clip.
+- Retain distinct repeated SNP count observations in donor analysis and apply
+  Python phase-aware deduplication after depth filtering. SNV outputs aggregate
+  retained observations per position; fitting preserves their phase identities.
+- Donor feature analysis follows archived Python phasing fallback for missing,
+  identical, or nonnumeric GT after depth filtering, without converting allele
+  strings. Provenance records requested and effective phase per included donor.
+- Donor dispersion fits retain overlapping SNP-feature observations as archived
+  Python does, with separate donor-feature identities in global fits. Unique-SNP
+  QC and donor eligibility counts retain their existing meanings.
 - Bulk SNP counting preserves archived Python input-order priority across paired
   mates, with read names consumed before checking aligned bases at each SNP.
-- Bulk region outputs remove one pseudocount per retained SNP and allele, restoring
-  the archived Python raw-count aggregation without changing the fitted statistics.
-- Single-cell analysis input preparation and differential comparisons now follow
-  the cells-by-variants AnnData layout throughout, including barcode groups and genotypes.
-- Differential pseudobulk sums widen to UInt64 rather than wrapping at 65,535.
+- Restore archived chi-square LRT p-values and first-SNP-anchored unphased
+  feature likelihoods in donor-local and legacy bulk analysis.
+- Default donor-local pseudocounts to 1, matching archived Python. Per-donor
+  nuisance fitting and independent donor tests/BH correction remain available.
+- Integrate the raw-count and single-cell layout compatibility fixes from #135.
+  Exact-SNV requests are rejected rather than switching the p-value method.
 
-### Tests
-- Frozen outputs from archived Python commit `41354fc` cover bulk, single-cell,
-  and differential analysis. Linear-reference array-shape adaptations are explicit;
-  the archive's ignored phased-linear flag is reported as a reference limitation.
-
-## [1.5.0] - 2026-07-14
+## [1.5.0] - 2026-07-21
 
 ### Added
 - `wasp2-count count-cohort` for one-final-BAM-per-donor bulk counting with explicit
   `donor_id`/`vcf_sample` mapping, atomic locked outputs, and SHA-256 provenance manifests.
 - Linear-dispersion and donor-specific-dispersion beta-binomial implementations in the Rust
   analysis API, including phased feature-level inference.
+- Donor-local bulk imbalance orchestration through `wasp2-analyze find-imbalance` with explicit
+  SNV or feature units, global or per-donor nuisance fitting, donor-local inference and BH
+  correction, manifest pinning, and non-overwriting result, dispersion, QC, and provenance files.
 - `--per-variant` / `--snv-solo` for independent SNV tests and explicit `--region_col` support
   for feature- or peak-level aggregation.
 - A Rust single-cell ATAC allele counter with chromosome-parallel execution and Python fallback.
 
 ### Changed
 - VCF conversion now selects biallelic SNVs by default; multi-allelic records require the
   explicit `--include-multiallelic` option.
-- Phased inference accepts only numeric VCF heterozygote orientation (`0|1` and `1|0`) for parity
-  with the archived Python implementation; allele-string genotypes do not imply cross-SNV phase.
+- Phased inference follows archived numeric first-allele parsing after table-wide
+  genotype checks; allele-string genotypes do not imply cross-SNV phase.
 - `--phased` and `--region_col` are forwarded to the Rust backend; unsupported `--groupby`
   requests fail closed and direct users to `--region_col <parent_column>`.
+- SNV analysis is unphased by contract, while feature analysis supports phased and unphased
+  likelihoods without pooling effect inference across donors.
 - Single-cell AnnData output now follows the scverse `obs=cells`, `var=variants` orientation.
 - Packaging derives the Python version from installed metadata and synchronizes Cargo, pixi,
   containers, Galaxy, Nextflow, documentation, and citation metadata from one release command.
@@ -48,6 +92,10 @@ All notable changes to WASP2 will be documented in this file.
 - Unphased multi-SNV likelihoods now use log-space dynamic programming to prevent underflow.
 - Linear inference uses per-SNV dispersion throughout the alternative optimization, phased
   donor-specific inference uses each observation's donor dispersion, and LRTs are non-negative.
+- Linear nuisance fits now standardize raw depth and persist the fitted center and scale with both
+  coefficients so fixed global fits are reused exactly in every donor-local run.
+- Ambiguous multi-allelic count rows are rejected instead of being silently coerced into
+  biallelic SNV identities.
 - Bulk counts use `UInt32`, and single-cell pseudobulk/observation sums widen safely instead of
   overflowing at 65,535; single-cell counting also works without an optional feature file.
 - Count-analysis TSV parsing is header-driven, so `pos0`, `GT`, sample, and donor metadata cannot

CITATION.cff

@@ -23,7 +23,7 @@ authors:
     affiliation: "Salk Institute for Biological Studies"
 
 version: "1.5.0"
-date-released: "2026-07-14"
+date-released: "2026-07-21"
 license: MIT
 
 repository-code: "https://github.com/mcvickerlab/WASP2"

Dockerfile

@@ -70,13 +70,15 @@ FROM python:3.11-slim-bookworm
 # Version: keep in sync with rust/Cargo.toml (single source of truth)
 # Run scripts/check-version-consistency.sh to verify
 ARG VERSION=1.5.0
+ARG REVISION=unknown
 
 LABEL org.opencontainers.image.source="https://github.com/mcvickerlab/WASP2"
 LABEL org.opencontainers.image.description="WASP2: Allele-specific analysis of NGS data with Rust acceleration"
 LABEL org.opencontainers.image.licenses="MIT"
 LABEL org.opencontainers.image.vendor="mcvickerlab"
 LABEL org.opencontainers.image.title="WASP2"
 LABEL org.opencontainers.image.version="${VERSION}"
+LABEL org.opencontainers.image.revision="${REVISION}"
 LABEL maintainer="Jeff Jaureguy <jeffpjaureguy@gmail.com>"
 
 # Install runtime deps + temporary build deps for pybedtools C++ extension

README.md

@@ -114,6 +114,77 @@ to retain peak identifiers for feature-level analysis.
 wasp2-analyze find-imbalance counts.tsv -o ai_results.tsv
 ```
 
+### Donor-local cohort analysis
+
+Analyze the locked count directory (or its `count_manifest.json` or
+`counts.tsv.gz`) with an explicit unit. SNV analysis rejects `--phased` and tests
+each variant independently within each donor:
+
+```bash
+wasp2-analyze find-imbalance cohort_counts \
+  --scope per-donor --unit snv --min-donor-observations 50 \
+  -o atac_snv_results.tsv
+```
+
+`--dispersion-scope per-donor` is the default. `--dispersion-scope global`
+fits one nuisance model over retained eligible donor-feature-SNV observations and reuses that exact
+fit for every included donor; testing and BH correction remain donor-local. Both
+`single` and `linear` dispersion models are supported. Linear fits record and
+reuse `linear_d1`, `linear_d2`, `linear_depth_center`, and
+`linear_depth_scale` together.
+New fits preserve the original adjusted-depth coefficients: `linear_depth_center`
+is `-2 * pseudocount` and `linear_depth_scale` is 1. This avoids rounding from
+converting the predictor to mean/std coordinates. Existing fixed parameter sets
+continue to use their supplied center and scale.
+
+Allelic-proportion fitting follows archived Python bounded Brent defaults for
+SNP, phased, and unphased bulk tests. An unphased region whose first retained SNP
+has equal reference and alternate counts has a symmetric likelihood: `mu` and
+`1 - mu` are equally likely. Numerical libraries may select opposite optima;
+the reported direction is not identifiable in that case.
+
+Archive validation reports strict numerical matching separately from statistical
+equivalence. A mirrored unphased proportion qualifies only with a balanced first
+retained SNP, complementary proportions, independently verified equal archived
+likelihoods, and passing counts, dispersion and all other statistical checks.
+Reports retain the raw proportion differences and strict failures.
+
+Feature analysis preserves one SNV in every overlapping feature. `peak` remains
+an alias for `feature`:
+
+```bash
+wasp2-count count-cohort donors.tsv variants.vcf.gz feature_counts \
+  --unit feature --regions peaks.bed
+
+wasp2-analyze find-imbalance feature_counts \
+  --scope per-donor --unit feature --model linear --phased \
+  -o atac_feature_results.tsv
+```
+
+The donor-local route defaults to pseudocount 1. It writes
+`atac_feature_results.tsv`, `atac_feature_results.dispersion.tsv`,
+`atac_feature_results.qc.tsv`, and `atac_feature_results.provenance.json` and
+refuses to overwrite any of them. For standalone legacy count tables only,
+`sample` is accepted and normalized to `donor_id`; providing both columns is an
+error. Use `--expected-manifest-sha256` to require an externally trusted locked
+bundle manifest. An input `snv_id` is preserved exactly after validating its
+donor and allele identity; coordinate-based IDs are synthesized only when the
+column is absent. Dispersion and provenance metadata report optimizer status
+only when Rust provides it; `parameters_returned` does not assert convergence.
+
+Phasing checks run after depth filtering. Missing GT, identical GT tokens, or
+absence of numeric `0|1`/`1|0` genotypes selects unphased inference. Otherwise,
+every retained GT must have a numeric first allele: mixed `0|0`/`1|1` tokens
+follow archived Python's first-allele behavior, while malformed mixed tokens
+(such as `0/1` or `A|G`) raise an error. For phased inference, prepare a table of
+truly phased heterozygous SNPs; otherwise omit `--phased`.
+
+Valid interior dispersion values are retained without clipping. Exact rho
+endpoints and probability-product underflow use documented finite extensions
+because archived Python is undefined there. The full public ATAC comparison
+and its reference adaptations are described in
+[the validation report](docs/BULK_PARITY_VALIDATION.md).
+
 ## Single-Cell Example
 
 ```bash