noise-outlier exclusion + robustness test (PCN-017/EDGE-026), platform from the CSV manifest (CLI-032)
Also: star-allele calls are checked for one manifest/build per run (PGX-011), duplications at the CN cap are
named xN (PGX-022), and tandem names attach supporting variants to their full copy.
fetch_sources.sh downloads PharmVar per-gene zips; pgx db build adds core alleles for genes CPIC does not
define (normalised against the FASTA, *1 = reference). Masks re-derived with the new variants. GeT-RM:
CYP3A4 15/15, CYP2A6 7/9; all comparable genes 154/168 (CYP1A2/CYP2C8 definitions differ from GeT-RM's).
Sub-region CNs (p5, intron 2, exon 9) resolve full copies and *36/*68/*13 units; candidates carry the
tandem names, background alleles, and CPIC activity sums. GeT-RM: CYP2D6 13/15 top (NA12878 *3/*68+*4 and
NA11832 *1/*68+*4 exact), 15/15 within feasible; all genes 132/144.
- DB build resolves the PGX-007 tag variants from the manifest's probes; coverage over all gene kinds.
- Variant-based genes reported per variant (rs.. REF/ALT | Reference) with the CPIC lowest-activity diplotype;
tagging genes as positive/negative with the population caveat; GSTM1/GSTT1/UGT2B17 *0 from targeted CN.
- Case-insensitive phenotype keys (MT-RNR1). GeT-RM check extended to VKORC1 and the CN genes:
131/144 gene calls (docs/validation/V2_pgx_star.md).
- supporting_variants.tsv, unresolved_variants.tsv, variants_not_covered.tsv with PGX-026 metadata; per-allele
supporting/unresolved sites for both alleles of every candidate; variant_change vs top, distant_alternative,
unexplained_variants, alleles_impacted_by_no_call, alleles_tested.
- simurg pgx coverage: callable / not_covered / masked / equivalent_to per allele and manifest.
- HGVS on GRCh38 RefSeq accessions (NC_...:g.).
CPIC lookups by allele function counts or activity values (duplications N x activity, *5/*0 = 0, total-score
fallback), Indeterminate when feasible solutions disagree or an allele is not in the guideline; per-sample
simurg.pgx/1 JSON and diplotype_annotations.tsv with PGX-026 metadata lines.
- scripts/pgx/fetch_sources.sh: CPIC API tables (CC0) into a dated folder with version and hashes.
- simurg/pgx/db.py: CPIC HGVS locations and values -> left-normalised VCF variants (SNV, del, ins, dup,
repeats, delins, chrM); IUPAC values are core-allele wildcards; alleles without definitions are flagged.
- pipeline/pgx_db.py: deterministic SQLite (content-hash name), probe_map via the SNV VCF site mapper;
GSA v4 PGx: 22 genes, 1354 alleles, 5002 probe links; two builds byte-identical.
The mosaic BAF split is now estimated over heterozygous calls against each probe's expected noise (the
reference heterozygous spread scaled to the sample), inverting the median of |N(+-d, sigma)|. A mosaic loss
lacking the heterozygous calls its PFB predicts gets MF 1 and is promoted (CYTO-020).
Mosaic fraction from the BAF split now subtracts the heterozygous BAF noise (median |BAF-0.5| of unsplit
hets = 0.6745 sigma_b), and a mosaic loss lacking its expected heterozygotes gets MF 1 and is promoted to a
constitutional DEL (CYTO-020). Adds scripts/v2/cn_tune.sh + cn_tune_report.py (GSE165845 rounds on the Pi).
The Pi now boots with cgroup_enable=memory cgroup_memory=1; the MemAvailable guard stays as a
second layer and as the fallback on systems without the memory controller.
scripts/pi_run.sh runs heavy jobs in their own process group and stops them when
MemAvailable drops below a floor (Raspberry Pi OS boots with cgroup_disable=memory, so a
systemd MemoryMax cap is not available); a 360-sample genotype run with cohort QC had
exhausted memory and the OOM killer took down the desktop session.
GSE165845 360 arrays: constitutional DEL/DUP calls 77/74 % in the NYGC SV truth, LOH
reproduced across replicates 74 %; mosaic calls are mostly false positives (12-23 %),
sensitivity for 25 kb-1 Mb events 7-13 % at QUAL >= 20. All samples had been flagged
SampleQuality by the summary logr_dev (cluster-file mismatch); the gate now uses the SD
of the LRR after the batch reference. Validation script: --ignore-filter.
Found by the CN cloud validation: mask_GRCh3{7,8}.bed and cytoBand_GRCh3{7,8}.txt.gz were
never committed nor shipped in the wheel. The ignore rule now covers only the top-level data/.
cn call now records per-chromosome and sample statistics (percent het, LRR stats, BAF
deviations, robust LogR dev, call rate) in the sidecar; cn annotate writes
<sample>.cyto.json (CYTO-033) with ISCN microarray/traditional strings (CYTO-034), bands
and RefSeq gene lists (CYTO-035), LOH summaries, weighted median CN (CYTO-026); annotation
bundle from UCSC cytoBand + ncbiRefSeq (CYTO-036, IN-054). Tests for ISCN strings, gene
overlap, sidecar-missing error and an end-to-end call+annotate.
216 implanted events in 22 synthetic GSA samples (4 replicates per kind x fraction x
size; reference model from 30 other samples): constitutional DEL/DUP 16/16 at 5-20 Mb,
mosaic loss/gain/CN-LOH at 15 % / 5 Mb 12/12.
beta r 0.99998 vs sesame 1.28.1 (max |dbeta| 0.011), detection calls 99.1% identical,
noob within 0.5% median; replicate r 0.997; all 21 control metrics pass. Spec deviations
recorded (5-field control lines, pOOBAH formula, minfi sheet keys). User guide section.
Methylation CSV manifest reader (MET-002/003), raw M/U assembly with bead masking
(MET-005), pOOBAH detection p-values, noob background correction, linear/nonlinear dye
bias, beta and M-values (MET-010..012), the 21 methylation control metrics (MET-006,
bs-conv red-bg set option), per-sample ctrl/cgs tables, summary tsv/xlsx/html with the
detection histogram and a beta PCA, pcs.tsv.gz, controls_raw/controls_qc.tsv; methylation
datasets load in qc report (MET-014); qc controls --array-type methylation; minfi-style
Sentrix_ID/Sentrix_Position sheet keys. Vendor EPIC manifests carry a 5th field on some
control lines: >= 4 fields accepted. Runs on the public minfiDataEPIC set (3 x NA12878).
Wheels (pure + manylinux_2_28 kernels), sdist, Docker image and conda packages built and
installed on fresh Ubuntu 24.04 machines; 225 unit tests pass against the installed
package; the tutorial runs with pip (with and without kernels); GTC/VCF content digests are
identical across pip, Docker and conda on both architectures and equal to the Pi's run.
docs/validation/P10.md; guide: install from wheels, Docker, conda; exit code 141.
Found by the P10 fresh-VM run: `simurg inspect … | head` raised BrokenPipeError, reported as
an internal error, and stopped examples/gsa_demo.sh (set -o pipefail). run_guarded now flushes
stdout inside the guard and maps BrokenPipeError to ExitCode.SIGPIPE (141) without a message.
p10.sh: SKIP_DOCKER/SKIP_CONDA, clearer digest report; p10t.sh runs build+pip+tutorial only.
scripts/build_dist.sh: sdist + pure wheel + manylinux_2_28 kernel wheels (maturin --zig,
both architectures). Dockerfile (multi-stage, kernels included, non-root user).
packaging/conda/{simurg,simurg-core}. scripts/p10/digest.py: platform-independent content
digests of GTC/VCF/summary outputs. Cloud task p10 installs the artifacts on a fresh
Ubuntu VM and runs the tutorial with pip, Docker and conda installs. trial.sh picks the
arm64/amd64 image from the instance type and can run instances in parallel (PARALLEL=1).
Licence files ship in both wheels; Rust crate notices added.
Graviton3 c7g.8xlarge, 360 GSE165845 arrays: 360/360 GTCs identical between the Rust
kernels and the pure-Python paths; 1 worker 78-80 samples/min (pure Python 29.6).
Scaling 1->8 workers efficiency 0.87-0.90; above that the 125 MB/s gp3 disk caps the
run (~620 samples/min) and, on a RAM disk, a fixed per-worker start-up cost limits a
360-sample batch to efficiency 0.48 at 32 workers (K-025). Worker state is pickled once.
docs/validation/P9.md, result JSONs, guide/README/CHANGELOG updates.
Rust ports of the hot paths found by profiling: fused per-locus GenCall calling
and scoring, BAF/LRR, libm maps with CPython math semantics, both nearest-neighbour
searches of the normalization fit, and the BPM LocusEntry parser. Address
coverage is memoised per IDAT layout. SIMURG_PURE_PYTHON=1 forces the pure paths.
On the Pi: genotype_sample 3.4 s -> 1.1 s, 36-sample GSA run 131 s -> 46 s, BPM
read 4x faster; 36 GSA + 6 Omni2.5 GTCs and QC tables byte-identical to the pure
path; tests/unit/test_core_kernels.py checks every kernel on edge cases.
bench/profile_genotype.py, bench/scaling.py, cloud task p9 (thread scaling).
docs/user_guide.md (install, inputs, complete run, commands, presets, QC, outputs,
accuracy, limitations, troubleshooting), examples/gsa_demo.sh (end-to-end run on the
public GSA demo, verified on the Pi), examples/sample_sheet.csv, CHANGELOG.md; version 1.0.0.
360 GSA arrays (120 1000 Genomes individuals x 3 replicates): concordance 0.9934 vs
NYGC 30x (0.9958 on sites in both sets), 0.9985 vs GIAB HG001 v4.2.1 (non-ref F1
0.997), replicate concordance 0.9966; allele-swap sites 0.016 %; 358/360 replicate
pairs detected as duplicates. Plan §16.5 metrics in docs/validation/P8.md (K-024).
The download-only phase (EBI ~15 MB/s) left the CPU idle long enough for the
idle alarm to terminate the first run. Streams are now cut to the manifest
positions (±100 bp) and the 120 individuals while genotyping runs.
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).
- 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
- 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
- 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
- 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
- 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
- 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)
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.
- 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)
- 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).
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.
- 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
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).
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.
- 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
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).
Phase 0 deliverables: full technical specification and roadmap
(SIMURGARRAY_MASTER_PLAN.md), bootstrap script for reference
open-source repos and tiered dataset fetch script with checksums.
kemal