From f11b39b5a307c99e34401b17ee94d72e28114a0e Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 20:38:24 +0100 Subject: [PATCH] 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. --- CHANGELOG.md | 33 ++++++++ README.md | 37 ++++----- docs/user_guide.md | 161 ++++++++++++++++++++++++++++++++++++++ examples/README.md | 9 +++ examples/gsa_demo.sh | 70 +++++++++++++++++ examples/sample_sheet.csv | 9 +++ src/simurg/__init__.py | 2 +- 7 files changed, 299 insertions(+), 22 deletions(-) create mode 100644 CHANGELOG.md create mode 100644 docs/user_guide.md create mode 100644 examples/README.md create mode 100644 examples/gsa_demo.sh create mode 100644 examples/sample_sheet.csv diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..fe71354 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,33 @@ +# Changelog + +## 1.0.0 — 2026-09-28 + +First release: the v1.0 goals of the master plan (§1.1), built and validated phase by phase (P0–P8). + +**Formats (P1–P2).** Readers for IDAT (plain or gzipped; encrypted v1 refused), BPM, CSV manifest, EGT, GTC v3–v5 +and sample sheets; `simurg inspect`; raw X/Y assembly into a canonical Parquet model (`simurg ingest`, resumable). + +**Signal and calling (P3–P4).** Per-bin affine normalization 1.2.0 and 1.1.2; GenCall-equivalent calling from an +EGT with GenCall score; θ/R, BAF, LRR, LogR deviation, GC10/GC50, intensity percentiles, sex (rules 1 and 2); +GTC v5/v3 writer (`simurg genotype`). Presets `current` and `autocall-1.6.2`. + +**QC (P5).** Control-probe metrics, threshold profiles and YAML, two-axis sample status and flags, +BAF-regression contamination, KING-robust relatedness and Mendel errors, locus summary, threshold suggestion, a +self-contained HTML QC report with xlsx/csv table (`simurg qc …`). + +**Outputs (P6).** Plus-strand harmonisation with an audit trail; per-sample VCF.gz + tabix or cohort BCF + CSI +(`GT:GS:BAF:LRR`, optional INFO and intensity fields, merged multi-probe sites, indels from SourceSeq, SNV+indel +merging, auxiliary MNV loci); PLINK 1 `.bed` and PLINK 2 `.pgen`; BAF/LRR bedGraph; final report and genotype matrix. + +**Testing (P7).** `simurg compare` (EXACT/NEAR/DISCREPANCY/UNKNOWN per field and locus class, HTML dashboard) and a +nightly differential run against gtc2vcf, GTCtoVCF and vendor GTCs. + +**Validation (P8).** 360 GSA arrays of 120 1000 Genomes individuals vs NYGC 30x WGS (concordance 0.9934) and +GIAB HG001 (0.9985, non-ref F1 0.997); replicate concordance 0.9966. Vendor-equivalence: 2138/2141 Omni2.5 GTCs +byte-identical, zero unexplained VCF/PLINK differences vs gtc2vcf, GTCtoVCF and plink2. + +Documentation: `docs/user_guide.md`, `examples/`, `docs/validation/P1.md … P8.md`, `docs/known_unknowns.md`, +`docs/messages.md`. + +Known limitations: needs an existing cluster file (training is v2); no CN/LOH/PGx yet; largest tested cohort +2141 samples; pure NumPy speed. diff --git a/README.md b/README.md index 7f70261..134725c 100644 --- a/README.md +++ b/README.md @@ -9,27 +9,21 @@ Personal research toolkit for Illumina Infinium genotyping BeadChip data ## Status -Phases **P0 (setup)**, **P1 (file parsers)**, **P2 (canonical model, raw X/Y)**, **P3 (normalization)** and **P4 (genotype -calling, GTC output)** are complete — see `docs/validation/`. Each is validated on 2141 vendor GTCs of the 1000 Genomes Omni2.5 set -(normalization bit-identical; genotypes/scores identical in 2138/2141 samples). -Phase **P5 (QC)** is complete: control-probe metrics, thresholds/profiles, two-axis status, BAF-regression contamination, -KING-robust relatedness, Mendelian errors, locus summary and a self-contained HTML report (`simurg qc report`). On 2141 -1000 Genomes samples all 717 pedigree parent–offspring pairs are recovered and trio Mendelian error rates are ~0.03 % -(`docs/validation/P5.md`). -Phase **P6 (outputs)** is complete: `simurg gtc-to-vcf` (plus-strand harmonisation with audit, merged multi-probe -sites, per-sample VCF.gz + tabix or a cohort BCF + CSI, optional INFO/intensity fields, SNV+indel merging, MNV -auxiliary loci) and `simurg export plink|bedgraph|finalreport|locus-summary` (PLINK 1 `.bed` and PLINK 2 `.pgen`). -On the GSA demo, genotypes match gtc2vcf and GTCtoVCF with zero unexplained differences and PLINK files match plink2 -(`docs/validation/P6.md`). -Phase **P7 (differential testing)** is complete: `simurg compare` (EXACT/NEAR/DISCREPANCY/UNKNOWN per field and -locus class, every discrepancy listed, HTML dashboard) and a nightly run on the Pi (GSA vs gtc2vcf, Omni2.5 subset vs -vendor GTCs, VCF regression) — `docs/validation/P7.md`. -Phase **P8 (scientific validation)** is complete: on GEO GSE165845 (360 GSA arrays, 120 1000 Genomes individuals × -3 replicates) genotype concordance is 0.9934 vs NYGC 30x WGS, 0.9985 vs GIAB HG001 v4.2.1 (non-ref F1 0.997), and -0.9966 between replicates — `docs/validation/P8.md`. -Available commands: `simurg version`, `simurg inspect [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`, -`simurg qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`, `simurg gtc-to-vcf`, -`simurg export plink|bedgraph|finalreport|locus-summary`, `simurg compare`. +**v1.0.0** (2026-09-28) — see `CHANGELOG.md`. IDAT → normalization → GenCall-equivalent genotype calling → GTC v5 → +sample/control QC with an HTML report → plus-strand VCF/BCF → PLINK 1/2, bedGraph, final report; plus a +differential-testing harness. Validation (`docs/validation/`): + +| Check | Result | +|---|---| +| 2141 vendor GTCs, 1000 Genomes Omni2.5 (`--preset autocall-1.6.2`) | normalization bit-identical; GTCs byte-identical in 2138/2141 | +| GSA demo vs gtc2vcf, GTCtoVCF, plink2 | zero unexplained differences | +| 1000 Genomes pedigrees (Omni2.5) | 717/717 parent–offspring pairs; Mendel error median 0.031 % | +| 360 GSA arrays (GSE165845) vs NYGC 30x WGS / GIAB HG001 | concordance 0.9934 / 0.9985; replicates 0.9966 | + +Commands: `simurg inspect`, `ingest`, `normalize`, `genotype`, `qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`, +`gtc-to-vcf`, `export plink|bedgraph|finalreport|locus-summary`, `compare`. **User guide: `docs/user_guide.md`**; +worked example: `examples/gsa_demo.sh`. Next (plan §18): P9 performance (compiled core), P10 packaging, then v2 +modules (cluster training, CN/LOH, PGx). ## Quick start @@ -44,6 +38,7 @@ bash scripts/fetch_data.sh p6 # GRCh38 FASTA (~0.9 GB download, 3.1 bash scripts/make_oracle_vcf.sh # gtc2vcf oracle VCF for the P6 tests uv run pytest tests/oracle # real-data acceptance tests (add -m "not slow" to skip the long ones) uv run simurg inspect data/p0/gsa_v4_pgx/idats_Demo_36/idats_Demo_36/207895980019/207895980019_R03C01_Grn.idat +bash examples/gsa_demo.sh results/gsa_demo 6 # end-to-end run on 6 demo samples ``` Vendor data files are downloaded for personal research use only and are never committed. See `THIRD_PARTY_NOTICES.md`. diff --git a/docs/user_guide.md b/docs/user_guide.md new file mode 100644 index 0000000..3a9d3e9 --- /dev/null +++ b/docs/user_guide.md @@ -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 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** (`__Grn.idat` + `_Red.idat`, optionally `.gz`) | `genotype`, `ingest`, `qc controls` | `--idats DIR` is searched recursively; the sample ID is `_`. 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//.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`: `.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`. diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000..0f9c037 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,9 @@ +# Examples + +| File | What it shows | +|---|---| +| `gsa_demo.sh` | End-to-end run on the public 36-sample GSA v4 PGx demo: `genotype` (GTC + QC) → `qc report` → `gtc-to-vcf` (per-sample VCF and cohort BCF) → `export plink/bedgraph/finalreport` → `compare`. Needs `scripts/fetch_data.sh p0` and `p6`. `bash examples/gsa_demo.sh results/gsa_demo 6` runs on the first 6 samples only. | +| `sample_sheet.csv` | Sample-sheet layout: key columns `SentrixBarcode_A`/`SentrixPosition_A`, reported sex, trio (`Parent1`/`Parent2`, used for Mendel errors), a replicate, and a free column (`Batch`) that is carried into the QC table. The sample identities in it are illustrative, not the real demo samples. | + +The QC threshold profiles shipped with SimurgArray (`src/simurg/qc/profiles/*.yaml`) are the starting point for a +`--qc-config` file; `simurg qc suggest-thresholds` derives cohort-specific values. See `docs/user_guide.md`. diff --git a/examples/gsa_demo.sh b/examples/gsa_demo.sh new file mode 100644 index 0000000..7dc1223 --- /dev/null +++ b/examples/gsa_demo.sh @@ -0,0 +1,70 @@ +#!/usr/bin/env bash +# End-to-end SimurgArray run on the public GSA v4 PGx demo set (36 samples): +# IDAT → GTC + sample QC → HTML QC report → VCF (per sample and cohort BCF) → PLINK / bedGraph / final report. +# +# Prerequisites (from the repository root): +# uv venv && uv pip install -e . +# bash scripts/fetch_data.sh p0 # demo IDATs + manifest + cluster file (~1.4 GB) +# bash scripts/fetch_data.sh p6 # GRCh38 FASTA for VCF output (~0.9 GB download) +# +# Usage: bash examples/gsa_demo.sh [OUT_DIR] [N_SAMPLES] (defaults: results/gsa_demo, all 36) +# Takes ~10 min for 36 samples on a Raspberry Pi 5 (2 worker processes); much less on a desktop. +set -euo pipefail +cd "$(dirname "$0")/.." +OUT=${1:-results/gsa_demo} +N=${2:-0} +THREADS=${THREADS:-2} +SIMURG=${SIMURG:-.venv/bin/simurg} + +D=data/p0/gsa_v4_pgx +BPM=$D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm +CSV=$D/GSA-PGx-48v4-0_20079540_E2_csv/GSA-PGx-48v4-0_20079540_E2.csv +EGT=$D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt +IDATS=$D/idats_Demo_36/idats_Demo_36 +FASTA=data/p6/GRCh38.primary.fa +for f in "$BPM" "$CSV" "$EGT" "$IDATS" "$FASTA"; do + [[ -e $f ]] || { echo "missing $f — run scripts/fetch_data.sh p0 and p6 first"; exit 1; } +done +mkdir -p "$OUT" + +# Optional subset: symlink the first N IDAT pairs into $OUT/idats (same /__.idat layout). +if [[ $N -gt 0 ]]; then + rm -rf "$OUT/idats" + find "$IDATS" -name '*_Grn.idat' | sort | head -n "$N" | while read -r g; do + b=$(basename "$(dirname "$g")"); mkdir -p "$OUT/idats/$b" + ln -s "$(realpath "$g")" "$OUT/idats/$b/"; ln -s "$(realpath "${g%_Grn.idat}_Red.idat")" "$OUT/idats/$b/" + done + IDATS=$OUT/idats +fi + +echo "== 1. inspect the inputs" +$SIMURG inspect "$BPM" "$EGT" | head -20 + +echo "== 2. genotype: GTC v5 per sample, sample_summary.tsv, control-probe QC, contamination, relatedness" +$SIMURG genotype --bpm "$BPM" --egt "$EGT" --csv "$CSV" --idats "$IDATS" --out "$OUT/gtc" \ + --qc-profile genotyping --threads "$THREADS" --force +column -t -s$'\t' "$OUT/gtc/sample_summary.tsv" | cut -c1-150 | head -8 + +echo "== 3. HTML QC report + QC table" +$SIMURG qc report --data "$OUT/gtc" --out "$OUT/qc_report" --output-format both + +echo "== 4. VCF: one bgzipped, tabix-indexed VCF per sample, then one cohort BCF" +$SIMURG gtc-to-vcf --bpm "$BPM" --csv "$CSV" --fasta "$FASTA" --gtcs "$OUT/gtc" --out "$OUT/vcf" \ + --threads "$THREADS" --force +$SIMURG gtc-to-vcf --bpm "$BPM" --csv "$CSV" --fasta "$FASTA" --gtcs "$OUT/gtc" --out "$OUT/bcf" --merge \ + --threads "$THREADS" --force + +echo "== 5. exports" +$SIMURG export plink --vcfs "$OUT/vcf" --gtcs "$OUT/gtc" --out "$OUT/plink" --format both --force +$SIMURG export bedgraph --bpm "$BPM" --gtcs "$OUT/gtc" --out "$OUT/bedgraph" --chrom-style ucsc --force +$SIMURG export finalreport --bpm "$BPM" --csv "$CSV" --egt "$EGT" --gtcs "$OUT/gtc" --out "$OUT/finalreport" --force + +echo "== 6. differential test: re-call with the older normalization and compare field by field" +$SIMURG genotype --bpm "$BPM" --egt "$EGT" --idats "$IDATS" --out "$OUT/gtc_norm112" --norm-version 1.1.2 \ + --no-controls --no-cohort-qc --threads "$THREADS" --quiet --force +$SIMURG compare "$OUT/gtc" "$OUT/gtc_norm112" --bpm "$BPM" --out "$OUT/compare" --title "norm 1.2.0 vs 1.1.2" \ + --threads "$THREADS" --force || true + +echo +echo "Done. Open $OUT/qc_report/*.html and $OUT/compare/dashboard.html in a browser." +find "$OUT" -maxdepth 2 -type f \( -name '*.html' -o -name '*.bcf' -o -name '*.pgen' -o -name '*.bed' \) | sort diff --git a/examples/sample_sheet.csv b/examples/sample_sheet.csv new file mode 100644 index 0000000..985b3b4 --- /dev/null +++ b/examples/sample_sheet.csv @@ -0,0 +1,9 @@ +[Header] +Investigator Name,Example +Project Name,SimurgArray demo +[Data] +Sample_ID,SentrixBarcode_A,SentrixPosition_A,Sample_Plate,Sample_Well,Sex,Parent1,Parent2,Replicate,Batch +child1,207895980019,R03C01,Plate1,A01,F,father1,mother1,,B1 +father1,207895980019,R03C02,Plate1,B01,M,,,,B1 +mother1,207895980019,R04C01,Plate1,C01,F,,,,B1 +child1_rep,207895980019,R05C01,Plate1,D01,F,,,child1,B2 diff --git a/src/simurg/__init__.py b/src/simurg/__init__.py index e9f6f70..70eae7c 100644 --- a/src/simurg/__init__.py +++ b/src/simurg/__init__.py @@ -4,5 +4,5 @@ Specification: SIMURGARRAY_MASTER_PLAN.md (algorithms, roadmap) and SIMURGARRAY_FUNCTIONAL_SPEC.md (externally visible behaviour, item IDs such as CLI-021). """ -__version__ = "0.1.0.dev0" +__version__ = "1.0.0" PRODUCT_NAME = "SimurgArray"