# Reproduction plan — Moncla H3Nx host-strategy paper

Goal: use the `*eon` suite to independently pressure-test the two critiques from `PEER_REVIEW.md`
that are directly answerable with our tools. This is a **staging skeleton** — Steven takes it from here.

Env: `PY=/home/sweaver/.conda/envs/axomeme/bin/python` · eon CLIs `hyphaeon`, `chronaeon` on PATH.
GPU inference should run on SLURM `datamonkey` partition (only one with /storage), per axomeme CLAUDE.md.

## Data (see `scripts/00_fetch_data.sh`)
Source of truth is the authors' repo `github.com/moncla-lab/h3nx-paper` (alignments, trees, host
labels, TreeSort outputs) + the GISAID acknowledgment table (Supp Table 2) for isolate provenance.
Human seasonal H3N2 and non-human H3Nx were subsampled differently (see paper Methods) — preserve the
authors' per-host subtree alignments/trees rather than re-subsampling, so we test THEIR inference, not a
new dataset. Target layout after fetch:
```
data/
  ha/  {north_am_avian,eurasian_avian,swine,human,equine,canine}.{fasta,nwk}
  na/  { ... same host splits ... }
  pb1/ { ... conserved-gene control ... }
  treesort/  per-host summary reassortment trees + support JSON
```

## R1 — M2: is "no avian adaptation" real, or a persistence-dependent-estimator artifact?
The paper's MK/Bhatt estimator needs temporal persistence; avian lineages turn over fast, so it can
read ≈0 even if episodic positive selection is present. `hyphaeon meme` does **site-level episodic
diversifying selection** and does NOT require temporal persistence — so it can distinguish "no selection"
from "no signal the slope-estimator can integrate."
```bash
# per host × gene (HA, NA; PB1 as the same conserved control the paper uses)
for gene in ha na pb1; do
  for host in north_am_avian eurasian_avian swine human equine canine; do
    hyphaeon meme -a data/$gene/$host.fasta -t data/$gene/$host.nwk \
      --filter -o results/meme_${gene}_${host}.json -c results/meme_${gene}_${host}.csv
  done
done
```
**Prediction under the paper's thesis:** avian HA/NA show few/no MEME-significant sites (matches "no adaptation").
**Prediction under M2:** avian HA/NA show episodic-selection sites comparable to mammals → the MK ≈0 was an
estimator-persistence artifact, and the "no directional selection in birds" claim needs softening.
Compare site counts (p≤0.05) and LRT distributions across hosts; PB1 is the negative control (should be low everywhere).

## R2 — M4/m5: independently reproduce host-specific molecular clocks
Reassortment *rates* are scaled by per-host clock rates (Fig 2 legend). If our clocks disagree with the
paper's TreeTime clocks, the avian:swine rate ratio (and thus "reassortment-dominant in birds") shifts.
```bash
for host in north_am_avian eurasian_avian swine human equine canine; do
  chronaeon date -a data/ha/$host.fasta -t data/ha/$host.nwk \
    --date-regex '\d{4}(-\d{2}-\d{2})?' --method all --clock-model auto \
    --bootstrap 1000 --ci-method residual-boot --loocv \
    -o results/clock_ha_$host.json -c results/clock_ha_$host.csv
done
# cross-check unsupervised multi-clock structure per host (are there hidden sub-clocks the single-rate scaling misses?)
chronaeon autoclock -a data/ha/north_am_avian.fasta -t data/ha/north_am_avian.nwk \
  -o results/autoclock_na_avian.json
```
Compare `clock_rate ± CI` to the paper's TreeTime per-host rates; propagate any delta into the reassortment-rate
comparison to see if the rank ordering (avian > swine) survives.

## Not staged here (need author artifacts / out of eon scope)
- M1 (recompute adaptive rates *including* reassortant strains): needs the authors' MK pipeline
  (`github.com/blab/adaptive-evolution`) + their reassortant exclusion list; wire up after data fetch.
- M3 (rate-matched persistence): operate on `data/treesort/` summary trees — a small script to downsample
  avian reassortment events to the swine rate then recompute persistence. Stub in `scripts/` later.

## Status
- [x] review subdir + PEER_REVIEW.md staged
- [x] reproduction skeleton + fetch script staged
- [x] GISAID data landed (data/{ha,na,pb1}/*.fasta+nwk); NA reframed to codon frame (scripts/01_reframe_na.py)
- [x] R1 first pass (array 1337086) — INVALID: HA analyzed in wrong frame (frame 0). Retracted.
- [x] **R1 corrected (array 1337114, HA reframed to frame 2)** — avian HA = 0 FDR sites; **M2 NOT supported**,
      critique withdrawn. Our --no-tree HA test is likely underpowered (≤1 site even human/swine). See `RESULTS.md`.
- [x] R2 single-clock (array 1337107) — INVALID: tMRCA 1332 CE nonsense on pooled avian sublineages.
- [x] **R2b chronaeon autoclock (array 1337115)** — multi-clock deconvolution; rates now sane; **M4 worry not
      supported but test inconclusive** (single-gene HA / TN93, not paper's per-segment TreeTime clocks).
- [ ] R3 rate-matched persistence (M3) — needs data/treesort/ (still empty) — SKIPPED per Steven
- [x] **R1c tree-based rerun (array 1337127)** — authors' topology pruned to alignment tips. Tree ≈ TN93
      (avian HA still 0 FDR sites); power question resolved, near-zero avian result is REAL. M2 definitively withdrawn.
- [ ] (would-be-decisive) R2 rerun on the paper's reference-segment TreeTime clocks + matched trees for M4

## Lessons (do not rediscover)
1. **Frame-check every alignment before MEME**, not just length%3. These GISAID HA & NA alignments both carry a
   2-column frame-2 offset; PB1 is frame-0. Translate the longest ungapped seq in all 3 frames and pick the one
   with ~0 internal stops (HA/NA N-termini: MKTII... / MNPNQK...). A mult3 alignment can still be out of frame.
2. **For clock estimates on diverse multi-lineage host sets, use `chronaeon autoclock`, not `chronaeon date`.**
   Single-clock pools divergent sublineages and produces garbage (tMRCA 1332 CE). autoclock deconvolves into
   clock communities with sane rates.
3. **Compute-node env:** system-py3.9 user-site needed pytz, packaging, pyparsing, cycler, kiwisolver, fonttools
   installed into ~/.local (shared fs) for the eon tools to import under srun on datamonkey nodes.

## Full audit (harness h3nx_audit_harness.js, run w0d6jmd4q) — 35 agents, 51 claims, 24 upheld findings
Report: ../FINAL_AUDIT_REPORT.md . Verdict: **accept with MINOR revision** — thesis survived hostile reproduction
against the authors' own code; every conclusion-undermining hypothesis fell in the paper's favor.
Open follow-ups the audit surfaced (reproducibility, not correctness):
- [ ] Segment-enrichment Monte-Carlo null (C24-C26, load-bearing) — generating code ABSENT from repo (verified:
      only flyway_binomial is present). Ask authors to archive it; attempt independent recompute from summary.nwk.
- [ ] Eurasian-avian (C14=0.1248) + canine (C17) reassortment rates — no shipped per-replicate log.csv; recompute
      from summary.nwk × clock rate to confirm.
- [ ] p-value transcription: paper text "9 x 10-7" for host-switch OR should be 9.34e-8 (verified vs authors' FET).
- [ ] C6 HA tMRCA CI text (1910-1917) vs shipped auspice JSON ([1908.85, 1916.20]) mismatch — reconcile.