198 lines
13 KiB
Markdown
198 lines
13 KiB
Markdown
# SimurgArray user guide (v1.1)
|
||
|
||
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
|
||
```
|
||
|
||
**Optional compiled kernels** (≈ 3–4× faster genotyping, 4× faster manifest loading; results are byte-identical):
|
||
|
||
```bash
|
||
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal # Rust (needs a C linker, e.g. gcc)
|
||
uv pip install ./rust/simurg-core
|
||
.venv/bin/python -c "from simurg._core import backend; print(backend())" # rust 1.1.0
|
||
```
|
||
|
||
Without them the same code runs in pure Python/NumPy; `SIMURG_PURE_PYTHON=1` forces that path.
|
||
|
||
**From release files** (`bash scripts/build_dist.sh` writes them to `dist/`: a pure wheel, an sdist and
|
||
`manylinux_2_28` kernel wheels for x86_64 and aarch64, glibc ≥ 2.28):
|
||
|
||
```bash
|
||
python3 -m venv simurg-env && simurg-env/bin/pip install simurg-1.1.0-py3-none-any.whl \
|
||
simurg_core-1.1.0-cp312-abi3-manylinux_2_28_$(uname -m).whl # the second wheel is optional
|
||
```
|
||
|
||
**Docker** (image ≈ 550 MB, kernels included, runs as a non-root user; data stay in the mounted folder):
|
||
|
||
```bash
|
||
docker build -t simurg:1.1.0 .
|
||
docker run --rm --user "$(id -u):$(id -g)" -v "$PWD":/work simurg:1.1.0 \
|
||
genotype --bpm M.bpm --egt C.egt --idats idats --out run/gtc
|
||
```
|
||
|
||
**conda** (recipes in `packaging/conda/`):
|
||
|
||
```bash
|
||
conda build -c conda-forge --output-folder conda-bld packaging/conda/simurg-core packaging/conda/simurg
|
||
conda create -n simurg -c ./conda-bld -c conda-forge simurg simurg-core
|
||
```
|
||
|
||
All four installs give identical results on x86_64 and aarch64 (`docs/validation/P10.md`).
|
||
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: with the compiled kernels a GSA sample takes ~1.1 s on a Raspberry Pi 5 core (pure Python ~3.4 s); see
|
||
`docs/validation/P9.md` for server numbers and thread scaling. Loading a manifest/cluster file is a fixed cost
|
||
per run (GSA ~5 s, Omni2.5 ~15 s on the Pi).
|
||
* 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, 141 the reader of
|
||
stdout went away (e.g. `simurg inspect … | head`; nothing is reported). Findings about vendor behaviour that SimurgArray reproduces
|
||
(or deliberately does not) are in `docs/known_unknowns.md`.
|