v1.0.0: user guide, examples, changelog
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.
This commit is contained in:
kemal committed
1 parent
a7ea7c3bdf
commit
f64bdf11f2
7 files changed
+299
-22
No files matched your search
@@ -0,0 +1,161 @@
|
||||
# SimurgArray user guide (v1.0)
|
||||
|
||||
SimurgArray turns Illumina Infinium genotyping BeadChip scans (IDAT files) into genotype calls, QC tables and
|
||||
reports, VCF/BCF and PLINK files. It is a personal, non-commercial **research** tool — not for clinical or
|
||||
diagnostic use. Behaviour is specified item by item in `SIMURGARRAY_FUNCTIONAL_SPEC.md` (IDs such as CLI-021 are
|
||||
cited below); the algorithms and their sources are in `SIMURGARRAY_MASTER_PLAN.md`.
|
||||
|
||||
Contents: [Install](#install) · [Inputs](#inputs) · [A complete run](#a-complete-run) · [Commands](#commands) ·
|
||||
[Presets and reproducibility](#presets-and-reproducibility) · [QC](#quality-control) · [Outputs](#outputs) ·
|
||||
[Accuracy](#accuracy) · [Limitations](#limitations) · [Troubleshooting](#troubleshooting)
|
||||
|
||||
## Install
|
||||
|
||||
Python ≥ 3.12. With [uv](https://docs.astral.sh/uv/):
|
||||
|
||||
```bash
|
||||
git clone <repository> simurg && cd simurg
|
||||
uv venv && uv pip install -e . # add ".[dev]" for tests and linters
|
||||
.venv/bin/simurg --version
|
||||
```
|
||||
|
||||
Only pure-Python dependencies are needed (numpy, pyarrow, typer, rich, pyyaml, openpyxl). VCF/BCF, tabix/CSI
|
||||
indexes and PLINK files are written by SimurgArray itself; bcftools/plink2 are not required (they are used only
|
||||
as test oracles).
|
||||
|
||||
## Inputs
|
||||
|
||||
| File | Needed by | Notes |
|
||||
|---|---|---|
|
||||
| **IDAT pairs** (`<barcode>_<RxxCxx>_Grn.idat` + `_Red.idat`, optionally `.gz`) | `genotype`, `ingest`, `qc controls` | `--idats DIR` is searched recursively; the sample ID is `<barcode>_<position>`. Encrypted (v1) IDATs are refused. Duplicate pairs: `--on-duplicate-idat error\|first\|newest` |
|
||||
| **BPM manifest** | all commands | defines loci, bead addresses, normalization bins, positions (its genome build decides the VCF build) |
|
||||
| **EGT cluster file** | `genotype`, `qc loci`, optional elsewhere | must describe the manifest's loci (`E_EGT_MANIFEST_MISMATCH` otherwise); a cluster file of a sibling manifest version works if it covers the loci (a warning is printed) |
|
||||
| **CSV manifest** | `gtc-to-vcf` (indels), `export finalreport` (strands); optional for `genotype` (cross-check) | same product/version as the BPM |
|
||||
| **Reference FASTA** | `gtc-to-vcf` | plain or bgzip; contig names must match the manifest (`1…22, X, Y, MT`) or use `--contig-alias auto\|FILE`. An mmap cache is built next to it on first use (~3 GB for GRCh38) |
|
||||
| **Sample sheet** (CSV) | optional everywhere | GenomeStudio style (`[Header]`, `[Manifests]`, `[Data]`) or a plain table. Key columns `SentrixBarcode_A`, `SentrixPosition_A`; recognised: `Sample_ID`, `Sample_Name`, `Sample_Plate`, `Sample_Well`, `Sample_Group`, `Sex`/`Gender`, `Parent1`, `Parent2`, `Replicate`. Other columns are carried into the QC table. Example: `examples/sample_sheet.csv` |
|
||||
|
||||
`simurg inspect FILE…` prints the header of any of these (add `--json` for scripts).
|
||||
|
||||
## A complete run
|
||||
|
||||
`examples/gsa_demo.sh` runs everything below on the public 36-sample GSA v4 PGx demo
|
||||
(`bash scripts/fetch_data.sh p0 && bash scripts/fetch_data.sh p6` first).
|
||||
|
||||
```bash
|
||||
# 1. IDAT → GTC v5 per sample + sample_summary.tsv + control-probe QC + contamination + relatedness
|
||||
simurg genotype --bpm M.bpm --egt C.egt --csv M.csv --idats idats/ --sheet sheet.csv \
|
||||
--out run/gtc --qc-profile genotyping --threads 0
|
||||
|
||||
# 2. self-contained HTML QC report (+ xlsx/csv QC table)
|
||||
simurg qc report --data run/gtc --out run/qc_report --output-format both
|
||||
|
||||
# 3. plus-strand VCF.gz + .tbi per sample, or one cohort BCF + .csi (--merge)
|
||||
simurg gtc-to-vcf --bpm M.bpm --csv M.csv --fasta GRCh38.fa --gtcs run/gtc --out run/vcf
|
||||
simurg gtc-to-vcf --bpm M.bpm --csv M.csv --fasta GRCh38.fa --gtcs run/gtc --out run/bcf --merge
|
||||
|
||||
# 4. exports
|
||||
simurg export plink --vcfs run/vcf --gtcs run/gtc --out run/plink --format both # .bed/.bim/.fam + .pgen/.pvar/.psam
|
||||
simurg export bedgraph --bpm M.bpm --gtcs run/gtc --out run/tracks --chrom-style ucsc
|
||||
simurg export finalreport --bpm M.bpm --csv M.csv --egt C.egt --gtcs run/gtc --out run/report # --matrix for loci × samples
|
||||
```
|
||||
|
||||
`--threads 0` uses all cores. Every command refuses to overwrite existing outputs unless `--force` is given.
|
||||
|
||||
## Commands
|
||||
|
||||
| Command | Purpose |
|
||||
|---|---|
|
||||
| `simurg inspect FILES… [--json]` | headers / tables of contents of IDAT, BPM, CSV, EGT, GTC, sample sheets (CLI-020) |
|
||||
| `simurg genotype` | IDAT → normalization → GenCall-equivalent calls → GTC; sample summary, controls, cohort QC. Key options: `--cutoff` (GenCall score no-call threshold, 0.15), `--preset`, `--norm-version`, `--sex-version`, `--gtc-version 3\|5`, `--autocall-date`, `--no-controls`, `--no-cohort-qc`, QC threshold options |
|
||||
| `simurg ingest` / `simurg normalize` | the first two stages on their own: raw per-sample intensities (Parquet), then fitted transforms and normalized X/Y/θ/R. `ingest --resume` keeps outputs whose inputs are unchanged |
|
||||
| `simurg qc controls` | control-probe metrics from IDATs without genotyping (`--array-type genotyping\|methylation`) |
|
||||
| `simurg qc samples\|contamination\|relatedness\|loci` | sample QC, BAF-regression contamination, KING-robust kinship (+ Mendel errors from `Parent1/Parent2`), per-locus summary — from existing GTCs |
|
||||
| `simurg qc report` | HTML report (views, linked selection, thresholds editor) and QC table from one or more dataset folders |
|
||||
| `simurg qc suggest-thresholds` | data-driven threshold YAML from a cohort's `controls_qc.tsv` |
|
||||
| `simurg gtc-to-vcf` | harmonised VCF/BCF. Options: `--merge`, `--unsquash-duplicates`, `--merge-snv-indel`, `--auxiliary-loci`, `--info-fields`, `--extra-format`, `--filter-gs`, `--drop-ambiguous`, `--ploidy-policy sex-aware\|diploid`, `--include-loci/--exclude-loci` |
|
||||
| `simurg export plink\|bedgraph\|finalreport\|locus-summary` | downstream formats |
|
||||
| `simurg compare A B` | field-level differential test of GTCs or VCFs (files or folders paired by name): EXACT / NEAR / DISCREPANCY / UNKNOWN with an HTML dashboard; `--fail-on-discrepancy` exits 3 |
|
||||
|
||||
Global logging options on every command: `--quiet`, `--log-level`, `--log-json` (JSON Lines on stderr), `--debug`.
|
||||
Run `simurg COMMAND --help` for the full list.
|
||||
|
||||
## Presets and reproducibility
|
||||
|
||||
* `--preset current` (default): normalization 1.2.0, sex rule 2, GTC v5 — the current vendor software behaviour.
|
||||
* `--preset autocall-1.6.2`: normalization 1.1.2, sex rule 1, GTC v3 with legacy sample statistics — reproduces the
|
||||
public 1000 Genomes Omni2.5 GTCs byte for byte (2138 of 2141; the 3 exceptions are in `docs/known_unknowns.md`, K-009).
|
||||
* Individual settings override the preset (`--norm-version`, `--sex-version`, `--gtc-version`).
|
||||
* The default numerics are bit-exact with the reference implementations (including their float32 and median quirks);
|
||||
where a "clean-math" choice exists it is an explicit option, e.g. `qc samples --median exact`.
|
||||
* Outputs are deterministic for the same inputs and options. `--autocall-date` fixes the one time-stamp written
|
||||
into GTCs. Every run writes `run_manifest.json` (version, git revision, argv, effective options and where each
|
||||
came from, environment, seed) and a log under `logs/`.
|
||||
|
||||
## Quality control
|
||||
|
||||
`genotype` computes per-sample call rate (all / autosomal), LogR deviation, GC10/GC50, intensity percentiles,
|
||||
estimated sex and heterozygosity; control-probe metrics (staining, extension, hybridisation, target removal,
|
||||
stringency, non-specific binding, non-polymorphic, restoration); BAF-regression contamination (≥ 10 samples) and
|
||||
KING-robust relatedness (duplicates, first/second/third degree, unexpected relatives, Mendel errors when the sheet
|
||||
has parents).
|
||||
|
||||
Each sample gets two statuses: `functional_status` (call rate, LogR deviation, …) and `control_status` (control
|
||||
probes), plus a `flags` list (`low_call_rate`, `high_lrr_sd`, `het_outlier`, `duplicate_of:…`, `sex_mismatch`, …).
|
||||
Thresholds come from a **profile** (`--qc-profile genotyping|genotyping_example|genotyping_ex|methylation_epic|methylation_msa|auto`,
|
||||
files in `src/simurg/qc/profiles/`), a YAML file (`--qc-config`), or single overrides
|
||||
(`--threshold-call-rate 0.97`, `--threshold KEY=VALUE`). Without a profile, flags are reported but no sample fails.
|
||||
`simurg qc suggest-thresholds` proposes cohort-specific values.
|
||||
|
||||
## Outputs
|
||||
|
||||
`simurg genotype --out run/gtc`:
|
||||
|
||||
| Path | Content |
|
||||
|---|---|
|
||||
| `samples/<id>/<id>.gtc` | GTC v5 (or v3): calls, scores, raw and normalized intensities, BAF, LRR, transforms, sample statistics — readable by BeadArrayFiles, Picard, gtc2vcf |
|
||||
| `sample_summary.tsv` / `.json` | one row per sample: status, call rates, LogR dev, sex (estimated/reported/check), GC10/50, intensity percentiles, het rate, contamination, control status, flags, sheet columns |
|
||||
| `controls_raw.tsv`, `controls_qc.tsv` | control-probe intensities and derived metrics |
|
||||
| `relatedness.tsv` | related / duplicate pairs: kinship, concordance, degree, expected relation |
|
||||
| `mendel.tsv` | trio Mendelian error counts (when the sheet has `Parent1/Parent2`) |
|
||||
| `run_manifest.json`, `logs/` | provenance |
|
||||
|
||||
`simurg gtc-to-vcf --out run/vcf`: `<id>.snv.vcf.gz` + `.tbi` per sample, or `cohort.snv.bcf` + `.csi` with
|
||||
`--merge`. FORMAT `GT:GS:BAF:LRR` (GS = GenCall score); alleles on the plus strand of the reference; multi-probe
|
||||
designs of one site merged (`--unsquash-duplicates` keeps one record per probe); indels resolved from the CSV
|
||||
`SourceSeq`. `audit/loci.parquet` records the strand decision of every locus, `audit/merges.parquet` every merge.
|
||||
|
||||
Exports: PLINK 1 (`.bed/.bim/.fam`) and PLINK 2 (`.pgen/.pvar/.psam`); gzipped BAF/LRR bedGraph per sample
|
||||
(`--smoothed-lrr K` adds a moving-median track); a long final report (`final_report.tsv.gz`: Top/Forward/Plus
|
||||
alleles, GC and GT scores, LRR, BAF, position) or a genotype matrix. File schemas are OUT-xxx items of the spec.
|
||||
|
||||
## Accuracy
|
||||
|
||||
Validated phase by phase (reports in `docs/validation/`):
|
||||
|
||||
| Check | Result |
|
||||
|---|---|
|
||||
| Normalization and calls vs 2141 vendor GTCs (1000 Genomes Omni2.5, `--preset autocall-1.6.2`) | bit-identical transforms; GTCs byte-identical in 2138/2141 |
|
||||
| GSA v4 demo vs gtc2vcf / GTCtoVCF / plink2 | 0 unexplained differences in GTC fields, VCF records and PLINK files |
|
||||
| Pedigree QC, 2141 Omni2.5 samples | 717/717 parent–offspring pairs; trio Mendel error median 0.031 % |
|
||||
| 360 GSA arrays (GSE165845) vs NYGC 30x WGS | genotype concordance 0.9934; replicate concordance 0.9966 |
|
||||
| NA12878 × 3 vs GIAB HG001 v4.2.1 | concordance 0.9985; non-reference F1 0.997 |
|
||||
|
||||
A nightly differential run (`scripts/nightly.py`) re-checks these against the oracles.
|
||||
|
||||
## Limitations
|
||||
|
||||
* Calls need an existing EGT cluster file; cluster training (v2), copy number / LOH / mosaicism and PGx
|
||||
diplotyping are not in v1.0.
|
||||
* The largest cohort tested is 2141 samples (plan goal ≥ 10 000). Memory grows with samples × loci in the cohort
|
||||
QC steps; on an 8 GB machine keep `--threads` low (≈ 2–3 GB per worker for GSA/Omni-size arrays).
|
||||
* Speed is that of NumPy: genotyping a GSA sample takes ~3 CPU-seconds on a server core (360 samples in 2 min on
|
||||
16 cores) and ~25 s on a Raspberry Pi 5 core; a compiled core is planned (P9).
|
||||
* Methylation arrays: control-probe QC only.
|
||||
|
||||
## Troubleshooting
|
||||
|
||||
Errors and warnings carry a stable code (`E_…`, `W_…`) — all are listed in `docs/messages.md` with their spec item.
|
||||
Add `--debug` for a traceback or `--log-json` for machine-readable logs. Exit codes (CLI-014): 0 success,
|
||||
1 internal error, 2 usage or pre-flight validation error (nothing processed), 3 some samples failed (or
|
||||
`compare --fail-on-discrepancy` found discrepancies), 4 no sample succeeded, 130/143 interrupted. Findings about vendor behaviour that SimurgArray reproduces
|
||||
(or deliberately does not) are in `docs/known_unknowns.md`.
|
||||
Reference in new issue
Block a user