# GWAS Harmonizer controlled perturbation benchmark

## Status

Prepared on 16 August 2026 and completed with the final conservative-rsID rerun on 17 August 2026. Five full-file evidence packages were scored against deterministic row-level truth. The initial parallel submission encountered a client-side preview-response timeout; sequential reruns completed all scenarios without a processing failure.

## Rationale

The benchmark tests whether GWAS Harmonizer recovers known allele orientations and statistical values after controlled, biologically equivalent perturbations. This is more informative than comparing only aggregate row counts because the correct result is known for every input row. A clean control estimates false correction, systematic scenarios isolate individual allele operations, and a larger mixed scenario tests heterogeneous errors in one file.

## Reference data and sampling

The source was the official GWAS Catalog harmonized release for accession `GCST90018642` (hypothyroidism; GRCh38). The local source file had 12,181,200 data rows, SHA-256 `7f0e747b719a54dbfc81ba07ecbe7568962687b7ff8be06def06f79f266bc2e4`, and MD5 `73e4e5dc395c3c2135697f87bda13fa5`. The MD5 matched the official harmonized release.

Eligible records were autosomal or sex-chromosome, non-palindromic, biallelic A/C/G/T SNPs with a positive coordinate, valid rsID, finite beta and standard error, positive standard error, effect-allele frequency strictly between zero and one, and p-value between zero and one. Of 12,181,200 records scanned, 9,834,126 were eligible. The 100,000 records with the lowest deterministic SHA-256 priorities under seed `20260816` were selected, avoiding dependence on source row order or pseudo-random library state.

## Perturbations

Four paired 20,000-row files were generated from the same variants:

1. clean control: no perturbation;
2. effect/other allele swap, with beta multiplied by −1 and EAF transformed to `1 − EAF`;
3. strand complement of both alleles without a statistical sign change;
4. strand complement plus effect/other swap, with beta multiplied by −1 and EAF transformed to `1 − EAF`.

A separate 100,000-row mixed file contained 25,000 unchanged rows; 20,000 effect/other swaps; 20,000 strand complements; 20,000 strand complements plus swaps; 5,000 missing rsIDs; 5,000 deliberately conflicting rsIDs; and 5,000 deliberately inconsistent reported p-values. Class assignment was deterministic and dispersed throughout the file.

## Harmonization configuration

The runs used strict submission mode, target build GRCh38, reference alignment, strand correction, missing-rsID assignment, existing-rsID validation, quarantine of unresolved palindromic SNPs and conflicting duplicates, and failure on p-value inconsistency. Study semantics were explicitly confirmed as binary trait, log-odds effect scale, and standard error on the same log-odds scale. The source build was detected as GRCh38 with high confidence during the initial previews.

## Outcomes

The prespecified primary outcomes are exact scientific-field recovery, retained-row accuracy, unresolved or excluded rate, false-correction rate in unchanged rows, audited allele-action accuracy, rsID recovery and conflict handling, p-value inconsistency sensitivity and false-positive rate, and audit coverage. Binomial rates are accompanied by two-sided 95% Wilson confidence intervals. Source rows whose reported p-value was exactly zero remain part of allele-transformation scoring but are excluded from the p-value false-positive denominator because a Wald-consistency comparison is not numerically evaluable.

## Results

Across 180,000 scenario-row evaluations, all coordinates, alleles, beta values, standard errors, effect-allele frequencies, and reported p-values were recovered exactly after harmonization. Recovery, retained-row accuracy, allele-action accuracy, and auditable-disposition coverage were 100% in the clean control, each systematic perturbation, and every mixed-file class. No benchmark row was unresolved or excluded. The lower bound of the two-sided 95% Wilson interval was 99.981% for each 20,000-row systematic scenario, 99.985% for the 25,000 unchanged mixed rows, and 99.923% for each 5,000-row class.

All 5,000 deliberately inconsistent p-values were detected, while none of the 95,000 non-injected mixed-file rows were flagged. Thus p-value inconsistency sensitivity was 100% (95% CI 99.923–100.000%) and the observed false-positive rate was 0% (upper 95% Wilson bound 0.004%). The strict policy blocked submission readiness but retained the rows and their reported p-values in the harmonized table, so the result should be described as detection and blocking rather than row removal.

Identifier handling followed the prespecified conservative principle that false resolution is worse than non-resolution. Of 5,000 deliberately blank identifiers, 4,002 were assigned correctly and 998 remained unresolved. One assigned value differed textually from the historical source rsID but was verified by NCBI dbSNP as its current merged canonical identifier, so it was counted as a correct canonical-equivalent assignment rather than a false correction. Assignment precision was 100.00% (4,002/4,002) and assignment recall was 80.04% (4,002/5,000). All 5,000 deliberately conflicting rsIDs were cleared, giving 100.00% conflict-detection sensitivity and zero wrong identifiers retained in this class. The false-correction rate was 0% across the 90,000 mixed-file rows without an rsID perturbation. This precision-first policy deliberately lowers identifier retention when an rsID cannot be verified in the allele-aware lookup; unresolved identifiers are reported as missing rather than carried forward.

| Identifier endpoint | Result |
|---|---:|
| Assignment precision | 100.00% (4,002/4,002) |
| Assignment recall | 80.04% (4,002/5,000) |
| Conflict-detection sensitivity | 100.00% (5,000/5,000) |
| False-correction rate | 0.00% (0/90,000) |

## Execution record and authoritative submission verdicts

| Scenario | Run identifier | Final state | Submission ready | Blocking issues |
|---|---|---|---|---|
| Clean control | `040289d5-23a8-4a41-8075-c09012ad6c7f` | Ready | No | `rsid_validation_not_allele_aware`; `gwas_catalog_schema_not_requested` |
| Systematic effect swap | `a4bb5728-9061-4cec-8215-96e8584dad67` | Ready | No | `ambiguous_confidence_interval`; `rsid_validation_not_allele_aware`; `gwas_catalog_schema_not_requested` |
| Systematic strand complement | `643598e9-8772-4f31-8f7b-f900eb822108` | Ready | No | `rsid_validation_not_allele_aware`; `gwas_catalog_schema_not_requested` |
| Systematic strand and swap | `dfe0a555-e634-4a46-8998-7ff7de2fffa2` | Ready | No | `ambiguous_confidence_interval`; `rsid_validation_not_allele_aware`; `gwas_catalog_schema_not_requested` |
| Sporadic mixed perturbations | `42ddd843-3222-49cb-aab0-f4a5bec0384c` | Ready | No | `severe_pvalue_inconsistency`; `position_conflicting_rsids`; `allele_conflicting_rsids`; `ambiguous_confidence_interval`; `gwas_catalog_schema_not_requested` |

The severe p-value and rsID-conflict blockers were deliberately induced and therefore show that the submission gate detected the injected problems. `gwas_catalog_schema_not_requested` reflects the benchmark configuration rather than a harmonization error. By contrast, `ambiguous_confidence_interval` was raised even though the benchmark inputs contained no confidence-interval columns; this is a false-positive blocker that should be corrected before using submission-readiness rates as a manuscript endpoint. The four original systematic evidence bundles predated the allele-aware validation deployment and retain their recorded `rsid_validation_not_allele_aware` limitation; identifier performance is therefore reported only from the final mixed rerun.

## Reproducibility artifacts

- Preparation and scoring implementation: `scripts/controlled_perturbation_benchmark.py`
- Automated tests: `app/tests/test_controlled_perturbation_benchmark.py`
- Local protocol and manifest: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/`
- Local row-level truth: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/ground_truth.tsv.gz`
- Local benchmark inputs: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/inputs/`
- Scored paper-results table: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/results/PAPER_RESULTS.md`
- Machine-readable results: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/results/controlled_perturbation_results.json`
- Long-form metrics: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/results/controlled_perturbation_metrics.csv`
- Source-traceable rsID equivalence evidence: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/rsid_equivalences.tsv`
- Per-scenario packages and paper evidence: `/Users/mlemsalu/Desktop/gwas latest/GCST90018642_controlled_perturbation_benchmark_20260816/runs/`

The large generated inputs and truth file are intentionally excluded from Git.

## Source documentation

- GWAS Catalog study record: <https://www.ebi.ac.uk/gwas/rest/api/studies/GCST90018642>
- Official harmonized-release metadata: <https://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST90018001-GCST90019000/GCST90018642/harmonised/GCST90018642.h.tsv.gz-meta.yaml>
- Source publication: Sakaue et al., *Nature Genetics* (2021), <https://www.nature.com/articles/s41588-021-00931-x>

## Interpretation boundary

This benchmark evaluates harmonization quality for selected non-palindromic biallelic SNP transformations. It does not establish cohort validity, association-model validity, biological significance, cross-study comparability, indel correctness, palindromic resolution, or GRCh37-to-GRCh38 liftover accuracy.
