diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..5ef958c --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,88 @@ +name: CI + +on: + push: + branches: [main, refactor/main] + pull_request: + workflow_dispatch: # allow manual runs from the Actions tab (any branch) + +concurrency: + group: ci-${{ github.ref }} + cancel-in-progress: true + +jobs: + test: + name: Unit + regression tests + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: "3.11" + cache: pip + + - name: Cache external tools (PLINK/PLINK2/KING) + uses: actions/cache@v4 + with: + path: ~/.genotools + key: genotools-deps-${{ runner.os }}-${{ hashFiles('genotools/dependencies.py') }} + + - name: Install genotools + test deps + run: | + python -m pip install --upgrade pip + pip install -e . + pip install pytest + + - name: Download PLINK / PLINK2 + # The genotools dependency resolver fetches the Linux binaries into + # ~/.genotools/misc/executables so the QC/GWAS steps (and the parity + # genotype comparator) can find them. + run: | + python -c "from genotools.dependencies import check_plink, check_plink2; print(check_plink()); print(check_plink2())" + + - name: Run unit + regression tests + # Parity tests under tests/regression skip cleanly here (no .venv-stable); + # they run in the dedicated parity job below. + run: pytest tests/unit tests/regression -q + + parity: + name: Old-vs-new parity + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + with: + # Full history so the pre-refactor baseline ref (origin/main) is + # available to build .venv-stable from. + fetch-depth: 0 + + - uses: actions/setup-python@v5 + with: + python-version: "3.11" + cache: pip + + - name: Cache external tools (PLINK/PLINK2/KING) + uses: actions/cache@v4 + with: + path: ~/.genotools + key: genotools-deps-${{ runner.os }}-${{ hashFiles('genotools/dependencies.py') }} + + - name: Install genotools (new) + test deps + run: | + python -m pip install --upgrade pip + pip install -e . + pip install pytest + + - name: Download PLINK / PLINK2 + run: | + python -c "from genotools.dependencies import check_plink, check_plink2; print(check_plink()); print(check_plink2())" + + - name: Fetch pre-refactor baseline ref + run: git fetch origin main:refs/remotes/origin/main + + - name: Build .venv-stable (pre-refactor baseline) + # umap-learn==0.5.3 (old dep) needs Python <=3.11, matched by this job. + run: PYTHON=python bash tests/scripts/setup_stable_venv.sh origin/main + + - name: Run old-vs-new parity tests + run: pytest tests/regression/test_parity.py -v diff --git a/REFACTOR_HARDENING.md b/REFACTOR_HARDENING.md index 6673d9d..720f7ff 100644 --- a/REFACTOR_HARDENING.md +++ b/REFACTOR_HARDENING.md @@ -24,13 +24,13 @@ Tier list is the "easy wins" audit against `origin/refactor/main`. | 1b — geno honors threshold + real pass/fail | ✅ Done | `qc/steps/variant_missingness.py:79`; failure propagates as exception. | | 1c — haplotype reports failure | ✅ Fixed here | Step raises correctly; **CLI now catches it under `--warn`** (see below). | | 1d — GWAS case/control guard | ✅ Done | `association.py:175-177` uses `.get(k, 0)`. | -| 1d — GWAS log capture | ❌ Open | No log aggregation in the new path; `_all_logs.log` is header-only. | +| 1d — GWAS log capture | ✅ Addressed here | Structured logging now routed to `{out}_all_logs.log` (see below); it is no longer header-only. Raw PLINK `.log` aggregation was intentionally **not** resurrected — the FilterResult/GWASResult `log` fields still carry PLINK stderr if per-step raw capture is ever wanted. | | 1e — stop destroying input files | ✅ Done | Guarded deletes; `shutil.copy2` not `os.rename`. | | 2 — parallelize per-ancestry groups | ❌ Open | Still a serial loop (`cli/runner.py:337-354`). Pure-function design makes it cheap to add. | | 3 — cache ancestry PCA | 🟡 Partial | `--model` inference reuses fitted PCA; no `(ref_panel, common_snps)` auto-cache; `get_raw_files` still re-extracts. | -| logging (structured) | 🟡 Built, unplugged | `core/logging.py` good but `setup_logging()` never called at runtime → `logger.info` dropped. | +| logging (structured) | ✅ Wired here | `_setup_logging()` now calls `core.logging.setup_logging(log_file={out}_all_logs.log)` after `upfront_check`; step `logger.info`/`error` (with `[step]` markers) land in the consolidated log. Test: `tests/regression/test_logging.py`. | | tests | ✅ / parity unverified | 343→346 unit tests pass; golden = new-vs-new (self-consistency), not old-vs-new. | -| CI | ❌ Missing | No `.github/workflows`. | +| CI | ✅ Added here | `.github/workflows/ci.yml`: (1) unit+regression on Python 3.11 with PLINK/PLINK2 auto-downloaded; (2) parity job builds `.venv-stable` and runs `test_parity.py`. | --- @@ -77,6 +77,101 @@ bash tests/scripts/setup_stable_venv.sh origin/main # install old baseline pytest tests/regression/test_parity.py -v # old vs new ``` +### Round 2 (`refactor/hardening-round-2`) + +1. **`psutil` declared** in `setup.py::install_requires` (imported by + `ancestry/model.py` + legacy `ancestry.py`; a clean `pip install .` used to + ImportError). Verified via package metadata. +2. **`--warn` terminal-step-failure test** — `tests/unit/test_cli/test_runner_regression.py` + now covers the *last* step raising under `--warn`, asserting the last-passed + output is promoted to `{out}.pgen/.pvar/.psam` (`_handle_final_step_failure`), + plus fail-fast without `--warn`. +3. **Structured logging wired** — `cli/runner.py::_setup_logging()` now calls + `core.logging.setup_logging(level="INFO", log_file={out}_all_logs.log, + console=False)` *after* `upfront_check` (which errors if that log already + exists). It opens the file in append mode, so the ASCII banner stays on top + and structured, step-tagged (`[callrate_prune]`, …) records follow. + `setup_logging` now also closes handlers before clearing to avoid FD leaks + across repeated in-process runs. Test: `tests/regression/test_logging.py`. + + **Consolidated-log design decision:** chose the *structured-file-handler* + approach over resurrecting legacy `concat_logs` raw-PLINK-`.log` aggregation. + Rationale: single mechanism, forward-looking (step-context aware), and it + reuses the already-built `core/logging.py`. Trade-off: `{out}_cleaned_logs.log` + remains an empty placeholder (legacy filename kept for compatibility); the + populated consolidated log is `{out}_all_logs.log`. Raw PLINK output is still + available per-step in `FilterResult.log`/`GWASResult.log` if finer capture is + wanted later. + +4. **GWAS `--glm` tokenization bug fixed** — `gwas/steps/association.py`. + `run_gwas` passed `config.glm_options` (`"hide-covar firth-fallback no-x-sex + cols=..."`) as a *single* argv token. The legacy code built an f-string that + `shell_do` split on whitespace; `run_command` passes list elements verbatim, + so PLINK2 got one giant invalid `--glm` argument and produced **no GWAS + output on any run** — silently swallowed under `--warn`. (The prior PR fixed + the summary crash, but GWAS never actually ran.) Now splits the modifiers + into tokens and keeps `allow-no-covars` in the `--glm` group. Test: + `tests/unit/test_gwas/test_steps_regression.py` runs a real PLINK2 `--glm`. + +5. **Parity harness extended** — `tests/regression/test_parity.py` + + `tests/regression/compare.py`: + - per-step old/new flag map (old `--case_control`/`--full_output`, new + `--case-control`/`--full-output`; single-word `haplotype`/`ld`/`pca`/`gwas` + identical); + - multi-word QC steps `case_control`, `haplotype`, `ld` (parametrized), + `--all_sample --all_variant` full pipeline — all **parity-equal** on + synthetic data (IDs + genotype content); + - GWAS (`--pca --gwas`): new comparator `compare_gwas_results()`/`compare_gwas()` + (+ `find_gwas_output`, `_lambda_gc`) in `compare.py`, unit-tested by + `tests/regression/test_compare_gwas.py`. + All parity tests pass with `.venv-stable` present and skip cleanly without it. + +6. **CI added** — `.github/workflows/ci.yml`: a unit+regression job (Python 3.11, + PLINK/PLINK2 auto-downloaded via the dependency resolver) and a parity job + that builds `.venv-stable` from `origin/main` and runs `test_parity.py`. + Both jobs cache `~/.genotools`; the workflow also supports manual dispatch. + **The first CI run flaked, then went green:** the pre-refactor code calls + `check_king()` at *import*, so on Linux the old CLI downloaded KING from the + flaky `kingrelatedness.com` on every start (a 17-min hang across the parity + tests). Fixed by pre-caching KING in `setup_stable_venv.sh` (Linux only) plus + `~/.genotools` caching. Details: [TESTING.md §7](TESTING.md#7-king-on-linux). + +7. **Handoff tooling & docs** — `TESTING.md` (regression/parity guide for the dev + running real-cohort testing) and `tests/scripts/run_parity.py` (turnkey + old-vs-new parity on a real cohort: copies the cohort to a workdir, runs both + CLIs across QC / full-pipeline / GWAS scenarios, prints a PASS/FAIL report and + exits non-zero on divergence). + +### ✅ Resolved (decision B): PCA region exclusion is an intentional fix + +The GWAS parity run surfaced a **real old-vs-new scientific difference** in PCA +pruning (`gwas/steps/pca.py` vs legacy `gwas.py`): + +- **New** uses `--exclude range {file}` — correctly drops the high-LD/MHC ranges + in the exclusion file (synthetic data: 9741 → 9644 variants; prune.in 3971). +- **Old** used `--exclude {file}` (no `range`) — PLINK2 treats a ranges file as a + variant-ID list, matches nothing, and excludes **zero** of those regions + (9741 variants; prune.in 3987). i.e. the old region exclusion was a no-op bug. + +Consequence: the PCA eigenvectors differ (eigenvalues differ ~0.01, values +beyond sign flips), so the GWAS covariates — and hence **every** GWAS p-value — +differ slightly between old and new. The genomic-inflation **lambda still agrees** +(synthetic: 1.0074 vs 1.0071) and the **set of tested variants is identical**. + +**Decision (maintainer, option B): ratified `--exclude range` as an intentional +correctness fix** — excluding MHC/high-LD regions before PCA is standard practice, +and the old behavior was a no-op bug. This is a **deliberate, accepted divergence +from the pre-refactor GWAS baseline**, not a regression: real-cohort GWAS +per-variant p-values will differ slightly from the old code, and that is expected. + +Implications for the harness: +- The GWAS parity test asserts the invariants that hold under B (GWAS runs, same + tested-variant set, lambda within tolerance) and intentionally does **not** + assert per-variant p-equality against the old baseline. +- A guard test locks the fix in so it can't silently revert to the old no-op: + `tests/unit/test_gwas/test_steps_regression.py::TestPcaExcludesHighLdRegions` + asserts PCA pruning removes every variant inside the exclusion ranges. + --- ## Remaining work (tracked, not yet done) @@ -85,15 +180,18 @@ Priority order for making the refactor mergeable to `main`: 1. **Prove parity on real data** — run `test_parity.py` on representative real cohorts (not just synthetic), across the full QC set and, ideally, ancestry + - GWAS. This is the gate before merging to `main`. (Harness is ready.) -2. **Wire structured logging** — call `core.logging.setup_logging()` from the CLI - so step `logger.info`/error messages aren't dropped; decide on one log path. -3. **Restore GWAS/consolidated log capture** (1d) — new path aggregates no PLINK - logs; `_all_logs.log` is header-only. Regression vs legacy `concat_logs`. -4. **Add CI** — run unit + regression suites on push; add a parity job that builds - `.venv-stable` and runs `test_parity.py`. -5. **Declare `psutil`** in `setup.py` `install_requires` (imported by - `ancestry/model.py`; a clean `pip install .` currently ImportErrors). + GWAS. This is the gate before merging to `main`. (Harness extended in round 2: + multi-word QC steps, `--all_sample --all_variant`, and GWAS+lambda. Per + decision B above, real-cohort GWAS per-variant p-values will differ slightly + from the old baseline *by design* — PCA now excludes MHC/high-LD regions — so + GWAS parity is asserted at the tested-variant-set + lambda level, not + per-variant. QC and full-pipeline parity remain exact.) +2. ✅ **Wire structured logging** — DONE in round 2 (see above). +3. ✅ **Restore GWAS/consolidated log capture** (1d) — DONE in round 2 via the + structured-file-handler approach (see above). +4. ✅ **Add CI** — DONE in round 2 (`.github/workflows/ci.yml`): unit+regression + job (auto-downloads PLINK/PLINK2) + parity job (builds `.venv-stable`). +5. ✅ **Declare `psutil`** — DONE in round 2. 6. **Decouple from legacy** — new CLI still imports `utils.py` (`gt_header`/`bfiles_to_pfiles`/`vcf_to_pfiles`/`upfront_check`) and loads `ancestry.py` via an importlib file-path hack. Blocks legacy removal (Phase 5/6). diff --git a/TESTING.md b/TESTING.md new file mode 100644 index 0000000..1a5b349 --- /dev/null +++ b/TESTING.md @@ -0,0 +1,227 @@ +# GenoTools Testing & Regression Guide + +This guide explains how testing works on the `refactor/main` branch — especially +the **old-vs-new parity** testing that gates the refactor before it merges to +`main`. It's written for a developer picking this up to run regression testing on +real cohorts. + +If you only read one thing: the refactor's contract is **identical scientific +results to the pre-refactor code**. "Regression testing" here means proving that, +first on synthetic data (automated), then on real cohorts (your job — see +[Running parity on real data](#running-parity-on-real-data)). + +--- + +## 1. Environment setup + +Use the repo `.venv` (Python 3.11) for the **new** code: + +```bash +python3.11 -m venv .venv +source .venv/bin/activate +pip install -e . +pip install pytest +``` + +PLINK/PLINK2 are downloaded automatically the first time a step needs them +(into `~/.genotools/misc/executables/`). You don't install them manually. + +--- + +## 2. Running the test suite + +```bash +.venv/bin/python -m pytest tests/unit tests/regression -q +``` + +- **`tests/unit/`** — fast, hermetic unit tests of the new modular code + (`core/`, `qc/`, `gwas/`, `ancestry/`, `cli/`). No external baseline needed. +- **`tests/regression/`** — heavier tests that shell out to real PLINK and, for + parity, to the pre-refactor CLI. Tests that need something absent (PLINK, the + `.venv-stable` baseline, golden files) **skip cleanly** rather than fail. + +Everything green today: **391 passed** with `.venv-stable` present. + +--- + +## 3. Two kinds of regression testing + +### 3a. Golden regression (self-consistency of the new code) + +`tests/regression/test_qc_steps.py` compares the new code's output to committed +**golden files** under `tests/regression/golden/`. These were generated by the +new code, so they only prove the new code stays self-consistent (a step still +prunes the same samples/variants it did when the golden was captured). + +Regenerate golden files if you *intentionally* change QC behavior: + +```bash +.venv/bin/python tests/scripts/generate_golden.py \ + --geno tests/data/synthetic/genotools_test \ + --out tests/regression/golden +``` + +### 3b. Old-vs-new parity (the real regression gate) + +`tests/regression/test_parity.py` runs the **pre-refactor CLI** and the **new +CLI** on the same input and asserts they agree. This is what proves the refactor +didn't change results. It needs a `.venv-stable` containing the old code. + +**Build the baseline once:** + +```bash +PYTHON=python3.11 bash tests/scripts/setup_stable_venv.sh origin/main +``` + +This creates `.venv-stable/` with the pre-refactor `genotools` installed from +`origin/main` (the last pre-refactor commit). Python 3.11 is required — the old +`umap-learn==0.5.3` doesn't build on 3.12+. On Linux the script also pre-caches +KING (see [KING on Linux](#king-on-linux)). + +**Run parity on synthetic data:** + +```bash +.venv/bin/python -m pytest tests/regression/test_parity.py -v +``` + +What it covers (all on `tests/data/synthetic/`): + +| Test | Steps | Assertion | +|------|-------|-----------| +| `test_old_vs_new_qc_parity` | callrate, sex, het, geno, hwe | identical IDs **and** genotype content | +| `test_old_vs_new_multiword_qc_parity` | case_control / haplotype / ld | identical IDs **and** genotype content | +| `test_old_vs_new_full_pipeline_parity` | all_sample + all_variant | identical IDs **and** genotype content | +| `test_old_vs_new_gwas_parity` | pca + gwas | same tested-variant set + lambda within tolerance (see [decision B](#decision-b-gwaspca-region-exclusion)) | + +Note the old CLI uses **underscore** flags (`--case_control`, `--full_output`, +`--all_sample`) while the new CLI uses **hyphen** flags (`--case-control`, +`--full-output`, `--all-sample`); the test carries a per-step flag map so both +command lines build correctly. + +--- + +## 4. What the comparators actually check (`tests/regression/compare.py`) + +- **`compare_pfiles`** — sample and variant **ID sets** only. +- **`compare_genotypes`** — exports both pfile sets to PLINK2 `.traw` and + compares the genotype matrix order-independently. Catches same-IDs-but-changed- + genotypes and flipped allele coding, which `compare_pfiles` alone would miss. +- **`compare_gwas_results` / `compare_gwas`** — align two `--glm` association + tables on variant ID, compare per-variant p-values and the derived genomic- + inflation **lambda**. `find_gwas_output` locates the `.glm.*` file for a run. + +So QC parity is checked on real genotype content, not just which IDs survived. + +--- + +## 5. Running parity on real data + +This is the main handoff task. Use the runner, which does the same thing as the +pytest parity suite but points at your cohort: + +```bash +.venv/bin/python tests/scripts/run_parity.py \ + --geno /path/to/real_cohort \ # pfile prefix, no extension + --out /path/to/parity_workdir \ + --scenarios all # or: qc full-pipeline gwas +``` + +- `--geno` is a pfile prefix (`.pgen`/`.pvar`/`.psam`). The cohort is **copied + into `--out` before running**, so your source data dir is never written to. +- It runs old (`.venv-stable`) and new (`python -m genotools`) for each scenario, + compares, prints a PASS/FAIL summary, and exits non-zero if anything diverges. +- Scenarios: `qc`, `case_control`, `haplotype`, `ld`, `full-pipeline`, `gwas` + (default: `qc full-pipeline gwas`). + +Example output: + +``` + [PASS] qc: identical IDs + genotypes + [PASS] full-pipeline: identical IDs + genotypes + [PASS] gwas: same tested-variant set; lambda old=1.00740 new=1.00712; per-variant p-mismatches=40498 (expected under decision B ...) +``` + +For ancestry parity you need a real reference panel + labels (not in the +synthetic data); run `genotools --ancestry --ref-panel ... --ref-labels ...` +under both CLIs and compare the predicted-label assignments. The default +`genotools` entry point still uses the **legacy** ancestry code, so ancestry +parity there is near-trivial; `genotools-new` (new `AncestryModel`) is +experimental and out of scope for this gate. + +--- + +## 6. Decision B: GWAS/PCA region exclusion + +GWAS parity is **not** asserted per-variant, on purpose. The new PCA pruning uses +PLINK2 `--exclude range` and correctly drops high-LD/MHC regions; the old code +used `--exclude` on a ranges file, which PLINK2 treats as a variant-ID list — +a **no-op** that excluded none of those regions. So the PCA eigenvectors (GWAS +covariates) differ slightly, and every GWAS p-value shifts a little. The +genomic-inflation **lambda still agrees** and the **tested-variant set is +identical**. + +This was ratified as an intentional correctness fix ("decision B" in +`REFACTOR_HARDENING.md`): excluding those regions before PCA is standard practice. +So on **real data, expect GWAS per-variant p-values to differ slightly from the +pre-refactor baseline** — that is correct, not a regression. Parity is asserted +at the tested-variant-set + lambda level. A guard test +(`tests/unit/test_gwas/test_steps_regression.py::TestPcaExcludesHighLdRegions`) +locks the exclusion in so it can't silently revert to the old no-op. + +QC and full-pipeline parity **are** exact (IDs + genotype content). + +--- + +## 7. KING on Linux + +The **pre-refactor** code calls `check_king()` at *module import* (top-level in +the old `qc.py`). On **Linux**, merely starting the old CLI therefore tries to +download KING from `kingrelatedness.com` — a slow, flaky host — even for QC/GWAS +runs that never use KING. (On macOS `check_king()` returns `None`, so this never +happens.) An un-cached download once hung CI for 17 minutes of retries. + +`setup_stable_venv.sh` now pre-caches KING on Linux so the old CLI finds it and +skips the download. If that host is unreachable when you build `.venv-stable`, +the script warns you; you can drop any executable named `king` at +`~/.genotools/misc/executables/king` (parity steps never actually invoke it). +The new code uses lazy KING init and is unaffected. + +--- + +## 8. Continuous integration (`.github/workflows/ci.yml`) + +Runs on pushes to `main`/`refactor/main`, on every PR, and on manual dispatch. +Two parallel jobs on Ubuntu / Python 3.11: + +1. **Unit + regression tests** — installs the new code, downloads PLINK/PLINK2, + runs `pytest tests/unit tests/regression`. The parity tests skip here (no + `.venv-stable`). +2. **Old-vs-new parity** — builds `.venv-stable` from `origin/main` and runs + `tests/regression/test_parity.py`. + +`~/.genotools` (the downloaded PLINK/PLINK2/KING) is cached between runs, so the +external download hosts are hit at most once per cache key. + +--- + +## 9. Gotchas / FAQ + +- **`.venv-stable` won't build on Python 3.12+** — use 3.11 + (`PYTHON=python3.11 bash tests/scripts/setup_stable_venv.sh`). +- **Parity tests all skip** — you haven't built `.venv-stable`, or PLINK2 isn't + available. That's the intended fallback, not a failure. +- **Old CLI hangs on Linux** — KING download; see [KING on Linux](#king-on-linux). +- **A stray `.pheno` appears** — GWAS writes `{input}.pheno` next to the + input (pre-existing behavior, kept for parity). Both the runner and the GWAS + tests copy the input to a temp/work dir first to avoid polluting source data; + do the same in any ad-hoc runs. +- **`--warn` vs `--no-warn`** — the new CLI continues past a failed step by + default (`--warn` behavior); pass `--no-warn` to fail fast. The old CLI's + equivalent is `--warn True`/absent. +- **Where are the run logs?** — `{out}_all_logs.log` is the consolidated, + step-tagged log; `{out}.json` has the structured metrics. + +--- + +See `REFACTOR_HARDENING.md` for the full audit scorecard, what's been hardened, +and remaining work before `refactor/main` → `main`. diff --git a/genotools/cli/runner.py b/genotools/cli/runner.py index 9acfecc..7fbf7ba 100644 --- a/genotools/cli/runner.py +++ b/genotools/cli/runner.py @@ -266,7 +266,16 @@ def _initialize_modules(self) -> None: self._ancestry = Ancestry() def _setup_logging(self) -> None: - """Set up log files.""" + """Set up the consolidated run log. + + Runs *after* upfront_check (which errors if ``{out}_all_logs.log`` + already exists), so it can safely (re)create the log. Beyond creating + the legacy-named files, this attaches the structured logging file + handler so every step's ``logger.info``/``error`` is aggregated into + ``{out}_all_logs.log``. Without this call ``setup_logging`` is never + invoked at runtime and the root logger's WARNING default drops all step + logs, leaving the consolidated log header-only. + """ assert self.state is not None out_path = str(self.state.out_path) @@ -289,6 +298,13 @@ def _setup_logging(self) -> None: with open(cleaned_logs, "w") as fp: pass + # Route all genotools step/runner logs into the consolidated log file. + # setup_logging() opens the file in append mode, so the banner header + # written above stays at the top and structured records follow it. + from ..core.logging import setup_logging + + setup_logging(level="INFO", log_file=Path(all_logs), console=False) + def _convert_input_format(self) -> None: """Convert input format to pfiles if needed.""" input_format = self.args.input.input_format diff --git a/genotools/core/logging.py b/genotools/core/logging.py index 93426d1..228b592 100644 --- a/genotools/core/logging.py +++ b/genotools/core/logging.py @@ -69,7 +69,11 @@ def setup_logging( logger = logging.getLogger("genotools") logger.setLevel(level) - # Clear any existing handlers + # Clear any existing handlers, closing them first so repeated + # setup_logging() calls in one process (e.g. successive pipeline runs) + # don't leak file descriptors. + for handler in logger.handlers[:]: + handler.close() logger.handlers.clear() # Create formatter with step context placeholder diff --git a/genotools/gwas/steps/association.py b/genotools/gwas/steps/association.py index e9c299f..1690e5b 100644 --- a/genotools/gwas/steps/association.py +++ b/genotools/gwas/steps/association.py @@ -229,13 +229,22 @@ def run_gwas( pheno_path = _create_phenotype_file(data, config.pheno_name) logger.debug(f"Created phenotype file: {pheno_path}") - # Build PLINK2 command + # Build PLINK2 command. glm_options is a space-separated string of --glm + # modifiers; each must be a separate argv token. The legacy code relied + # on shell_do splitting the whole command on whitespace, but run_command + # passes list elements verbatim, so an un-split string reaches PLINK2 as + # one giant (invalid) --glm argument. allow-no-covars is itself a --glm + # modifier, so it must stay in the --glm group when no covariates apply. plink2 = get_plink2() + glm_args = config.glm_options.split() + if covar_path is None: + glm_args.append("allow-no-covars") + cmd = [ str(plink2), "--pfile" if data.format == "pfile" else "--bfile", str(data.path), - "--glm", config.glm_options, + "--glm", *glm_args, "--pheno-name", config.pheno_name, "--pheno", str(pheno_path), ] @@ -256,9 +265,6 @@ def run_gwas( if config.covariate_variance_standardize: cmd.append("--covar-variance-standardize") - else: - # No covariates - need allow-no-covars flag - cmd[cmd.index("--glm") + 1] = config.glm_options + " allow-no-covars" cmd.extend(["--out", str(out_path)]) diff --git a/setup.py b/setup.py index 68ebd83..4f5c7a8 100755 --- a/setup.py +++ b/setup.py @@ -54,6 +54,7 @@ 'pandas>=2.0.3', 'Pillow>=9.3.0', 'plotly>=5.11.0', + 'psutil>=5.9.0', 'requests>=2.28.1', 'scikit_learn>=1.3.0', 'scipy>=1.9.3', diff --git a/tests/regression/compare.py b/tests/regression/compare.py index f2962a9..1a00b0d 100644 --- a/tests/regression/compare.py +++ b/tests/regression/compare.py @@ -405,6 +405,175 @@ def _export_traw(prefix: Path, plink2_exec: str, out_prefix: Path) -> Path: return traw +# --------------------------------------------------------------------------- +# GWAS result comparison +# +# GWAS output is NOT a pfile: it is a PLINK2 --glm association table +# ({out}.PHENO1.glm.logistic.hybrid or .glm.linear). GWAS was broken until the +# last hardening PR, so old-vs-new GWAS parity had never been checked. These +# helpers align the two result tables on variant ID and compare per-variant +# p-values (and the genomic-inflation lambda derived from them) within a small +# tolerance. +# --------------------------------------------------------------------------- + + +def find_gwas_output(out_prefix: Path) -> Optional[Path]: + """Locate a PLINK2 --glm output file for a run output prefix. + + Matches both logistic (``.glm.logistic.hybrid``) and linear + (``.glm.linear``) results, regardless of phenotype name. + + Args: + out_prefix: Run output prefix (without extension). + + Returns: + Path to the glm results file, or None if not found. + """ + out_prefix = Path(out_prefix) + patterns = ["*.glm.logistic.hybrid", "*.glm.linear"] + for pattern in patterns: + matches = sorted(out_prefix.parent.glob(f"{out_prefix.name}{pattern}")) + if matches: + return matches[0] + return None + + +def _lambda_gc(pvalues) -> float: + """Compute the genomic inflation factor (lambda GC) from p-values. + + lambda_gc = median(chi2) / expected_median, with chi2 the 1-df chi-squared + quantiles of (1 - p). Matches the approach used by the production inflation + code so a comparison here is meaningful. + """ + from scipy.stats import chi2 + + p = pd.to_numeric(pd.Series(pvalues), errors="coerce").to_numpy() + p = p[~pd.isna(p)] + if len(p) == 0: + return float("nan") + chi2_stats = chi2.isf(p, 1) + expected_median = chi2.ppf(0.5, 1) + import numpy as np + + return float(np.nanmedian(chi2_stats) / expected_median) + + +def compare_gwas_results( + expected: Path, + actual: Path, + p_tolerance: float = 1e-6, + lambda_tolerance: float = 1e-3, +) -> ComparisonResult: + """Compare two PLINK2 --glm association tables, order-independently. + + Aligns on variant ID (and TEST, when present), then compares per-variant + p-values within ``p_tolerance`` and the derived lambda GC within + ``lambda_tolerance``. ``mismatched_variants`` lists variants whose p-value + differs beyond tolerance. + + Args: + expected: Path to the expected .glm.* file. + actual: Path to the actual .glm.* file. + p_tolerance: Absolute tolerance for per-variant p-value comparison. + lambda_tolerance: Absolute tolerance for lambda GC comparison. + + Returns: + ComparisonResult over the association tables. + """ + import numpy as np + + exp = pd.read_csv(expected, sep=r"\s+", dtype={"#CHROM": str}, low_memory=False) + act = pd.read_csv(actual, sep=r"\s+", dtype={"#CHROM": str}, low_memory=False) + + # Restrict to the additive test if a TEST column is present. + if "TEST" in exp.columns: + exp = exp[exp["TEST"] == "ADD"] + if "TEST" in act.columns: + act = act[act["TEST"] == "ADD"] + + exp = exp.set_index("ID") + act = act.set_index("ID") + + exp_ids, act_ids = set(exp.index), set(act.index) + missing = exp_ids - act_ids + extra = act_ids - exp_ids + common = sorted(exp_ids & act_ids) + + mismatched_variants: list[str] = [] + if common: + e = exp.loc[common] + a = act.loc[common] + e_p = pd.to_numeric(e["P"], errors="coerce").to_numpy() + a_p = pd.to_numeric(a["P"], errors="coerce").to_numpy() + both_nan = pd.isna(e_p) & pd.isna(a_p) + # Treat matching NaNs as equal; compare finite p-values within tolerance. + diff = ~both_nan & ~(np.abs(e_p - a_p) <= p_tolerance) + mismatched_variants = [common[i] for i in np.where(diff)[0]] + + lambda_exp = _lambda_gc(exp["P"]) if "P" in exp.columns else float("nan") + lambda_act = _lambda_gc(act["P"]) if "P" in act.columns else float("nan") + lambda_ok = ( + (pd.isna(lambda_exp) and pd.isna(lambda_act)) + or abs(lambda_exp - lambda_act) <= lambda_tolerance + ) + + equal = ( + len(missing) == 0 + and len(extra) == 0 + and len(mismatched_variants) == 0 + and lambda_ok + ) + + message = ( + f"variants: missing {len(missing)}, extra {len(extra)}, " + f"p-mismatches {len(mismatched_variants)}; " + f"lambda expected {lambda_exp:.5f} vs actual {lambda_act:.5f} " + f"(ok={lambda_ok})" + ) + + return ComparisonResult( + equal=equal, + sample_diff=0, + variant_diff=len(missing) + len(extra) + len(mismatched_variants), + mismatched_samples=[], + mismatched_variants=sorted(missing | extra) + mismatched_variants, + message=message, + ) + + +def compare_gwas(expected_prefix: Path, actual_prefix: Path, **kwargs) -> ComparisonResult: + """Compare GWAS output for two run prefixes. + + Locates each run's .glm.* file and compares them via + :func:`compare_gwas_results`. + + Args: + expected_prefix: Expected run output prefix (without extension). + actual_prefix: Actual run output prefix (without extension). + **kwargs: Forwarded to compare_gwas_results (tolerances). + + Returns: + ComparisonResult, non-equal with an explanatory message if either + glm file is missing. + """ + exp_glm = find_gwas_output(Path(expected_prefix)) + act_glm = find_gwas_output(Path(actual_prefix)) + + if exp_glm is None or act_glm is None: + return ComparisonResult( + equal=False, + sample_diff=0, + variant_diff=0, + mismatched_samples=[], + mismatched_variants=[], + message=( + f"GWAS output missing: expected={exp_glm}, actual={act_glm}" + ), + ) + + return compare_gwas_results(exp_glm, act_glm, **kwargs) + + def compare_genotypes( expected_prefix: Path, actual_prefix: Path, diff --git a/tests/regression/test_compare_gwas.py b/tests/regression/test_compare_gwas.py new file mode 100644 index 0000000..5f18013 --- /dev/null +++ b/tests/regression/test_compare_gwas.py @@ -0,0 +1,111 @@ +# Copyright 2023 The GenoTools Authors. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================== + +"""Tests for the GWAS association-table comparator used by the parity harness. + +GWAS output is a PLINK2 --glm table, not a pfile, so it needs its own +comparator: align on variant ID, compare per-variant p-values and the derived +genomic-inflation lambda. These tests cover that comparator directly (no +.venv-stable / PLINK required). +""" + +from pathlib import Path + +from tests.regression.compare import ( + compare_gwas, + compare_gwas_results, + find_gwas_output, +) + +_HEADER = "#CHROM\tPOS\tID\tREF\tALT\tA1\tTEST\tOBS_CT\tOR\tP\n" + + +def _write_glm(path: Path, pvalues: list[float]) -> None: + """Write a minimal PLINK2 .glm.logistic.hybrid file with ADD-test rows.""" + lines = [_HEADER.rstrip("\n")] + for i, p in enumerate(pvalues): + lines.append(f"1\t{i + 1}\trs{i}\tA\tG\tG\tADD\t500\t1.0\t{p}") + path.write_text("\n".join(lines) + "\n") + + +_PVALS = [0.5, 0.1, 0.9, 0.3, 0.7, 0.05, 0.6, 0.2, 0.8, 0.4, 0.55] + + +def test_identical_glm_files_are_equal(tmp_path: Path) -> None: + a = tmp_path / "a.PHENO1.glm.logistic.hybrid" + b = tmp_path / "b.PHENO1.glm.logistic.hybrid" + _write_glm(a, _PVALS) + _write_glm(b, _PVALS) + + result = compare_gwas_results(a, b) + + assert result.equal, result.message + assert result.variant_diff == 0 + + +def test_perturbed_pvalue_is_detected(tmp_path: Path) -> None: + a = tmp_path / "a.PHENO1.glm.logistic.hybrid" + b = tmp_path / "b.PHENO1.glm.logistic.hybrid" + _write_glm(a, _PVALS) + perturbed = list(_PVALS) + perturbed[3] = 0.31 # beyond default p_tolerance + _write_glm(b, perturbed) + + result = compare_gwas_results(a, b) + + assert not result.equal + assert "rs3" in result.mismatched_variants + + +def test_missing_variant_is_detected(tmp_path: Path) -> None: + a = tmp_path / "a.PHENO1.glm.logistic.hybrid" + b = tmp_path / "b.PHENO1.glm.logistic.hybrid" + _write_glm(a, _PVALS) + _write_glm(b, _PVALS[:-1]) # drop one variant + + result = compare_gwas_results(a, b) + + assert not result.equal + assert result.variant_diff >= 1 + + +def test_lambda_divergence_flagged_even_when_ids_match(tmp_path: Path) -> None: + """Same variant IDs but a wildly different p-distribution → lambda mismatch.""" + a = tmp_path / "a.PHENO1.glm.logistic.hybrid" + b = tmp_path / "b.PHENO1.glm.logistic.hybrid" + _write_glm(a, [0.5] * 20) # lambda ~ 0 + _write_glm(b, [1e-8] * 20) # hugely inflated + result = compare_gwas_results(a, b, p_tolerance=1.0) # ignore per-variant p + assert not result.equal, result.message + + +def test_find_gwas_output_locates_glm(tmp_path: Path) -> None: + prefix = tmp_path / "run" + glm = tmp_path / "run.PHENO1.glm.logistic.hybrid" + _write_glm(glm, _PVALS) + + found = find_gwas_output(prefix) + assert found == glm + + +def test_compare_gwas_reports_missing_output(tmp_path: Path) -> None: + a_prefix = tmp_path / "a" + b_prefix = tmp_path / "b" + _write_glm(tmp_path / "a.PHENO1.glm.logistic.hybrid", _PVALS) + # b has no glm output + + result = compare_gwas(a_prefix, b_prefix) + assert not result.equal + assert "missing" in result.message.lower() diff --git a/tests/regression/test_logging.py b/tests/regression/test_logging.py new file mode 100644 index 0000000..87950c5 --- /dev/null +++ b/tests/regression/test_logging.py @@ -0,0 +1,79 @@ +# Copyright 2023 The GenoTools Authors. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================== + +"""Consolidated run-log integration tests. + +core/logging.py::setup_logging() was never called at runtime, so the root +logger's WARNING default dropped every step's logger.info/error and the +consolidated {out}_all_logs.log came out header-only. These tests run the real +CLI and assert the consolidated log is actually populated with structured, +step-tagged records. + +They require PLINK2 and the synthetic test data; they skip cleanly otherwise. +""" + +import subprocess +import sys +from pathlib import Path + +REPO_ROOT = Path(__file__).resolve().parents[2] + + +def _run_cli(args: list[str]) -> subprocess.CompletedProcess: + return subprocess.run( + [sys.executable, "-m", "genotools", *args], + cwd=str(REPO_ROOT), + capture_output=True, + text=True, + timeout=300, + ) + + +def test_run_log_contains_step_markers(test_geno_path: Path, tmp_path: Path) -> None: + """The consolidated {out}_all_logs.log must contain the steps' structured logs.""" + out = tmp_path / "logtest" + + res = _run_cli( + [ + "--pfile", str(test_geno_path), + "--out", str(out), + "--callrate", + "--geno", + ] + ) + assert res.returncode == 0, ( + f"Pipeline failed.\nSTDOUT:\n{res.stdout}\nSTDERR:\n{res.stderr}" + ) + + all_logs = Path(f"{out}_all_logs.log") + assert all_logs.exists(), "consolidated {out}_all_logs.log was not created" + + content = all_logs.read_text() + + # The banner header alone is short; a populated log is materially longer. + header_only_len = 400 + assert len(content) > header_only_len, ( + "consolidated log looks header-only (structured step logs were dropped):\n" + f"{content!r}" + ) + + # Structured step-context markers written by core.logging.step_context. + assert "[callrate_prune]" in content, "no callrate step marker in run log" + assert "[geno_prune]" in content, "no geno step marker in run log" + + # Actual step messages emitted via logger.info must be present. + assert "filtering complete" in content.lower(), ( + "step completion messages missing from run log" + ) diff --git a/tests/regression/test_parity.py b/tests/regression/test_parity.py index b638914..21e0a7c 100644 --- a/tests/regression/test_parity.py +++ b/tests/regression/test_parity.py @@ -19,7 +19,8 @@ golden files were generated by the NEW code, so they only prove self-consistency. These tests run the PRE-REFACTOR CLI (from .venv-stable) and the NEW CLI on the same input and assert the outputs match on sample/variant IDs AND genotype -content (via compare_genotypes, not just ID sets). +content (via compare_genotypes, not just ID sets), plus GWAS association tables +and genomic-inflation lambda. They skip cleanly when .venv-stable is absent. To enable: @@ -28,74 +29,112 @@ See REFACTOR_HARDENING.md for the parity workflow. """ +import shutil import subprocess import sys from pathlib import Path +import pandas as pd import pytest from tests.regression.compare import ( compare_pfiles, compare_genotypes, + find_gwas_output, + _lambda_gc, _resolve_plink2, ) REPO_ROOT = Path(__file__).resolve().parents[2] -# QC step flags spelled identically in BOTH the old and new CLIs. The -# persist-output flag differs (old: --full_output, new: --full-output) and is -# appended per-CLI below. -COMMON_QC_FLAGS = ["--callrate", "--sex", "--het", "--geno", "--hwe"] +# Per-step flag spelling: token -> (old_flag, new_flag). Single-word steps +# (haplotype, ld, pca, gwas) are spelled the same in both CLIs; the multi-word +# ones (case_control) and the persist-output flag differ (underscore vs hyphen). +_STEP_FLAGS = { + "callrate": ("--callrate", "--callrate"), + "sex": ("--sex", "--sex"), + "het": ("--het", "--het"), + "geno": ("--geno", "--geno"), + "case_control": ("--case_control", "--case-control"), + "haplotype": ("--haplotype", "--haplotype"), + "hwe": ("--hwe", "--hwe"), + "ld": ("--ld", "--ld"), + "all_sample": ("--all_sample", "--all-sample"), + "all_variant": ("--all_variant", "--all-variant"), + "pca": ("--pca", "--pca"), + "gwas": ("--gwas", "--gwas"), +} + +# Explicit values appended after a flag (identical for both CLIs). The old CLI's +# --ld indexes args['ld'][0..2], so bare --ld would IndexError there; pass all +# three values so both command lines are unambiguous and identical. +_STEP_ARGS = {"ld": ["50", "5", "0.5"]} + +_OLD_PERSIST = "--full_output" # old CLI spelling +_NEW_PERSIST = "--full-output" # new CLI spelling def _run(cmd: list[str], cwd: str | None = None) -> subprocess.CompletedProcess: - return subprocess.run(cmd, cwd=cwd, capture_output=True, text=True) + return subprocess.run(cmd, cwd=cwd, capture_output=True, text=True, timeout=600) -def test_old_vs_new_qc_parity( - stable_venv_genotools, - test_geno_path: Path, - tmp_path: Path, -) -> None: - """Old and new CLIs must produce identical QC output (IDs + genotypes).""" - if stable_venv_genotools is None: - pytest.skip( - "No .venv-stable found. Build it with: " - "bash tests/scripts/setup_stable_venv.sh" - ) +def _old_cmd(stable_bin: Path, geno: Path, out: Path, tokens: list[str]) -> list[str]: + cmd = [str(stable_bin), "--pfile", str(geno), "--out", str(out)] + for token in tokens: + cmd.append(_STEP_FLAGS[token][0]) + cmd.extend(_STEP_ARGS.get(token, [])) + cmd.append(_OLD_PERSIST) + return cmd - plink2 = _resolve_plink2() - if plink2 is None: - pytest.skip("plink2 not available for genotype-content comparison") +def _new_cmd(geno: Path, out: Path, tokens: list[str]) -> list[str]: + cmd = [sys.executable, "-m", "genotools", "--pfile", str(geno), "--out", str(out)] + for token in tokens: + cmd.append(_STEP_FLAGS[token][1]) + cmd.extend(_STEP_ARGS.get(token, [])) + cmd.append(_NEW_PERSIST) + return cmd + + +def _run_both( + stable_bin: Path, geno: Path, tmp_path: Path, tokens: list[str] +) -> tuple[Path, Path]: + """Run old and new CLIs with the same steps; return (old_out, new_out). + + Asserts each CLI produced pfile output (fails loudly with captured output + otherwise). + """ old_out = tmp_path / "old" new_out = tmp_path / "new" - old_cmd = [ - str(stable_venv_genotools), - "--pfile", str(test_geno_path), - "--out", str(old_out), - *COMMON_QC_FLAGS, - "--full_output", # old CLI spelling - ] - new_cmd = [ - sys.executable, "-m", "genotools", - "--pfile", str(test_geno_path), - "--out", str(new_out), - *COMMON_QC_FLAGS, - "--full-output", # new CLI spelling - ] - - old_res = _run(old_cmd) - new_res = _run(new_cmd, cwd=str(REPO_ROOT)) + old_res = _run(_old_cmd(stable_bin, geno, old_out, tokens)) + new_res = _run(_new_cmd(geno, new_out, tokens), cwd=str(REPO_ROOT)) assert old_out.with_suffix(".pgen").exists(), ( - f"Old CLI produced no output.\nSTDOUT:\n{old_res.stdout}\nSTDERR:\n{old_res.stderr}" + f"Old CLI produced no output for {tokens}.\n" + f"STDOUT:\n{old_res.stdout}\nSTDERR:\n{old_res.stderr}" ) assert new_out.with_suffix(".pgen").exists(), ( - f"New CLI produced no output.\nSTDOUT:\n{new_res.stdout}\nSTDERR:\n{new_res.stderr}" + f"New CLI produced no output for {tokens}.\n" + f"STDOUT:\n{new_res.stdout}\nSTDERR:\n{new_res.stderr}" ) + return old_out, new_out + +def _require_parity_env(stable_venv_genotools) -> tuple[Path, str]: + """Skip unless both .venv-stable and plink2 are available.""" + if stable_venv_genotools is None: + pytest.skip( + "No .venv-stable found. Build it with: " + "bash tests/scripts/setup_stable_venv.sh" + ) + plink2 = _resolve_plink2() + if plink2 is None: + pytest.skip("plink2 not available for genotype-content comparison") + return stable_venv_genotools, plink2 + + +def _assert_pfile_and_genotype_parity(old_out: Path, new_out: Path, plink2: str) -> None: id_result = compare_pfiles(old_out, new_out) assert id_result.equal, f"Sample/variant IDs differ: {id_result.message}" @@ -104,3 +143,114 @@ def test_old_vs_new_qc_parity( f"Genotype content differs between old and new: {geno_result.message}\n" f"Mismatched variants (first 10): {geno_result.mismatched_variants[:10]}" ) + + +def test_old_vs_new_qc_parity( + stable_venv_genotools, + test_geno_path: Path, + tmp_path: Path, +) -> None: + """Same-spelled QC steps must produce identical output (IDs + genotypes).""" + stable_bin, plink2 = _require_parity_env(stable_venv_genotools) + old_out, new_out = _run_both( + stable_bin, test_geno_path, tmp_path, ["callrate", "sex", "het", "geno", "hwe"] + ) + _assert_pfile_and_genotype_parity(old_out, new_out, plink2) + + +@pytest.mark.parametrize("step", ["case_control", "haplotype", "ld"]) +def test_old_vs_new_multiword_qc_parity( + stable_venv_genotools, + test_geno_path: Path, + tmp_path: Path, + step: str, +) -> None: + """Multi-word variant QC steps must be parity-equal despite flag-spelling + differences (old --case_control vs new --case-control, etc.).""" + stable_bin, plink2 = _require_parity_env(stable_venv_genotools) + old_out, new_out = _run_both(stable_bin, test_geno_path, tmp_path, [step]) + _assert_pfile_and_genotype_parity(old_out, new_out, plink2) + + +def test_old_vs_new_full_pipeline_parity( + stable_venv_genotools, + test_geno_path: Path, + tmp_path: Path, +) -> None: + """--all_sample --all_variant must be parity-equal end to end.""" + stable_bin, plink2 = _require_parity_env(stable_venv_genotools) + old_out, new_out = _run_both( + stable_bin, test_geno_path, tmp_path, ["all_sample", "all_variant"] + ) + _assert_pfile_and_genotype_parity(old_out, new_out, plink2) + + +def _read_add_results(glm_path: Path) -> pd.DataFrame: + """Read a PLINK2 --glm table, restricted to the additive (ADD) test.""" + df = pd.read_csv(glm_path, sep=r"\s+", dtype={"#CHROM": str}, low_memory=False) + if "TEST" in df.columns: + df = df[df["TEST"] == "ADD"] + return df + + +def test_old_vs_new_gwas_parity( + stable_venv_genotools, + test_geno_path: Path, + tmp_path: Path, +) -> None: + """Old and new --pca --gwas must run and agree on the association output. + + Asserts the strongest invariants that hold today: + * both CLIs produce a --glm results file (the new CLI produced NONE before + the --glm-tokenization fix), + * the same set of variants is tested, + * the genomic-inflation lambda agrees within a small tolerance. + + Per-variant p-values are NOT asserted equal: the new PCA pruning excludes + high-LD/MHC regions (`--exclude range`) that the old code left in + (`--exclude`, a no-op on a ranges file), so the PCA covariates -- and hence + every GWAS p-value -- differ slightly. That divergence is a *ratified, + intentional* correctness fix (decision B in REFACTOR_HARDENING.md): excluding + those regions before PCA is standard practice, so GWAS parity is asserted at + the tested-variant-set + lambda level, not per-variant, by design. + """ + stable_bin, _ = _require_parity_env(stable_venv_genotools) + + # Copy the input into tmp: GWAS writes {input}.pheno next to the input, so + # running against the tracked synthetic data would pollute it. + geno = tmp_path / "input" + for ext in (".pgen", ".pvar", ".psam"): + shutil.copy2(test_geno_path.with_suffix(ext), geno.with_suffix(ext)) + + old_out = tmp_path / "old" + new_out = tmp_path / "new" + old_res = _run(_old_cmd(stable_bin, geno, old_out, ["pca", "gwas"])) + new_res = _run(_new_cmd(geno, new_out, ["pca", "gwas"]), cwd=str(REPO_ROOT)) + + old_glm = find_gwas_output(old_out) + new_glm = find_gwas_output(new_out) + assert old_glm is not None, ( + f"Old CLI produced no GWAS output.\n" + f"STDOUT:\n{old_res.stdout}\nSTDERR:\n{old_res.stderr}" + ) + assert new_glm is not None, ( + f"New CLI produced no GWAS output.\n" + f"STDOUT:\n{new_res.stdout}\nSTDERR:\n{new_res.stderr}" + ) + + old_add = _read_add_results(old_glm) + new_add = _read_add_results(new_glm) + + old_ids = set(old_add["ID"]) + new_ids = set(new_add["ID"]) + assert old_ids == new_ids, ( + "GWAS tested different variant sets: " + f"old-only {len(old_ids - new_ids)}, new-only {len(new_ids - old_ids)}" + ) + + lambda_old = _lambda_gc(old_add["P"]) + lambda_new = _lambda_gc(new_add["P"]) + assert abs(lambda_old - lambda_new) <= 0.05, ( + f"Genomic-inflation lambda differs beyond tolerance: " + f"old={lambda_old:.5f}, new={lambda_new:.5f}" + ) diff --git a/tests/scripts/run_parity.py b/tests/scripts/run_parity.py new file mode 100755 index 0000000..c0f729f --- /dev/null +++ b/tests/scripts/run_parity.py @@ -0,0 +1,275 @@ +#!/usr/bin/env python +# Copyright 2023 The GenoTools Authors. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================== + +"""Run old-vs-new GenoTools parity on a real cohort. + +The pytest parity suite (tests/regression/test_parity.py) runs on the small +synthetic dataset. This script is the same idea aimed at a *real* cohort: it runs +the pre-refactor CLI (from .venv-stable) and the new CLI on the same input and +compares their outputs. + +Prerequisites +------------- +1. Build the pre-refactor baseline once: + PYTHON=python3.11 bash tests/scripts/setup_stable_venv.sh origin/main +2. Run this script with the new code installed in the current environment + (``pip install -e .``). + +Usage +----- + python tests/scripts/run_parity.py --geno /path/to/cohort --out /path/to/workdir + python tests/scripts/run_parity.py --geno cohort --out work --scenarios qc gwas + +``--geno`` is a pfile prefix (without .pgen/.pvar/.psam). The cohort is copied +into ``--out`` before running, so the source data dir is never written to. + +Scenarios +--------- + qc callrate, sex, het, geno, hwe (IDs + genotype content) + case_control case_control (IDs + genotype content) + haplotype haplotype (IDs + genotype content) + ld ld 50 5 0.5 (IDs + genotype content) + full-pipeline all_sample + all_variant (IDs + genotype content) + gwas pca + gwas (tested-variant set + lambda) + +GWAS parity is asserted at the tested-variant-set + genomic-inflation-lambda +level, NOT per-variant p-value: the new PCA pruning intentionally excludes +high-LD/MHC regions the old code left in (ratified "decision B" in +REFACTOR_HARDENING.md), so every GWAS p-value differs slightly by design. The +per-variant mismatch count is still reported, for information. + +Exit status is non-zero if any selected scenario fails its parity check. +""" + +from __future__ import annotations + +import argparse +import platform +import shutil +import subprocess +import sys +from pathlib import Path + +REPO_ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(REPO_ROOT)) + +from tests.regression.compare import ( # noqa: E402 + compare_pfiles, + compare_genotypes, + compare_gwas_results, + find_gwas_output, + _lambda_gc, + _resolve_plink2, +) + +# token -> (old_flag, new_flag). Single-word steps are spelled the same. +_STEP_FLAGS = { + "callrate": ("--callrate", "--callrate"), + "sex": ("--sex", "--sex"), + "het": ("--het", "--het"), + "geno": ("--geno", "--geno"), + "case_control": ("--case_control", "--case-control"), + "haplotype": ("--haplotype", "--haplotype"), + "hwe": ("--hwe", "--hwe"), + "ld": ("--ld", "--ld"), + "all_sample": ("--all_sample", "--all-sample"), + "all_variant": ("--all_variant", "--all-variant"), + "pca": ("--pca", "--pca"), + "gwas": ("--gwas", "--gwas"), +} +_STEP_ARGS = {"ld": ["50", "5", "0.5"]} # old CLI requires all three ld values + +SCENARIOS = { + "qc": ["callrate", "sex", "het", "geno", "hwe"], + "case_control": ["case_control"], + "haplotype": ["haplotype"], + "ld": ["ld"], + "full-pipeline": ["all_sample", "all_variant"], + "gwas": ["pca", "gwas"], +} +GWAS_SCENARIOS = {"gwas"} +DEFAULT_SCENARIOS = ["qc", "full-pipeline", "gwas"] + + +def _build_cmd(base: list[str], geno: Path, out: Path, tokens: list[str], new: bool) -> list[str]: + cmd = list(base) + ["--pfile", str(geno), "--out", str(out)] + for token in tokens: + cmd.append(_STEP_FLAGS[token][1 if new else 0]) + cmd.extend(_STEP_ARGS.get(token, [])) + cmd.append("--full-output" if new else "--full_output") + return cmd + + +def _run(cmd: list[str], cwd: str | None = None, timeout: int = 3600) -> subprocess.CompletedProcess: + return subprocess.run(cmd, cwd=cwd, capture_output=True, text=True, timeout=timeout) + + +def _copy_pfileset(src_prefix: Path, dst_prefix: Path) -> None: + dst_prefix.parent.mkdir(parents=True, exist_ok=True) + for ext in (".pgen", ".pvar", ".psam"): + shutil.copy2(src_prefix.with_suffix(ext), dst_prefix.with_suffix(ext)) + + +def _check_king_hint(stable_bin: Path) -> None: + """Warn if the old baseline will hang fetching KING at import on Linux.""" + if platform.system() != "Linux": + return + king = Path.home() / ".genotools" / "misc" / "executables" / "king" + if not king.exists(): + print( + "WARNING: on Linux the pre-refactor CLI downloads KING at import " + "(kingrelatedness.com) and may hang. Re-run setup_stable_venv.sh " + "(it pre-caches KING) or place a 'king' binary at " + f"{king}.\n", + file=sys.stderr, + ) + + +def _parity_qc(old_out: Path, new_out: Path, plink2: str) -> tuple[bool, str]: + if not old_out.with_suffix(".pgen").exists(): + return False, "old CLI produced no output" + if not new_out.with_suffix(".pgen").exists(): + return False, "new CLI produced no output" + ids = compare_pfiles(old_out, new_out) + if not ids.equal: + return False, f"IDs differ: {ids.message}" + geno = compare_genotypes(old_out, new_out, plink2_exec=plink2) + if not geno.equal: + return False, f"genotypes differ: {geno.message}" + return True, "identical IDs + genotypes" + + +def _parity_gwas(old_out: Path, new_out: Path) -> tuple[bool, str]: + old_glm = find_gwas_output(old_out) + new_glm = find_gwas_output(new_out) + if old_glm is None or new_glm is None: + return False, f"missing GWAS output (old={old_glm}, new={new_glm})" + + # Full comparison for the informational per-variant mismatch count. + detail = compare_gwas_results(old_glm, new_glm) + + import pandas as pd + + def _ids_and_lambda(path: Path) -> tuple[set, float]: + df = pd.read_csv(path, sep=r"\s+", dtype={"#CHROM": str}, low_memory=False) + if "TEST" in df.columns: + df = df[df["TEST"] == "ADD"] + return set(df["ID"]), _lambda_gc(df["P"]) + + old_ids, lam_old = _ids_and_lambda(old_glm) + new_ids, lam_new = _ids_and_lambda(new_glm) + + if old_ids != new_ids: + return False, ( + f"different tested-variant sets (old-only {len(old_ids - new_ids)}, " + f"new-only {len(new_ids - old_ids)})" + ) + if abs(lam_old - lam_new) > 0.05: + return False, f"lambda differs: old={lam_old:.5f}, new={lam_new:.5f}" + + return True, ( + f"same tested-variant set; lambda old={lam_old:.5f} new={lam_new:.5f}; " + f"per-variant p-mismatches={len(detail.mismatched_variants)} " + f"(expected under decision B: PCA excludes MHC/high-LD)" + ) + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("--geno", required=True, type=Path, help="Real cohort pfile prefix (no extension)") + parser.add_argument("--out", required=True, type=Path, help="Working directory for the two runs") + parser.add_argument( + "--stable-bin", + type=Path, + default=REPO_ROOT / ".venv-stable" / "bin" / "genotools", + help="Path to the pre-refactor genotools (default: .venv-stable/bin/genotools)", + ) + parser.add_argument( + "--scenarios", + nargs="+", + choices=[*SCENARIOS.keys(), "all"], + default=DEFAULT_SCENARIOS, + help=f"Scenarios to run (default: {' '.join(DEFAULT_SCENARIOS)})", + ) + parser.add_argument("--timeout", type=int, default=3600, help="Per-CLI timeout in seconds") + args = parser.parse_args() + + scenarios = list(SCENARIOS) if "all" in args.scenarios else args.scenarios + + if not args.geno.with_suffix(".pgen").exists(): + parser.error(f"input not found: {args.geno}.pgen") + if not args.stable_bin.exists(): + parser.error( + f"pre-refactor CLI not found at {args.stable_bin}. Build it with: " + "PYTHON=python3.11 bash tests/scripts/setup_stable_venv.sh origin/main" + ) + plink2 = _resolve_plink2() + if plink2 is None: + parser.error("plink2 not found (needed for genotype comparison)") + + _check_king_hint(args.stable_bin) + + args.out.mkdir(parents=True, exist_ok=True) + # Copy the cohort so GWAS's {input}.pheno side-effect never touches source data. + geno = args.out / "input" + print(f"Copying cohort {args.geno} -> {geno}.*") + _copy_pfileset(args.geno, geno) + + old_base = [str(args.stable_bin)] + new_base = [sys.executable, "-m", "genotools"] + + results: list[tuple[str, bool, str]] = [] + for name in scenarios: + tokens = SCENARIOS[name] + scen_dir = args.out / name + scen_dir.mkdir(parents=True, exist_ok=True) + old_out = scen_dir / "old" + new_out = scen_dir / "new" + + print(f"\n=== scenario: {name} ({' '.join(tokens)}) ===") + old_res = _run(_build_cmd(old_base, geno, old_out, tokens, new=False), timeout=args.timeout) + new_res = _run(_build_cmd(new_base, geno, new_out, tokens, new=True), cwd=str(REPO_ROOT), timeout=args.timeout) + + if name in GWAS_SCENARIOS: + ok, msg = _parity_gwas(old_out, new_out) + else: + ok, msg = _parity_qc(old_out, new_out, plink2) + + if not ok and (old_res.returncode != 0 or new_res.returncode != 0): + msg += ( + f"\n old rc={old_res.returncode} new rc={new_res.returncode}" + f"\n old stderr tail: {old_res.stderr.strip()[-300:]}" + f"\n new stderr tail: {new_res.stderr.strip()[-300:]}" + ) + results.append((name, ok, msg)) + print(f" {'PASS' if ok else 'FAIL'}: {msg}") + + print("\n" + "=" * 70) + print("PARITY SUMMARY") + print("=" * 70) + for name, ok, msg in results: + print(f" [{'PASS' if ok else 'FAIL'}] {name}: {msg.splitlines()[0]}") + + failed = [n for n, ok, _ in results if not ok] + if failed: + print(f"\n{len(failed)} scenario(s) FAILED parity: {', '.join(failed)}") + return 1 + print(f"\nAll {len(results)} scenario(s) passed parity.") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/scripts/setup_stable_venv.sh b/tests/scripts/setup_stable_venv.sh index 6900e04..8888a91 100755 --- a/tests/scripts/setup_stable_venv.sh +++ b/tests/scripts/setup_stable_venv.sh @@ -16,6 +16,10 @@ set -euo pipefail REF="${1:-origin/main}" +# Interpreter used to create .venv-stable. Override to pin a version that the +# pre-refactor deps (e.g. umap-learn==0.5.3 / numba) still support, e.g. +# PYTHON=python3.11 bash tests/scripts/setup_stable_venv.sh +PYTHON="${PYTHON:-python3}" ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/../.." && pwd)" WORKTREE="$(mktemp -d)/genotools-stable-src" @@ -34,7 +38,7 @@ fi git -C "${ROOT}" worktree add --detach "${WORKTREE}" "${REF}" trap 'git -C "${ROOT}" worktree remove --force "${WORKTREE}" 2>/dev/null || true' EXIT -python3 -m venv "${ROOT}/.venv-stable" +"${PYTHON}" -m venv "${ROOT}/.venv-stable" "${ROOT}/.venv-stable/bin/pip" install --quiet --upgrade pip "${ROOT}/.venv-stable/bin/pip" install "${WORKTREE}" # Runtime deps the old code needs that its setup.py omits / that fresh venvs on @@ -42,6 +46,33 @@ python3 -m venv "${ROOT}/.venv-stable" # by umap-learn==0.5.3) and psutil (imported by the old ancestry.py). "${ROOT}/.venv-stable/bin/pip" install --quiet "setuptools<81" psutil +# Pre-cache KING (Linux only). The pre-refactor genotools runs check_king() at +# *module import* (top-level in qc.py), so on Linux merely starting the old CLI +# downloads KING from kingrelatedness.com -- a slow, flaky host -- even for QC/ +# GWAS steps that never use KING. That download hung CI for 17 min of retries. +# Fetch the binary once here so the old CLI finds it and skips the download. +# (macOS: check_king() returns None, so there is nothing to do.) +if [ "$(uname -s)" = "Linux" ]; then + KING_DIR="${GENOTOOLS_DEP_DIR:-${HOME}/.genotools/misc/executables}" + if [ -x "${KING_DIR}/king" ]; then + echo ">> KING already cached at ${KING_DIR}/king" + else + echo ">> Pre-caching KING for the old baseline (Linux)..." + mkdir -p "${KING_DIR}" + if curl -fsSL --retry 3 --retry-delay 5 --max-time 180 \ + -o /tmp/king.tar.gz \ + https://www.kingrelatedness.com/executables/Linux-king232.tar.gz; then + tar -xzf /tmp/king.tar.gz -C "${KING_DIR}" && chmod +x "${KING_DIR}/king" + echo ">> KING cached at ${KING_DIR}/king" + else + echo "!! WARNING: could not download KING (kingrelatedness.com unreachable)." + echo "!! The old baseline calls check_king() at import, so on Linux the old" + echo "!! CLI may hang trying to fetch it. Retry, or drop a 'king' binary at" + echo "!! ${KING_DIR}/king (any executable works -- parity steps never call it)." + fi + fi +fi + echo ">> Done. Old baseline installed at:" echo " ${ROOT}/.venv-stable/bin/genotools" echo ">> Run parity tests with:" diff --git a/tests/unit/test_cli/test_runner_regression.py b/tests/unit/test_cli/test_runner_regression.py index f291818..6e3aa67 100644 --- a/tests/unit/test_cli/test_runner_regression.py +++ b/tests/unit/test_cli/test_runner_regression.py @@ -63,6 +63,29 @@ def fake_single_step( return fake_single_step +def _write_traceable_pfiles(prefix: Path) -> None: + """Create pfiles whose contents encode their own prefix. + + This lets a test tell *which* step's output ended up promoted to the final + output path (the bytes are copied verbatim by shutil.copy2). + """ + for ext in (".pgen", ".pvar", ".psam"): + Path(f"{prefix}{ext}").write_text(f"{prefix.name}{ext}") + + +def _traceable_step_factory(failing_step: str): + """Like _fake_step_factory but passing steps write traceable pfiles.""" + def fake_single_step( + step: str, step_input: str, step_output: str, legacy_args: Dict[str, Any] + ) -> Dict[str, Any]: + if step == failing_step: + raise QCError(f"{step} failed") + _write_traceable_pfiles(Path(step_output)) + return {"pass": True, "step": step, "metrics": {}, "output": {}} + + return fake_single_step + + class TestWarnModeContinuesAfterRaisingStep: """Regression: refactored steps RAISE on failure, but _run_qc_pipeline had no try/except, so under --warn a raising step aborted the entire run instead of @@ -99,3 +122,59 @@ def test_no_warn_still_fails_fast( geno_path=str(geno), out_path=str(out), ) + + +class TestWarnModeFinalStepFailurePromotesLastPassedOutput: + """Regression: when the LAST step raises under --warn, the pipeline must still + produce final output pfiles at {out}.pgen/.pvar/.psam by promoting the output + of the last step that passed (handled by _handle_final_step_failure). Only a + *middle* raising step was previously covered; a terminal failure left no final + output at all. + """ + + def test_terminal_failure_promotes_last_passed_output( + self, tmp_path: Path, monkeypatch: pytest.MonkeyPatch + ) -> None: + runner, geno, out = _make_runner(tmp_path, warn_only=True) + # het is the LAST step and raises; callrate + sex pass. + monkeypatch.setattr( + runner, "_run_single_step", _traceable_step_factory("het") + ) + + out_dict = runner._run_qc_pipeline( + steps=["callrate", "sex", "het"], + geno_path=str(geno), + out_path=str(out), + ) + + pf = out_dict["pass_fail"] + assert pf["callrate"]["status"] is True + assert pf["sex"]["status"] is True + assert pf["het"]["status"] is False # terminal failure recorded, not raised + + # Final output must exist despite the terminal failure... + for ext in (".pgen", ".pvar", ".psam"): + final = Path(f"{out}{ext}") + assert final.exists(), f"final output {final.name} was not produced" + + # ...and it must be the LAST PASSED step's output (sex), not the raw input + # or the callrate output. sex ran on the callrate output, so its prefix is + # "{out}_callrate_sex". + expected_prefix = f"{out}_callrate_sex" + assert Path(f"{out}.pgen").read_text() == f"{Path(expected_prefix).name}.pgen" + + def test_terminal_failure_without_warn_raises( + self, tmp_path: Path, monkeypatch: pytest.MonkeyPatch + ) -> None: + """Sanity check: without --warn a terminal failure still fails fast.""" + runner, geno, out = _make_runner(tmp_path, warn_only=False) + monkeypatch.setattr( + runner, "_run_single_step", _traceable_step_factory("het") + ) + + with pytest.raises(QCError): + runner._run_qc_pipeline( + steps=["callrate", "sex", "het"], + geno_path=str(geno), + out_path=str(out), + ) diff --git a/tests/unit/test_gwas/test_steps_regression.py b/tests/unit/test_gwas/test_steps_regression.py index 4243f17..4109751 100644 --- a/tests/unit/test_gwas/test_steps_regression.py +++ b/tests/unit/test_gwas/test_steps_regression.py @@ -21,7 +21,19 @@ from genotools.core.genotypes import GenotypeData from genotools.gwas.config import GWASConfig -from genotools.gwas.steps.association import _process_logistic_results +from genotools.gwas.steps.association import _process_logistic_results, run_gwas + +SYNTHETIC = ( + Path(__file__).resolve().parents[2] / "data" / "synthetic" / "genotools_test" +) + + +def _plink2_available() -> bool: + import shutil + + if shutil.which("plink2"): + return True + return (Path.home() / ".genotools" / "misc" / "executables" / "plink2").exists() def _write_binary_cohort(tmp_path: Path, n: int = 10) -> GenotypeData: @@ -73,3 +85,105 @@ def test_binary_gwas_summary_does_not_crash(self, tmp_path: Path) -> None: assert result.phenotype_type == "binary" assert result.inflation is not None assert result.inflation.lambda_1000 is not None # cases+controls present + + +class TestGwasGlmOptionsAreTokenized: + """Regression: run_gwas passed the multi-modifier --glm options as a SINGLE + argv token. shell_do (old code) split the command on whitespace; run_command + (new code) treats each list element as one argument, so PLINK2 received + 'hide-covar firth-fallback no-x-sex cols=...' as one token and rejected it + ('Invalid --glm argument ...'), producing no GWAS output. This broke every + GWAS run in the new CLI, silently under --warn. + """ + + @pytest.mark.skipif(not _plink2_available(), reason="plink2 not available") + def test_run_gwas_produces_output(self, tmp_path: Path) -> None: + """A real PLINK2 --glm run must succeed and write a results file.""" + if not SYNTHETIC.with_suffix(".pgen").exists(): + pytest.skip("synthetic test data not found") + + # Copy the pfileset into tmp: _create_phenotype_file writes {input}.pheno + # next to the input, so running on the tracked data would pollute it. + import shutil + + local = tmp_path / "cohort" + for ext in (".pgen", ".pvar", ".psam"): + shutil.copy2(SYNTHETIC.with_suffix(ext), local.with_suffix(ext)) + + data = GenotypeData.from_path(local) + out = tmp_path / "gwas_out" + + # No covariates -> exercises the allow-no-covars branch of the command. + result = run_gwas(data, GWASConfig(), out) + + assert result.success is True + assert result.output_path is not None and result.output_path.exists() + assert result.n_variants_tested > 0 + + +class TestPcaExcludesHighLdRegions: + """Ratified behavior (REFACTOR_HARDENING.md, decision B): PCA pruning excludes + the high-LD/MHC regions via PLINK2 `--exclude range`. The pre-refactor code + used `--exclude `, which PLINK2 treats as a variant-ID list -- a + no-op that excluded none of those regions. Excluding them before PCA is + standard practice, so this behavior was intentionally kept; this test guards + against silently regressing to the old no-op. + """ + + def _parse_ranges(self) -> list[tuple[str, int, int]]: + from genotools.gwas.config import get_exclusion_regions + + ranges = [] + for line in get_exclusion_regions("hg38").strip().splitlines(): + parts = line.split() + ranges.append((str(parts[0]), int(parts[1]), int(parts[2]))) + return ranges + + def _count_in_ranges(self, pvar_path: Path, ranges) -> int: + import pandas as pd + + df = pd.read_csv( + pvar_path, + sep="\t", + comment="#", + header=None, + usecols=[0, 1], + names=["chrom", "pos"], + dtype={"chrom": str}, + ) + total = 0 + for chrom, start, end in ranges: + total += df[ + (df["chrom"] == chrom) + & (df["pos"] >= start) + & (df["pos"] <= end) + ].shape[0] + return total + + @pytest.mark.skipif(not _plink2_available(), reason="plink2 not available") + def test_pca_pruning_drops_excluded_regions(self, tmp_path: Path) -> None: + if not SYNTHETIC.with_suffix(".pgen").exists(): + pytest.skip("synthetic test data not found") + + import shutil + + from genotools.gwas.config import PCAPruningConfig + from genotools.gwas.steps.pca import run_pca_pruning + + local = tmp_path / "cohort" + for ext in (".pgen", ".pvar", ".psam"): + shutil.copy2(SYNTHETIC.with_suffix(ext), local.with_suffix(ext)) + + ranges = self._parse_ranges() + + # Non-vacuous: the input actually contains variants inside the excluded + # ranges (also confirms chromosome naming matches between pvar and ranges). + assert self._count_in_ranges(local.with_suffix(".pvar"), ranges) > 0 + + data = GenotypeData.from_path(local) + out = tmp_path / "pruned" + result = run_pca_pruning(data, PCAPruningConfig(), out, build="hg38") + assert result.success is True + + # Ratified: no excluded-region variant survives pruning. + assert self._count_in_ranges(out.with_suffix(".pvar"), ranges) == 0