Commit Graph
28 Commits
Author SHA1 Message Date
kemal 88f74c2fe7 cloud: optional AZ for spot launches 2026-09-28 17:44:26 +01:00
kemal 8cc1e52f69 P8: GSE165845 validation scripts (prepare, §16.5 metrics) and cloud task
scripts/p8/prepare.py lays out the 360 GEO IDATs under their header barcode/position;
scripts/p8/validate_p8.py computes concordance, NRD, no-call rate, per-class and per-MAF
sensitivity/specificity/precision/F1 vs NYGC 30x and GIAB HG001, replicate concordance,
systematic/allele-swap site errors, per-chip batch table and duplicate detection.
Cloud userdata can skip the Omni2.5 archive (OMNI=0).
2026-09-28 17:43:31 +01:00
kemal b699d2c188 P7: differential-testing harness (simurg compare) and nightly run on the Pi
- compare/core: EXACT/NEAR/DISCREPANCY/UNKNOWN with the plan §16.3 tolerances (abs, float32 ULP), NaN = NaN,
  locus classes (assay, variant type, chromosome group, intensity-only); fields absent on both sides are
  not applicable
- compare/records: GTC vs GTC (per-locus arrays, transforms per bin, controls, sample statistics; metadata
  listed, not counted) and VCF vs VCF (records keyed incl. ID; SNV/MNV/indel, merged, multi-allelic classes)
- simurg compare A B (files or folders paired by name): summary.json, discrepancies.tsv.gz (never truncated),
  self-contained dashboard.html; --tolerance, --fail-on-discrepancy (exit 3)
- scripts/nightly.py + install_nightly.sh (crontab 02:30): unit tests, GSA vs gtc2vcf (GTC + VCF), VCF
  regression vs golden, Omni2.5 subset vs vendor GTCs; first run PASS on every step
- tests/unit/test_compare.py; docs/validation/P7.md, README
2026-09-28 17:00:37 +01:00
kemal 4d190d559e P6 complete: optional VCF fields, SNV+indel merging, auxiliary MNV loci, PLINK 2 .pgen
- gtc-to-vcf --info-fields ASSAY_TYPE,NORM_BIN,GENTRAIN,CLUSTER_SEP,ILMN_STRAND,PROBE_IDS (--egt) and
  --extra-format X,Y,NORMX,NORMY,THETA,R, in per-sample VCFs and the --merge BCF (OUT-012, plan §8.3);
  equal to gtc2vcf's values (THETA/R missing where both normalized intensities are 0, K-023)
- --merge-snv-indel (OUT-021, EDGE-013): alleles extended to the longest REF; SNV probes read the base at POS
- --auxiliary-loci (OUT-035, EDGE-015): rules (a)-(e), E_AUX_INVALID; rule (d) uses the last plus-strand base
  for minus-strand designs (spec amended, as GTCtoVCF)
- export plink --format bed|pgen|both: PLINK 2 fixed-width .pgen/.pvar/.psam; plink2 reads it back to a .bed
  byte-identical to ours
- cohort store keeps raw/normalized intensities for the BCF path; tests (unit + gtc2vcf oracle); docs
2026-09-28 14:29:42 +01:00
kemal 5b85e4e35f P6 validated against GTCtoVCF too: zero unexplained genotype differences
- GTCtoVCF (Apache-2.0, 36 GSA samples, merged sites): 23.66 M identical calls; every difference named
  (intensity-only probes merged by the oracle, its opposite-strand indel bug EDGE-012, normalisation-based
  grouping, OUT-018 ploidy); indel representations identical, including 144 gtc2vcf cannot resolve
- validate_p6_vcf.py: rule-based categories (--bpm for [I/D]/[D/I] detection)
- docs: validation P6 complete, README status
2026-09-28 13:20:44 +01:00
kemal 0b97462565 P6: gtc-to-vcf (VCF/BCF + indexes), harmonisation, exports (plink, bedgraph, finalreport)
- formats/fasta: plain/BGZF FASTA, .fai/.gzi built identical to htslib, index check, build detection,
  contig aliases, mmap genome cache (IN-060..064)
- formats/bgzf, formats/bcf: pure-Python BGZF, tabix .tbi, BCF 2.2 and CSI (OUT-010, OUT-038, OUT-040)
- harmonize/sites: plus-strand mapping with audit rows, indel resolution from SourceSeq (port of GTCtoVCF,
  Apache-2.0) and left normalisation (OUT-014..021, OUT-030, OUT-042)
- harmonize/combine: OUT-031 genotype combination (ploidy-aware, Infinium I/II signal classes, OUT-037)
- pipeline/vcf: `simurg gtc-to-vcf` (per-sample VCF.gz+tbi, --merge cohort BCF+csi, audit parquet)
- pipeline/export: `simurg export bedgraph|finalreport|plink|locus-summary` (OUT-043, OUT-050..052, OUT-071)
- validation: 36 GSA samples vs gtc2vcf: 23.66 M identical calls, 0 unexplained differences; PLINK identical
  to plink2 on 637,397 variants; bcftools validates VCF/BCF and our indexes answer region queries correctly
- fetch_data.sh p6 (Ensembl GRCh38 with checksum), make_oracle_vcf.sh, validate_p6_vcf.py; K-022/K-023
2026-09-28 12:17:53 +01:00
kemal 80b93e9f45 P5 validated on 2141 1000G samples; control-metric fixes (K-021)
- Omni2.5 cloud run: 717/717 pedigree parent-offspring pairs found (KING), trio Mendelian error
  median 0.031 %, sex vs pedigree 99.6 %, P4 preset re-check 2138/2141 byte-identical (K-015 closed)
- qc/controls: skip bead addresses absent from the IDAT; control missing only if none is on the chip
- QC-008 amended: signal / zero background = +Inf (passes), 0/0 = NaN (flagged); Inf in TSV/report
- plotly bundle stored gzip-compressed (repo grep stays clean: the DOM event name contains the
  forbidden substring)
- docs: validation P5 part 2, summaries, known_unknowns K-021
2026-09-28 10:25:32 +01:00
kemal fc01ab0a8d P5b: self-contained HTML QC report and QC table (simurg qc report)
- pipeline/qc_report: dataset discovery/labels/types (QC-041..043), robust TSV/CSV reading (QC-061),
  join + unmatched warning (QC-040), privacy options (QC-046), xlsx (4 sheets, conditional fills) / csv
  table (QC-045), byte-deterministic HTML with gzip+base64 embedded data (QC-044, QC-047)
- qc/report: template, CSS and JS app — Control Dashboard, Automated QC (table with sort/filter/search/
  paging/column menu/export, metric scatters with threshold lines), Sample Heatmaps (chip/plate,
  continuous/binary), Trend Analysis (grouped summary, Tukey box plots, counts, sex, PCA), thresholds modal
  (apply / download / upload YAML), linked selection across plots and table (QC-048..060)
- bundled plotly.js 2.35.2 (MIT); openpyxl dependency (MIT)
- tests: discovery/labels/errors/determinism/xlsx + headless-browser self-test
- known_unknowns K-020: remaining report deviations
2026-09-28 09:20:18 +01:00
kemal 81a7e262b3 P5: cloud task for cohort QC validation; configurable gp3 throughput; batched contamination map
- scripts/cloud/tasks/p5.sh: P4 preset re-check, genotype (v5) with QC, validate_p5_large.py
- scripts/validate_p5_large.py: identify samples via het pattern vs 1000G VCF; sex, parent-offspring,
  Mendel, contamination/control distributions
- trial.sh: DISK_MBPS / DISK_IOPS
- qc/contamination: bounded map batches (executor.map submits eagerly)
2026-09-28 09:10:02 +01:00
kemal a4c9c2d7d8 P5a: QC — control metrics, thresholds, status, contamination, relatedness, locus summary
- qc/controls: canonical control names (QC-002), bkg/offset (QC-003), 13 genotyping metrics (QC-004/005),
  missing/zero-divisor handling (QC-007/008); controls_raw.tsv / controls_qc.tsv (QC-009/011) incl. compat CSVs
- qc/thresholds + profiles: YAML schema, validation codes, precedence, shipped profiles, auto (QC-020..025)
- qc/status: functional PASS/FAIL/NA, control CLEAR/FLAGGED, flags (QC-030..033)
- qc/contamination: BAF regression; bit-compatible with gtc2vcf BAFregress (36/36, K-017), exact mode too
- qc/relatedness: KING-robust, duplicates, Mendelian errors, pedigree from sheet (IN-018)
- qc/locus_summary: OUT-072, bit-identical to BeadArrayFiles except ChiTest (K-018)
- qc/suggest: suggest-thresholds (QC-026)
- pipeline/qc: disk-backed cohort store; CLI qc controls|samples|contamination|relatedness|loci|suggest-thresholds
- genotype: control tables (OUT-070), het_rate, QC columns, cohort QC
- tests: unit + oracle (BAFregress plugin, BeadArrayFiles locus summary)
2026-09-28 08:35:10 +01:00
kemal f469c97917 P4 validated on 2141 vendor GTCs; legacy AutoCall 1.6.2 preset
- vendor comparison (1000 Genomes Omni2.5, AutoCall 1.6.2.2): every
  computed field identical in 2138/2141 samples; the 3 others trace back
  to the single-locus raw-intensity rule (K-009); sex rule 1 matches all
- AutoCall 1.6.2.x leaves GC10/GC50/#calls/#no-calls at 0 in GTC v3:
  new `--preset autocall-1.6.2` (1.1.2 normalization, sex rule 1, GTC v3,
  legacy statistics); the sample summary keeps the real statistics
- NaN sign canonicalised to the x86 pattern (K-014); report P4.md
2026-09-28 01:21:49 +01:00
kemal 7fb3967c3f P4: unit tests for calling/BAF/LRR/sex/CLI; vendor validation script
- tests/unit/test_genotype.py: fuzzy memberships incl. NaN rules, calls
  at cluster centres, intensity threshold, BAF/LRR interpolation and
  NaN pattern, LogR-dev locus set, sex rules v1/v2, TOP-strand base
  calls, invariant dates, end-to-end `simurg genotype` (summary schema,
  pass-through columns, thread/layout determinism, overwrite refusal)
- tests/oracle/test_p4_genotype.py: 36 GSA samples x 2 versions
  byte-identical to the oracle
- scripts/validate_p4_large.py + cloud task p4: compare every vendor
  GTC field by field and as whole files; reports which sex rule matches
2026-09-28 00:10:47 +01:00
kemal 354567e1d8 P4 (WIP): GenCall-equivalent calling, BAF/LRR, sex, GTC v5, simurg genotype
- genotype/clusters.py: EGT aligned to manifest order (float32 mean R)
- genotype/gencall.py: flair, six-case likelihood, priors/posteriors,
  fuzzy sub-scores and score map, vectorised with the reference's
  float/double evaluation order and libm transcendentals
- genotype/baf_lrr.py, sex.py (rules v1/v2), sample.py (IDAT -> Gtc)
- pipeline/inputs.py shared by ingest and genotype
- pipeline/genotype.py + CLI: GTC per sample, sample_summary.tsv/json
  (OUT-060..063), run manifest
- GTCs from IDATs are byte-identical to the oracle for 36 GSA samples x
  2 normalization versions (NaN sign canonicalised to the x86 pattern)
2026-09-28 00:07:37 +01:00
kemal f7d883826d P3 validated: normalization bit-identical on 2141 vendor GTCs
Normalization 1.1.2 reproduces all 109,191 bins of the 1000 Genomes
Omni2.5 GTCs (AutoCall 1.6.2.2) with 0 ULP difference, reserved fields
included; 1.2.0 and 1.1.2 are bit-identical to the no-FMA gtc2vcf oracle
on 36 GSA samples. AWS run: 33 min, ~$0.15, no resources left.
2026-09-27 23:41:43 +01:00
kemal 7deabd907c P3: vectorised nearest-neighbour search; task-driven cloud trials
- fit.py: exact and binned nearest-neighbour searches use one distance
  matrix per searched bin (first-minimum rule kept); fit 2.6 s -> 1.4 s
  per GSA sample, still bit-identical to the oracle
- scripts/validate_p3_large.py: refits 1.2.0 and 1.1.2 transforms for
  every vendor GTC and reports ULP distances per field
- cloud trial runs scripts/cloud/tasks/<TASK>.sh (TASK=p3 default;
  p1p2 = previous validation + ingest)
2026-09-27 22:57:03 +01:00
kemal c0d556d4af P3: normalization fit (1.2.0 and 1.1.2), apply, polar, simurg normalize
- signal/normalization/numerics.py: vendor percentile, trimmean, MAD sigma
  (old/new), bisquare robust fits, float32/float64 typing and sequential
  sums reproduced exactly (ported from gtc2vcf, MIT)
- signal/normalization/fit.py: outlier removal, offset, rotation, shear,
  scale, nn12/rr12, single-channel model for small Infinium I bins,
  exact and binned nearest-neighbour sweeps; vendor quirks marked QUIRK
- signal/normalization/apply.py: affine transform and theta/R
- pipeline/normalize.py + CLI `simurg normalize`
- oracle GTCs rebuilt with -ffp-contract=off (aarch64 gcc fuses a*b+c by
  default, shifting the last bits); the old FMA oracle differed by <=4 ULP

Result: transforms bit-identical to the oracle for all 36 GSA samples x
48 bins, for both normalization versions (3456/3456 bins, 0 ULP).
2026-09-27 22:51:06 +01:00
kemal c8cd2f741c P2 validated on 2141 vendor GTCs (1000 Genomes Omni2.5)
Raw X/Y and controls assembled from IDATs equal the vendor GTC arrays for
2138/2141 samples x 2.45M loci; 3 single-locus differences are recorded
as K-009 (probable per-channel missing rule in AutoCall 1.6.2.2). Channel
order X=Red confirmed (U1). `simurg ingest` processed all 2141 samples in
19 min (I/O-bound). AWS run cost ~$0.34 and left no resources.
2026-09-27 22:19:28 +01:00
kemal 7b83155d56 P2 validation: oracle test and large-scale raw X/Y check
- tests/oracle/test_p2_raw_xy.py: IDAT-assembled raw X/Y and controls
  equal the arrays in oracle GTCs
- validate_p1_large.py: phase-P2 comparison against vendor GTCs with
  diagnostics (channel swap, bead-count disagreements)
- cloud trial: also runs `simurg ingest` on all samples (time/memory),
  disk 250 GB
2026-09-27 20:57:42 +01:00
kemal 4c0ea5e04f P2 (WIP): raw X/Y assembly, canonical Parquet model, simurg ingest
- signal/raw_xy.py: vectorised Infinium I/II raw X/Y and control assembly
  (plan §7.1) with a per-chip address-map cache; IN-027 coverage gate
- core/schemas.py: Arrow schemas (loci, sample_intensity, controls),
  deterministic atomic Parquet writing
- pipeline/ingest.py + CLI: discovery, sheet, IN-010 naming, overwrite
  check before any work, spawn-based workers, per-sample isolation,
  --resume by input hashes, run manifest, exit codes 0/3/4
- matches gtc2vcf raw X/Y and controls exactly on 18 GSA samples;
  outputs byte-identical across thread counts and layouts
2026-09-27 20:54:06 +01:00
kemal fb1ff51b76 P1 validated on 1000 Genomes Omni2.5 (2141 vendor GTCs + IDATs)
All 2141 vendor-written GTC v3 files match BeadArrayFiles field by field
and re-encode byte-identically; all IDAT pairs parse; the 2.45M-locus BPM,
dbv-7 EGT and CSV manifest agree with BeadArrayFiles and each other.
Run on a self-destructing AWS spot instance (42 min, ~$0.20).
2026-09-27 20:22:35 +01:00
kemal b012f74016 cloud trial: set HOME/PATH for cloud-init, unit tests non-fatal
The first run stopped after 1 minute (HOME unbound under set -u) and
terminated itself as designed.
2026-09-27 19:36:49 +01:00
kemal e786b6b64b Lint fixes in validate_p1_large.py 2026-09-27 19:30:37 +01:00
kemal 6bea6c4930 Add self-destructing AWS trial and large-scale P1 validation script
scripts/cloud/trial.sh launches one spot instance that validates the
parsers on the 1000 Genomes Omni2.5 set (2141 samples, vendor GTCs)
streamed from the public S3 bucket, uploads the summary and terminates
itself. Cost guards: spot price cap, hard time limit (shutdown =
terminate), EXIT-trap shutdown, CloudWatch idle alarm that terminates
from outside, disk deleted with the instance, monthly budget alert, and
an `audit` command listing anything still billable.

scripts/validate_p1_large.py compares every GTC with BeadArrayFiles,
checks byte-identical round-trips, IDAT pairs and BPM/EGT/CSV.
2026-09-27 19:30:27 +01:00
kemal e518ea03d5 P1: file parsers, sample sheet, discovery, CLI inspect
Readers for IDAT v3, BPM (locus entries v4-v8), CSV manifest, EGT
(data blocks 5-9) and GTC v3-v5 (plus a layout-preserving GTC writer);
sample-sheet validation (IN-001..IN-014), IDAT discovery and pairing
(IN-020..IN-022), message catalogue and exit codes (CLI-014/016),
JSON Lines logging (CLI-007/008), atomic writes (CLI-011), run manifest,
and `simurg inspect`.

Validation (docs/validation/P1.md): GTC, BPM and EGT readers agree
exactly with Illumina BeadArrayFiles on GSA v4 / GDA v1 files; GTCs
round-trip byte-identically; BPM and CSV manifests are consistent.
2026-09-27 16:56:39 +01:00
kemal 705138d70b P0: project setup (packaging, CI, licences, data/oracle scripts)
- pyproject (uv/hatch, ruff, mypy strict, pytest markers), MIT LICENSE,
  THIRD_PARTY_NOTICES with verbatim BSD-2 (BeadArrayFiles) and MIT (gtc2vcf)
- CI workflow for Gitea/GitHub Actions: lint, type check, unit tests
- fetch_data.sh p0 now fetches the GSA v4 PGx demo IDATs (the GDA "Demo 36"
  bundle is a GenomeStudio project), unpacks archives, fixes pipefail bug
- make_oracle_gtcs.sh builds bcftools + gtc2vcf and writes oracle GTCs
- plan §17: corrected P0 dataset table
2026-09-27 16:56:39 +01:00
kemal 854cf95029 Add sanitized Phase-0 research archive under docs/research
Research notes R1-R4, documentation-study files A-E (source of the
functional spec's `src:` tags) and PGx coverage CSVs, rewritten to use
neutral naming and current project decisions (plan §2).
2026-09-27 16:27:34 +01:00
kemal 336a6b8226 Add functional specification (416 items) and link it from plan/docs
SIMURGARRAY_FUNCTIONAL_SPEC.md defines externally visible behaviour:
CLI, inputs, genotyping outputs, QC and report, cytogenetics CNV/LOH,
targeted PGx copy number, star alleles, methylation control QC, edge-case
catalogue, performance targets, open questions and a traceability table.
2026-09-27 16:13:40 +01:00
kemal 631c1d030d SimurgArray: master plan, repo scaffolding and data/ref fetch scripts
Phase 0 deliverables: full technical specification and roadmap
(SIMURGARRAY_MASTER_PLAN.md), CLAUDE.md, bootstrap script for reference
open-source repos and tiered dataset fetch script with checksums.
2026-09-27 15:20:06 +01:00