From 51747674ba8fcae25a3f2a0059a3d2fb2960d386 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 14:29:42 +0100 Subject: [PATCH] P6 complete: optional VCF fields, SNV+indel merging, auxiliary MNV loci, PLINK 2 .pgen MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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 --- README.md | 3 +- SIMURGARRAY_FUNCTIONAL_SPEC.md | 6 +- docs/known_unknowns.md | 2 +- docs/validation/P6.md | 12 ++- src/simurg/cli/main.py | 37 ++++++- src/simurg/formats/bcf.py | 12 ++- src/simurg/harmonize/combine.py | 17 +-- src/simurg/harmonize/sites.py | 127 ++++++++++++++++++++++- src/simurg/pipeline/export.py | 90 +++++++++++++--- src/simurg/pipeline/qc.py | 45 ++++++-- src/simurg/pipeline/vcf.py | 178 ++++++++++++++++++++++++++++++-- tests/oracle/test_p6_vcf.py | 71 +++++++++++++ tests/unit/test_vcf.py | 61 +++++++++++ 13 files changed, 611 insertions(+), 50 deletions(-) diff --git a/README.md b/README.md index 0419982..5328132 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,8 @@ KING-robust relatedness, Mendelian errors, locus summary and a self-contained HT 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) and `simurg export plink|bedgraph|finalreport|locus-summary`. +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`). Available commands: `simurg version`, `simurg inspect [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`, diff --git a/SIMURGARRAY_FUNCTIONAL_SPEC.md b/SIMURGARRAY_FUNCTIONAL_SPEC.md index b2219f3..03f33fb 100644 --- a/SIMURGARRAY_FUNCTIONAL_SPEC.md +++ b/SIMURGARRAY_FUNCTIONAL_SPEC.md @@ -1280,8 +1280,10 @@ about **algorithms or numerics**, the plan wins. Each such conflict found during - (a) exactly one ALT (not multi-allelic); - (b) it is an MNV: `len(REF) = len(ALT) ≥ 2` and REF ≠ ALT; - (c) the probe is the **only** assay at the site (no duplicates); - - (d) the first bases, on the plus strand, of REF and ALT match the probe's designed alleles, one each (for - example, `ATGC[AT/CG]GTAA` designed as `ATGC[A/C]NNNN`: A ↔ AT, C ↔ CG), which fixes the allele mapping; + - (d) the bases the probe interrogates match the probe's designed alleles, one each, which fixes the allele + mapping: the **first** base of REF and ALT on the plus strand for plus-strand designs (for example, + `ATGC[AT/CG]GTAA` designed as `ATGC[A/C]NNNN`: A ↔ AT, C ↔ CG), the **last** base for minus-strand designs, + because the probe extends from the other side (as GTCtoVCF); - (e) REF equals the FASTA. The aux record replaces the probe's REF/ALT, and its genotype is mapped through (d). The final report's Top alleles diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index e67f915..71714e6 100644 --- a/docs/known_unknowns.md +++ b/docs/known_unknowns.md @@ -4,7 +4,7 @@ Plan §20 lists the open questions; this file records what the implementation ha | ID | Date | Topic | Finding | Status | |---|---|---|---|---| -| K-023 | 2026-09-28 | P6 options not yet implemented | `--auxiliary-loci` (OUT-035), `--info-fields`/`--extra-format` (OUT-012; GENTRAIN/CLUSTER_SEP would also need the EGT), `--merge-snv-indel` (OUT-021), PLINK 2 `.pgen`. Everything else of CLI-022/CLI-023 works. | open | +| K-023 | 2026-09-28 | P6 optional outputs | Implemented and validated the same day: `--info-fields` (with `--egt`), `--extra-format`, `--merge-snv-indel`, `--auxiliary-loci`, PLINK 2 `.pgen`. Two decisions: OUT-035 rule (d) uses the MNV's last plus-strand base for minus-strand designs (spec amended); `--extra-format` writes NORMX/NORMY missing where the raw intensities are both 0, and THETA/R missing where both normalized intensities are 0 (148 values on 3 GSA samples; gtc2vcf writes 0). | resolved | | K-022 | 2026-09-28 | VCF oracles' limits (GSA v4 PGx, GRCh38) | gtc2vcf leaves 144 indel probes unresolved (`REF=I`/`D`) that SimurgArray resolves from `SourceSeq`; both leave 49 indels unresolved. GTCtoVCF (python-2.7 era) runs under Python 3 only with two local patches (pyvcf3 `Format` field; skip instead of crash on 3 CYP2D6 indels whose `SourceSeq` has the dash on the right, `[ALLELE/-]`, which SimurgArray handles), and needs the alt-contig loci filtered when the FASTA has no alts. Neither oracle applies ploidy (males' X/Y/MT hets stay diploid). | documented | | K-021 | 2026-09-28 | Control metrics on legacy Omni2.5 IDATs | (1) The Omni2.5 manifest lists two bead addresses per control (e.g. 27630314 and 127630314, two bead-pool versions); only one is on a given chip. Absent addresses were NaN and poisoned the bead mean → 1794/2141 samples falsely flagged on both staining metrics. Absent addresses are now skipped. (2) 1290 samples have Biotin (Bgnd) green = 0 exactly (background-subtracted), so staining green is signal/0; QC-008 now defines that as +Inf (pass) and only 0/0 as NaN (flag). | resolved (spec QC-008 amended) | | K-020 | 2026-09-28 | QC report scope (P5b) | Implemented: views A–E, linked selection, thresholds modal (apply/download/upload YAML), chip/plate heatmaps, trend table and Tukey box plots, PCA card, xlsx/csv table, privacy options, `--split-by dataset`. Not yet: client-side xlsx export (CSV only), density mode for > 12 000 points (WebGL scatter is used), `--split-by plate|run`, X/Y median R for the sex card (the summary does not carry them; the card shows the sex check instead), methylation datasets (v2). | open (deviations listed) | diff --git a/docs/validation/P6.md b/docs/validation/P6.md index 6fe7f41..be98fbb 100644 --- a/docs/validation/P6.md +++ b/docs/validation/P6.md @@ -74,7 +74,13 @@ has the dash on the right (patched to skip instead of crash), which SimurgArray (float32), otherwise unchanged (> 500 sites of each kind checked per sample). * EDGE-023: raw LRR bedGraph values equal the GTC values after `%.6g` formatting for every single-probe position. -### Not yet implemented (listed in known_unknowns K-023) +### Optional outputs (completed 2026-09-28) -`--auxiliary-loci` (OUT-035), `--info-fields` and `--extra-format` (OUT-012), `--merge-snv-indel` (OUT-021), -PLINK 2 `.pgen` output. +| Option | Check | Result | +|---|---|---| +| `--info-fields ASSAY_TYPE,NORM_BIN,GENTRAIN,CLUSTER_SEP,ILMN_STRAND,PROBE_IDS` (`--egt` for the cluster scores) | INFO vs gtc2vcf `ASSAY_TYPE`, `GenTrain_Score`, `Cluster_Sep`, 3 samples | identical | +| `--extra-format X,Y,NORMX,NORMY,THETA,R` | vs gtc2vcf FORMAT tags, 1 965 810 probe × sample values | X, Y, NORMX, NORMY identical (within the 6-digit print of `bcftools query`); THETA/R identical except 148 values where both normalized intensities are 0: missing for SimurgArray (plan §7.4, as the caller sees them), 0 for gtc2vcf | +| both, in `--merge` BCF | BCF vs per-sample VCFs, 1 970 064 records × all fields | 0 differences | +| `--merge-snv-indel` (OUT-021, EDGE-013) | 1 GSA sample; `bcftools norm -c e -f` | 677 SNV+indel positions merged (635 796 → 635 119 records), 0 REF mismatches; unit test of the combined genotypes | +| `--auxiliary-loci` (OUT-035, EDGE-015) | unit tests for every rule (a)–(e) and a minus-strand design; 3 real probes (plus and minus designs) redefined as MNVs | records replaced, genotypes preserved, 0 REF mismatches | +| `export plink --format pgen|both` (plan §8.4) | `plink2 --pfile` on our `.pgen` → `.bed` | byte-identical to our own `.bed` (637 397 variants), no plink2 warnings | diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 678a510..9e9fe82 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -397,6 +397,22 @@ def gtc_to_vcf_cmd( merge: Annotated[ bool, typer.Option("--merge", help="One multi-sample BCF (+ .csi) instead of per-sample VCFs.") ] = False, + auxiliary_loci: Annotated[ + Path | None, typer.Option("--auxiliary-loci", help="VCF of true MNV definitions for manifest probes (OUT-035).") + ] = None, + merge_snv_indel: Annotated[ + bool, typer.Option("--merge-snv-indel", help="One record for an SNV and an indel at the same position.") + ] = False, + info_fields: Annotated[ + str, + typer.Option( + "--info-fields", help="Comma list: ASSAY_TYPE,NORM_BIN,GENTRAIN,CLUSTER_SEP,ILMN_STRAND,PROBE_IDS" + ), + ] = "", + extra_format: Annotated[ + str, typer.Option("--extra-format", help="Comma list: X,Y,NORMX,NORMY,THETA,R (single-probe records)") + ] = "", + egt: Annotated[Path | None, typer.Option("--egt", help="Cluster file (INFO GENTRAIN, CLUSTER_SEP).")] = None, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, @@ -417,7 +433,12 @@ def gtc_to_vcf_cmd( out=out, bpm=bpm, fasta=fasta, gtcs=gtcs, sheet=sheet, csv=csv, bgzip=not no_bgzip, genome_cache=genome_cache, exclude_loci=exclude_loci, include_loci=include_loci, unsquash_duplicates=unsquash_duplicates, inf1_nocall_informative=inf1_nocall_informative, drop_ambiguous=drop_ambiguous, filter_gs=filter_gs, - contig_alias=contig_alias, ploidy_policy=ploidy_policy, threads=threads, merge=merge, force=force, + contig_alias=contig_alias, ploidy_policy=ploidy_policy, threads=threads, merge=merge, + merge_snv_indel=merge_snv_indel, + auxiliary_loci=auxiliary_loci, + force=force, + info_fields=[x.strip().upper() for x in info_fields.split(",") if x.strip()], + extra_format=[x.strip().upper() for x in extra_format.split(",") if x.strip()], egt=egt, argv=sys.argv[1:], ) # fmt: skip run_guarded(rep, lambda: run_gtc_to_vcf(opts, rep)) @@ -485,6 +506,7 @@ def export_plink_cmd( sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (reported sex).")] = None, gtcs: Annotated[Path | None, typer.Option("--gtcs", help="GTC folder (estimated sex).")] = None, multiallelic: Annotated[str, typer.Option("--multiallelic", help="split|drop")] = "split", + fmt: Annotated[str, typer.Option("--format", help="bed (PLINK 1), pgen (PLINK 2) or both")] = "bed", prefix: Annotated[str, typer.Option("--prefix", help="Output file prefix.")] = "cohort", force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, @@ -495,8 +517,17 @@ def export_plink_cmd( from simurg.pipeline.export import PlinkOptions, run_export_plink rep = _plain_reporter("export plink", log_json, debug, quiet) - opts = PlinkOptions(out=out, vcfs=vcfs, sheet=sheet, gtcs=gtcs, multiallelic=multiallelic, prefix=prefix, - force=force, argv=sys.argv[1:]) # fmt: skip + opts = PlinkOptions( + out=out, + vcfs=vcfs, + sheet=sheet, + gtcs=gtcs, + multiallelic=multiallelic, + format=fmt, + prefix=prefix, + force=force, + argv=sys.argv[1:], + ) run_guarded(rep, lambda: run_export_plink(opts, rep)) diff --git a/src/simurg/formats/bcf.py b/src/simurg/formats/bcf.py index f63f4f0..d56f459 100644 --- a/src/simurg/formats/bcf.py +++ b/src/simurg/formats/bcf.py @@ -71,17 +71,21 @@ def encode_record( filter_id: int, n_sample: int, fmt: Sequence[tuple[int, bytes, bytes]], + info: Sequence[tuple[int, bytes]] = (), ) -> bytes: - """``fmt``: (format dictionary index, type-descriptor bytes, per-sample payload) in FORMAT order.""" + """``fmt``: (format dictionary index, type-descriptor bytes, per-sample payload) in FORMAT order; + ``info``: (dictionary index, typed value) pairs.""" shared = bytearray() shared += struct.pack(" tuple[dict[str, int], dict[str, int name = line[13:].split(",", 1)[0].split(">", 1)[0] contigs.setdefault(name, len(contigs)) return strings, contigs + + +def typed_floats(values: Sequence[float]) -> bytes: + return _type_byte(len(values), FLOAT) + float_bytes(values) diff --git a/src/simurg/harmonize/combine.py b/src/simurg/harmonize/combine.py index be20d0d..fbb1347 100644 --- a/src/simurg/harmonize/combine.py +++ b/src/simurg/harmonize/combine.py @@ -22,13 +22,16 @@ CALL_TEXT = {NC: "NC", AA: "AA", AB: "AB", BB: "BB"} def signal_class(p: ProbeSite, allele: str) -> str | None: """OUT-031 step 1.""" a, b = p.alleles - if p.kind == "snv" and p.infinium == 2 and len(allele) == 1: - grp = RED if allele in RED else GREEN - if a in grp: - return "A" - if b in grp: - return "B" - return None + if p.kind == "snv": # SNV probes read the single base at POS; alleles may be extended by --merge-snv-indel + z, a0, b0 = allele[:1], a[:1], b[:1] + if p.infinium == 2: + grp = RED if z in RED else GREEN + if a0 in grp: + return "A" + if b0 in grp: + return "B" + return None + return "A" if z == a0 else "B" if z == b0 else None if allele == a: return "A" if allele == b: diff --git a/src/simurg/harmonize/sites.py b/src/simurg/harmonize/sites.py index b9fa189..df1ad62 100644 --- a/src/simurg/harmonize/sites.py +++ b/src/simurg/harmonize/sites.py @@ -78,7 +78,7 @@ class Site: @property def id(self) -> str: - return (";" if self.kind == "indel" else ",").join(p.name for p in self.probes) + return ("," if self.kind == "snv" else ";").join(p.name for p in self.probes) def allele_index(self, allele: str) -> int: return 0 if allele == self.ref else 1 + self.alts.index(allele) @@ -285,7 +285,13 @@ def map_locus( strand == "-", audit.reason), audit # fmt: skip -def group_sites(probes: Sequence[ProbeSite], contig_order: dict[str, int], *, unsquash: bool = False) -> list[Site]: +def group_sites( + probes: Sequence[ProbeSite], + contig_order: dict[str, int], + *, + unsquash: bool = False, + merge_snv_indel: bool = False, +) -> list[Site]: """OUT-030/OUT-033/OUT-034: probes sharing CHROM/POS/REF (and, for indels, the ALT) form one site; ALT alleles in order of first appearance over BPM order; records sorted by contig (.fai order), POS, REF, ALT, ID.""" groups: dict[tuple[str, int, str, str, str], list[ProbeSite]] = {} @@ -303,5 +309,122 @@ def group_sites(probes: Sequence[ProbeSite], contig_order: dict[str, int], *, un sites += [Site(contig, pos, ref, [x for x in p.alleles if x != ref], [p], kind) for p in ps] else: sites.append(Site(contig, pos, ref, alts, ps, kind)) + if merge_snv_indel: + sites = _merge_snv_indel(sites) sites.sort(key=lambda s: (contig_order.get(s.contig, 1 << 30), s.pos, s.ref, ",".join(s.alts), s.id)) return sites + + +def _merge_snv_indel(sites: list[Site]) -> list[Site]: + """OUT-021 ``--merge-snv-indel``: SNV and indel sites at the same CHROM/POS become one multi-allelic record. + Every allele is extended to the longest REF (VCF representation); SNV probes still read the base at POS.""" + from dataclasses import replace + + by_pos: dict[tuple[str, int], list[Site]] = {} + for st in sites: + by_pos.setdefault((st.contig, st.pos), []).append(st) + out: list[Site] = [] + for group in by_pos.values(): + kinds = {st.kind for st in group} + longest = max((st.ref for st in group), key=len) + if not ({"snv", "indel"} <= kinds) or not all(longest.startswith(st.ref) for st in group): + out += group + continue + probes = [] + for st in group: + suffix = longest[len(st.ref) :] + probes += [ + replace(p, ref=longest, alleles=(p.alleles[0] + suffix, p.alleles[1] + suffix)) for p in st.probes + ] + probes.sort(key=lambda p: p.index) + alts: list[str] = [] + for p in probes: + alts += [x for x in p.alleles if x != longest and x not in alts] + out.append(Site(group[0].contig, group[0].pos, longest, alts, probes, "mixed")) + return out + + +# ------------------------------------------------------------------------------------------------ auxiliary loci +@dataclass +class AuxRecord: + name: str + chrom: str + pos: int + ref: str + alts: list[str] + + +def apply_auxiliary( + probes: list[ProbeSite], + aux: Sequence[AuxRecord], + fetch: Fetch, + contig_of: Callable[[str], str | None], + strand_of: Callable[[int], str], +) -> tuple[list[ProbeSite], list[str]]: + """OUT-035: replace the plus-strand site of MNV probes by the true variant. Returns (probes, errors); every + violated rule is one error message (E_AUX_INVALID). ``strand_of(locus index)`` gives the RefStrand. + + Rule (d): the base the probe interrogates is the MNV's first base on the plus strand for plus-strand designs and + its last base for minus-strand designs (the probe extends from the other side) [OSS: GTCtoVCF, Apache-2.0].""" + from dataclasses import replace + + by_name = {p.name: k for k, p in enumerate(probes)} + at_site: dict[tuple[str, int], int] = {} + for p in probes: + at_site[(p.contig, p.pos)] = at_site.get((p.contig, p.pos), 0) + 1 + out = list(probes) + errors: list[str] = [] + for r in aux: + where = f"auxiliary record {r.name} ({r.chrom}:{r.pos})" + k = by_name.get(r.name) + if k is None: + errors.append(f"{where}: not a probe of the manifest (or the probe was skipped)") + continue + p = probes[k] + if len(r.alts) != 1: + errors.append(f"{where}: (a) exactly one ALT allele is required") + continue + alt = r.alts[0] + if not (len(r.ref) == len(alt) >= 2 and r.ref != alt): + errors.append(f"{where}: (b) not a multi-nucleotide variant (len(REF) = len(ALT) >= 2, REF != ALT)") + continue + if at_site[(p.contig, p.pos)] > 1: + errors.append(f"{where}: (c) the probe is not the only assay at its site") + continue + minus = strand_of(p.index) == "-" + + def pick(x: str, minus: bool = minus) -> str: + return x[-1] if minus else x[0] + + designed = [a[-1] if minus else a[0] for a in p.alleles] if p.kind == "snv" else [] + mapping = {} + for d in designed: + hits = [x for x in (r.ref, alt) if pick(x) == d] + if len(hits) == 1: + mapping[d] = hits[0] + if len(mapping) != 2 or mapping[designed[0]] == mapping[designed[1]]: + errors.append(f"{where}: (d) REF/ALT do not match the probe's alleles {p.alleles[0]}/{p.alleles[1]}") + continue + contig = contig_of(r.chrom) + if contig is None or fetch(contig, r.pos - 1, r.pos - 1 + len(r.ref)) != r.ref.upper(): + errors.append(f"{where}: (e) REF does not match the reference genome") + continue + out[k] = replace( + p, contig=contig, pos=r.pos, ref=r.ref.upper(), alleles=(mapping[designed[0]], mapping[designed[1]]), + reason="aux_locus", + ) # fmt: skip + return out, errors + + +def read_auxiliary(path: str) -> list[AuxRecord]: + import gzip + + opener = gzip.open if path.endswith(".gz") else open + out = [] + with opener(path, "rt") as fh: + for line in fh: + if line.startswith("#") or not line.strip(): + continue + f = line.rstrip("\n").split("\t") + out.append(AuxRecord(f[2], f[0], int(f[1]), f[3].upper(), [a.upper() for a in f[4].split(",")])) + return out diff --git a/src/simurg/pipeline/export.py b/src/simurg/pipeline/export.py index b0acb58..74ad017 100644 --- a/src/simurg/pipeline/export.py +++ b/src/simurg/pipeline/export.py @@ -8,6 +8,7 @@ import gzip import io import math import re +import struct from collections.abc import Iterator, Sequence from dataclasses import dataclass, field from pathlib import Path @@ -348,21 +349,38 @@ class PlinkOptions: sheet: Path | None = None gtcs: Path | None = None multiallelic: str = "split" + format: str = "bed" prefix: str = "cohort" force: bool = False argv: list[str] = field(default_factory=list) +PGEN_FROM_BED = bytes([2, 3, 1, 0]) # PLINK 1 code (A1 = ALT) → PLINK 2 hardcall (0 ref/ref, 1 het, 2 alt/alt, 3 NA) + + +def _pack(codes: Sequence[int], nbytes: int) -> bytes: + buf = bytearray(nbytes) + for j, c in enumerate(codes): + buf[j >> 2] |= c << (2 * (j & 3)) + return bytes(buf) + + def run_export_plink(opts: PlinkOptions, rep: Reporter) -> int: - """OUT-043: SimurgArray per-sample VCFs (identical record order) → PLINK 1 .bed/.bim/.fam.""" + """OUT-043 / plan §8.4: SimurgArray per-sample VCFs (identical record order) → PLINK 1 .bed/.bim/.fam and/or + PLINK 2 .pgen/.pvar/.psam (fixed-width hardcall storage, mode 0x02 [DOC: pgen format]). Multi-allelic sites are + split (``bcftools norm -m-``) or dropped.""" if opts.multiallelic not in ("split", "drop"): raise SimurgError("E_CONFIG_KEY", "--multiallelic must be split or drop") + if opts.format not in ("bed", "pgen", "both"): + raise SimurgError("E_CONFIG_KEY", "--format must be bed, pgen or both") rm = RunManifest(command="export plink", argv=opts.argv) vcfs = sorted(p for p in opts.vcfs.glob("*.snv.vcf*") if not p.name.endswith((".tbi", ".csi"))) if not vcfs: raise SimurgError("E_NO_INPUT", f"no *.snv.vcf(.gz) files in {opts.vcfs}") base = opts.out / opts.prefix - check_outputs_absent([base.with_suffix(x) for x in (".bed", ".bim", ".fam")], force=opts.force) + want_bed, want_pgen = opts.format in ("bed", "both"), opts.format in ("pgen", "both") + targets = ([".bed", ".bim", ".fam"] if want_bed else []) + ([".pgen", ".pvar", ".psam"] if want_pgen else []) + check_outputs_absent([base.with_suffix(x) for x in targets], force=opts.force) names, sex = [], {} for v in vcfs: opener = gzip.open if str(v).endswith(".gz") else open @@ -388,12 +406,23 @@ def run_export_plink(opts: PlinkOptions, rep: Reporter) -> int: nbytes = (n + 3) // 4 opts.out.mkdir(parents=True, exist_ok=True) n_var = dropped = 0 + pvar_header = "".join( + ln for ln in _vcf_header(vcfs[0]) if ln.startswith("##contig") or ln.startswith("##fileformat") + ) with ( - atomic_open(base.with_suffix(".bed"), "wb") as bed_raw, - atomic_open(base.with_suffix(".bim"), "w", "utf-8") as bim_raw, + atomic_open(base.with_suffix(".bed"), "wb") if want_bed else _Null() as bed_raw, + atomic_open(base.with_suffix(".bim"), "w", "utf-8") if want_bed else _Null() as bim_raw, + atomic_open(base.with_suffix(".pgen"), "wb") if want_pgen else _Null() as pgen_raw, + atomic_open(base.with_suffix(".pvar"), "w", "utf-8") if want_pgen else _Null() as pvar_raw, ): bed, bim = cast(IO[bytes], bed_raw), cast(IO[str], bim_raw) - bed.write(PLINK_MAGIC) + pgen, pvar = cast(IO[bytes], pgen_raw), cast(IO[str], pvar_raw) + if want_bed: + bed.write(PLINK_MAGIC) + if want_pgen: + # header byte 12: nonref-flag storage 1 = every REF is the reference base (from the FASTA) + pgen.write(b"l\x1b\x02" + struct.pack(" int: gts = [r[9].split(":", 1)[0] for r in recs] for k, alt in enumerate(alts, start=1): codes = [bed_code(g, k) for g in gts] - buf = bytearray(nbytes) - for j, c in enumerate(codes): - buf[j >> 2] |= c << (2 * (j & 3)) - bed.write(bytes(buf)) vid = r0[2] if len(alts) == 1 else f"{r0[2]}:{alt}" - bim.write(f"{r0[0]}\t{vid}\t0\t{r0[1]}\t{alt}\t{r0[3]}\n") + if want_bed: + bed.write(_pack(codes, nbytes)) + bim.write(f"{r0[0]}\t{vid}\t0\t{r0[1]}\t{alt}\t{r0[3]}\n") + if want_pgen: + pgen.write(_pack([PGEN_FROM_BED[c] for c in codes], nbytes)) + pvar.write(f"{r0[0]}\t{r0[1]}\t{vid}\t{r0[3]}\t{alt}\n") n_var += 1 - fam = "".join(f"{s}\t{s}\t0\t0\t{ {'M': 1, 'F': 2}.get(sex.get(s, 'U'), 0) }\t-9\n" for s in names) - with atomic_open(base.with_suffix(".fam"), "w", "utf-8") as fam_raw: - cast(IO[str], fam_raw).write(fam) + if want_pgen: + pgen.seek(3) + pgen.write(struct.pack(" None: + return None + + def __exit__(self, *a: object) -> None: + return None + + +def _vcf_header(path: Path) -> list[str]: + opener = gzip.open if str(path).endswith(".gz") else open + out = [] + with opener(path, "rt") as fh: + for line in fh: + if not line.startswith("##"): + break + out.append(line) + return out diff --git a/src/simurg/pipeline/qc.py b/src/simurg/pipeline/qc.py index fbb1d60..44ee319 100644 --- a/src/simurg/pipeline/qc.py +++ b/src/simurg/pipeline/qc.py @@ -308,27 +308,45 @@ class SampleInfo: class CohortStore: """Disk-backed (samples × loci) matrices: genotypes (uint8), BAF and GenCall scores (float32, optional).""" - def __init__(self, root: Path, n_samples: int, n_loci: int, *, baf: bool, scores: bool, lrr: bool = False) -> None: + def __init__( + self, + root: Path, + n_samples: int, + n_loci: int, + *, + baf: bool, + scores: bool, + lrr: bool = False, + intensities: bool = False, + ) -> None: self.root, self.shape = root, (n_samples, n_loci) root.mkdir(parents=True, exist_ok=True) self.gt = np.lib.format.open_memmap(root / "gt.npy", "w+", np.uint8, self.shape) self.baf = np.lib.format.open_memmap(root / "baf.npy", "w+", np.float32, self.shape) if baf else None self.score = np.lib.format.open_memmap(root / "score.npy", "w+", np.float32, self.shape) if scores else None self.lrr = np.lib.format.open_memmap(root / "lrr.npy", "w+", np.float32, self.shape) if lrr else None + # raw and normalized intensities (needs root/lookup.npy: normalization-transform index per locus) + self.rx = self.ry = self.nx = self.ny = None + if intensities: + self.rx = np.lib.format.open_memmap(root / "rx.npy", "w+", np.uint16, self.shape) + self.ry = np.lib.format.open_memmap(root / "ry.npy", "w+", np.uint16, self.shape) + self.nx = np.lib.format.open_memmap(root / "nx.npy", "w+", np.float32, self.shape) + self.ny = np.lib.format.open_memmap(root / "ny.npy", "w+", np.float32, self.shape) def close(self) -> None: - for name in ("gt", "baf", "score", "lrr"): + for name in ("gt", "baf", "score", "lrr", "rx", "ry", "nx", "ny"): a = getattr(self, name) if a is not None: a.flush() self.gt = self.baf = self.score = self.lrr = None # type: ignore[assignment] + self.rx = self.ry = self.nx = self.ny = None shutil.rmtree(self.root, ignore_errors=True) -def _load_one(task: tuple[int, str, str, str, bool, bool, bool]) -> SampleInfo: +def _load_one(task: tuple[int, str, str, str, bool, bool, bool, bool]) -> SampleInfo: from simurg.formats.gtc import read_gtc - i, sid, path, root, want_baf, want_scores, want_lrr = task + i, sid, path, root, want_baf, want_scores, want_lrr, want_int = task info = SampleInfo(sid, path) try: g = read_gtc(path) @@ -356,6 +374,15 @@ def _load_one(task: tuple[int, str, str, str, bool, bool, bool]) -> SampleInfo: sc = np.load(Path(root) / "score.npy", mmap_mode="r+") sc[i] = g.genotype_scores sc.flush() + if want_int: + from simurg.signal.normalization.apply import normalize + + lookup = np.load(Path(root) / "lookup.npy") + nx, ny = normalize(g.raw_x, g.raw_y, lookup, g.normalization_transforms) + for name, arr in (("rx", g.raw_x), ("ry", g.raw_y), ("nx", nx), ("ny", ny)): + mm = np.load(Path(root) / f"{name}.npy", mmap_mode="r+") + mm[i] = arr + mm.flush() except SimurgError as e: info.error = f"{e.code}: {e.message}" except Exception as e: @@ -372,11 +399,15 @@ def load_cohort( scores: bool, workers: int, lrr: bool = False, + lookup: npt.NDArray[np.int32] | None = None, ) -> tuple[CohortStore, list[SampleInfo]]: - """Read GTCs (sample_id, path) into a :class:`CohortStore` with a process pool.""" - store = CohortStore(tmp, len(samples), n_loci, baf=baf, scores=scores, lrr=lrr) + """Read GTCs (sample_id, path) into a :class:`CohortStore` with a process pool. With ``lookup`` (the manifest's + normalization-transform index per locus) raw and normalized intensities are stored too.""" + store = CohortStore(tmp, len(samples), n_loci, baf=baf, scores=scores, lrr=lrr, intensities=lookup is not None) + if lookup is not None: + np.save(tmp / "lookup.npy", lookup) store.gt.flush() - tasks = [(i, sid, p, str(tmp), baf, scores, lrr) for i, (sid, p) in enumerate(samples)] + tasks = [(i, sid, p, str(tmp), baf, scores, lrr, lookup is not None) for i, (sid, p) in enumerate(samples)] if workers <= 1 or len(tasks) <= 1: infos = [_load_one(t) for t in tasks] else: diff --git a/src/simurg/pipeline/vcf.py b/src/simurg/pipeline/vcf.py index 80fc012..d7dd37a 100644 --- a/src/simurg/pipeline/vcf.py +++ b/src/simurg/pipeline/vcf.py @@ -57,6 +57,11 @@ class VcfOptions: ploidy_policy: str = "sex-aware" threads: int = 0 merge: bool = False + merge_snv_indel: bool = False + auxiliary_loci: Path | None = None + info_fields: list[str] = field(default_factory=list) + extra_format: list[str] = field(default_factory=list) + egt: Path | None = None force: bool = False argv: list[str] = field(default_factory=list) @@ -82,11 +87,33 @@ class SiteTable: site_id: list[str] = field(default_factory=list) site_alleles: list[list[str]] = field(default_factory=list) single_ab: npt.NDArray[np.int16] | None = None # (n_single, 2) allele indices of the probe's A and B alleles + site_info: list[str] = field(default_factory=list) # INFO text per site ('.' without --info-fields) + site_info_typed: list[list[tuple[str, str, list[Any]]]] = field(default_factory=list) # (key, kind, values) + extra_format: list[str] = field(default_factory=list) + single_lookup: npt.NDArray[np.int32] | None = None # normalization-transform index of single-probe loci inf1_nc_informative: bool = False ploidy_policy: str = "sex-aware" filter_gs: float | None = None +INFO_DEFS = { # OUT-012, one value per contributing probe (BPM order) + "ASSAY_TYPE": ("Integer", "Assay type per probe: 0 Infinium II, 1 Infinium I red, 2 Infinium I green"), + "NORM_BIN": ("Integer", "Normalization bin per probe"), + "GENTRAIN": ("Float", "GenTrain score per probe (cluster file)"), + "CLUSTER_SEP": ("Float", "Cluster separation score per probe (cluster file)"), + "ILMN_STRAND": ("String", "Design strand per probe (TOP, BOT, PLUS, MINUS)"), + "PROBE_IDS": ("String", "Names of the contributing probes"), +} +EXTRA_DEFS = { # plan §8.3; single-probe records only + "X": ("Integer", "Raw intensity, X (red) channel"), + "Y": ("Integer", "Raw intensity, Y (green) channel"), + "NORMX": ("Float", "Normalized X intensity"), + "NORMY": ("Float", "Normalized Y intensity"), + "THETA": ("Float", "Normalized theta = (2/pi) atan2(NORMY, NORMX)"), + "R": ("Float", "Normalized R = NORMX + NORMY"), +} + + def _site_class(contig: str, pos: int, par_flag: bool, build: str) -> int: c = contig[3:] if contig.lower().startswith("chr") else contig if c in ("MT", "M"): @@ -196,8 +223,21 @@ def build_sites( probes.append(p) for code, count in sorted(skips.items()): rep.warn(code, f"{count} locus/loci: {code}") + if opts.auxiliary_loci is not None: + from simurg.harmonize.sites import apply_auxiliary, read_auxiliary + + aux = read_auxiliary(str(opts.auxiliary_loci)) + strands = mv.ref_strands + probes, errors = apply_auxiliary(probes, aux, genome.fetch, cmap.get, lambda i: (strands[i] or "").strip()) + if errors: + raise SimurgError("E_AUX_INVALID", f"{len(errors)} invalid auxiliary record(s):\n " + "\n ".join(errors)) + rep.info(f"{len(aux)} auxiliary locus definition(s) applied") + names = {r.name for r in aux} + for a_row in audit: + if a_row.locus in names: + a_row.reason = "aux_locus" order = {e.name: k for k, e in enumerate(genome.entries)} - sites = group_sites(probes, order, unsquash=opts.unsquash_duplicates) + sites = group_sites(probes, order, unsquash=opts.unsquash_duplicates, merge_snv_indel=opts.merge_snv_indel) used = {s.contig for s in sites} contigs = [e.name for e in genome.entries if e.name in used] # static per-site text and single-probe genotype tables @@ -222,6 +262,7 @@ def build_sites( else: merged_sites.append(s) merged_rows.append(k) + site_info, site_info_typed = _site_info(opts, bpm, sites) header = vcf_header(opts, genome, build, contigs, rm) table = SiteTable( prefix=prefix, @@ -240,6 +281,10 @@ def build_sites( site_id=[s.id for s in sites], site_alleles=[[s.ref, *s.alts] for s in sites], single_ab=np.array(single_ab, np.int16).reshape(-1, 2), + site_info=site_info, + site_info_typed=site_info_typed, + extra_format=list(opts.extra_format), + single_lookup=bpm.normalization_lookups()[np.array(single_locus, np.int64)] if opts.extra_format else None, inf1_nc_informative=opts.inf1_nocall_informative, ploidy_policy=opts.ploidy_policy, filter_gs=opts.filter_gs, @@ -251,6 +296,49 @@ def _sha(path: Path | None) -> str: return sha256_file(path) if path else "" +def _site_info( + opts: VcfOptions, bpm: Any, sites: list[Site] +) -> tuple[list[str], list[list[tuple[str, str, list[Any]]]]]: + """Static INFO per site for --info-fields (OUT-012).""" + if not opts.info_fields: + return ["."] * len(sites), [[] for _ in sites] + gentrain = sep = None + if {"GENTRAIN", "CLUSTER_SEP"} & set(opts.info_fields): + from simurg.formats.egt import read_egt + + assert opts.egt is not None + egt = read_egt(opts.egt) + eidx = egt.name_index() + rows = np.array([eidx.get(n, -1) for n in bpm.names], np.int64) + gentrain = np.where(rows >= 0, egt.gentrain_score[np.maximum(rows, 0)], np.nan) + sep = np.where(rows >= 0, egt.cluster_sep[np.maximum(rows, 0)], np.nan) + bins = bpm.norm_bin_ids() + texts, typed = [], [] + for s in sites: + idx = [p.index for p in s.probes] + vals: dict[str, list[Any]] = { + "ASSAY_TYPE": [int(bpm.assay_types[i]) for i in idx], + "NORM_BIN": [int(bins[i]) for i in idx], + "GENTRAIN": [float(gentrain[i]) for i in idx] if gentrain is not None else [], + "CLUSTER_SEP": [float(sep[i]) for i in idx] if sep is not None else [], + "ILMN_STRAND": [bpm.ilmn_strands[i] or "." for i in idx], + "PROBE_IDS": [p.name for p in s.probes], + } + parts, tv = [], [] + for key in opts.info_fields: + v = vals[key] + kind = INFO_DEFS[key][0] + text = ",".join( + "." if (isinstance(x, float) and math.isnan(x)) else (f"{x:.6g}" if isinstance(x, float) else str(x)) + for x in v + ) + parts.append(f"{key}={text}") + tv.append((key, kind, v)) + texts.append(";".join(parts)) + typed.append(tv) + return texts, typed + + def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm: RunManifest) -> str: """OUT-011 (without the ##SAMPLE and #CHROM lines, which are per sample).""" @@ -271,14 +359,18 @@ def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm "##simurg_options=" f"unsquash_duplicates={str(opts.unsquash_duplicates).lower()};exclude_loci={lst(opts.exclude_loci)};" f"include_loci={lst(opts.include_loci)};inf1_nocall_informative={str(opts.inf1_nocall_informative).lower()};" - f"auxiliary_loci=none;ploidy_policy={opts.ploidy_policy};contig_alias={opts.contig_alias}" + f"auxiliary_loci={lst(opts.auxiliary_loci)};ploidy_policy={opts.ploidy_policy};contig_alias={opts.contig_alias}" + (";drop_ambiguous=true" if opts.drop_ambiguous else "") + + (";merge_snv_indel=true" if opts.merge_snv_indel else "") + (";unknown_sex_policy=diploid_x" if opts.ploidy_policy == "sex-aware" else "") ) lines.append(f"##simurg_run={rm.run_id}") for c in contigs: lines.append(f"##contig=") lines.append('##FILTER=') + for key in opts.info_fields: + kind, desc = INFO_DEFS[key] + lines.append(f'##INFO=') if opts.filter_gs is not None: lines.append(f'##FILTER=') lines += [ @@ -289,6 +381,9 @@ def vcf_header(opts: VcfOptions, genome: Any, build: str, contigs: list[str], rm 'ALT (REF=A, ALT=B relative to the reference genome)">', '##FORMAT=', ] + for key in opts.extra_format: + kind, desc = EXTRA_DEFS[key] + lines.append(f'##FORMAT=') return "\n".join(lines) + "\n" @@ -338,6 +433,22 @@ def write_sample(task: tuple[str, str, str, str, bool]) -> dict[str, Any]: b = np.where(t.single_flip, (np.float32(1.0) - b).astype(np.float32), b) baf[rows] = b.astype(np.float64) lrr[rows] = g.logr_ratios[loc].astype(np.float64) + extra_cols: list[list[str]] = [] + if t.extra_format: + from simurg.signal.normalization.apply import normalize, polar + + rx, ry = g.raw_x[loc], g.raw_y[loc] + nx, ny = normalize(rx, ry, t.single_lookup, g.normalization_transforms) # type: ignore[arg-type] + zero = (rx == 0) & (ry == 0) # no intensity: normalized values are missing (as in the GTC semantics) + nx = np.where(zero, np.float32(np.nan), nx).astype(np.float32) + ny = np.where(zero, np.float32(np.nan), ny).astype(np.float32) + th, rr = polar(nx, ny) + src: dict[str, npt.NDArray[Any]] = {"X": rx, "Y": ry, "NORMX": nx, "NORMY": ny, "THETA": th, "R": rr} + for key in t.extra_format: + col = np.full(n, "", dtype=object) + vals = src[key] + col[rows] = [str(int(v)) for v in vals] if key in ("X", "Y") else [_g8(float(v)) for v in vals] + extra_cols.append(col.tolist()) het_haploid = int(((calls == 2) & (pl == 1)).sum()) if het_haploid: res["reasons"]["het_in_haploid"] += het_haploid @@ -373,10 +484,16 @@ def write_sample(task: tuple[str, str, str, str, bool]) -> dict[str, Any]: raw.write(header.encode()) # type: ignore[arg-type] for k in range(n): filt = "LowGS" if low is not None and not math.isnan(gsl[k]) and gsl[k] < low else "PASS" + info = t.site_info[k] if t.site_info else "." if merged_mask[k]: - body = f"{filt}\t.\tGT:GS\t{gt[k]}:{_f4(gsl[k])}\n" + body = f"{filt}\t{info}\tGT:GS\t{gt[k]}:{_f4(gsl[k])}\n" else: - body = f"{filt}\t.\tGT:GS:BAF:LRR\t{gt[k]}:{_f4(gsl[k])}:{_g8(bafl[k])}:{_g8(lrrl[k])}\n" + extra_keys = "".join(":" + key for key in t.extra_format) + extra_vals = "".join(":" + c[k] for c in extra_cols) + body = ( + f"{filt}\t{info}\tGT:GS:BAF:LRR{extra_keys}\t" + f"{gt[k]}:{_f4(gsl[k])}:{_g8(bafl[k])}:{_g8(lrrl[k])}{extra_vals}\n" + ) line = t.prefix[k] + body.encode() if bgzip: vs, ve = w.write_record(line) @@ -413,10 +530,20 @@ def run_gtc_to_vcf(opts: VcfOptions, rep: Reporter) -> int: (opts.sheet, "sample sheet"), (opts.exclude_loci, "locus list"), (opts.include_loci, "locus list"), + (opts.auxiliary_loci, "auxiliary loci VCF"), ] ) if opts.ploidy_policy not in ("sex-aware", "diploid"): raise SimurgError("E_CONFIG_KEY", "--ploidy-policy must be sex-aware or diploid") + bad = [x for x in opts.info_fields if x not in INFO_DEFS] + [x for x in opts.extra_format if x not in EXTRA_DEFS] + if bad: + raise SimurgError( + "E_CONFIG_KEY", + f"unknown field(s) {', '.join(bad)}; --info-fields takes {','.join(INFO_DEFS)}; " + f"--extra-format takes {','.join(EXTRA_DEFS)}", + ) + if {"GENTRAIN", "CLUSTER_SEP"} & set(opts.info_fields) and opts.egt is None: + raise SimurgError("E_CONFIG_KEY", "--info-fields GENTRAIN/CLUSTER_SEP need the cluster file (--egt)") sheet = parse_samplesheet(opts.sheet) if opts.sheet else None root = opts.gtcs or (sheet.root_folder if sheet else None) if root is None: @@ -534,25 +661,31 @@ def write_cohort_bcf( from simurg.formats.bcf import ( FLOAT, INT8, + INT32, BcfWriter, CsiBuilder, _type_byte, encode_record, header_dictionaries, + typed_floats, + typed_ints, + typed_string, ) from simurg.formats.bpm import read_bpm from simurg.pipeline.qc import default_tmp, load_cohort + from simurg.signal.normalization.apply import polar - n_loci = read_bpm(opts.bpm).num_loci + bpm = read_bpm(opts.bpm) workers = max(1, opts.threads or (os.cpu_count() or 1)) store, infos = load_cohort( [(sid, path) for sid, _, path in samples], - n_loci, + bpm.num_loci, default_tmp(out.parent), baf=True, scores=True, lrr=True, workers=workers, + lookup=bpm.normalization_lookups() if t.extra_format else None, ) try: for s in infos: @@ -591,6 +724,9 @@ def write_cohort_bcf( sc = np.asarray(store.score[:, sel]) if len(sel) else np.zeros((n, 0), np.float32) # type: ignore[index] bf = np.asarray(store.baf[:, sel]) if len(sel) else np.zeros((n, 0), np.float32) # type: ignore[index] lr = np.asarray(store.lrr[:, sel]) if len(sel) else np.zeros((n, 0), np.float32) # type: ignore[index] + if t.extra_format and len(sel): + ex_rx, ex_ry = np.asarray(store.rx[:, sel]), np.asarray(store.ry[:, sel]) # type: ignore[index] + ex_nx, ex_ny = np.asarray(store.nx[:, sel]), np.asarray(store.ny[:, sel]) # type: ignore[index] for k in rows: pl = ploidy[k] if single_pos[k] >= 0: @@ -616,6 +752,25 @@ def write_cohort_bcf( b = np.where(t.single_flip[j], (np.float32(1.0) - b).astype(np.float32), b) fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs)), (k_baf, fdesc, _f32_payload(b)), (k_lrr, fdesc, _f32_payload(lr[:, c]))] # fmt: skip + if t.extra_format: + zero = (ex_rx[:, c] == 0) & (ex_ry[:, c] == 0) + enx = np.where(zero, np.float32(np.nan), ex_nx[:, c]).astype(np.float32) + eny = np.where(zero, np.float32(np.nan), ex_ny[:, c]).astype(np.float32) + th, rr = polar(enx, eny) + xsrc = { + "X": ex_rx[:, c], + "Y": ex_ry[:, c], + "NORMX": enx, + "NORMY": eny, + "THETA": th, + "R": rr, + } + for key in t.extra_format: + if key in ("X", "Y"): + payload = xsrc[key].astype("= 0 else f_pass, n, fmt) # fmt: skip + f_low if low and f_low >= 0 else f_pass, n, fmt, info) # fmt: skip vs, ve = w.write(rec) csi.add(contigs[t.site_contig[k]], int(t.site_beg[k]), int(t.site_end[k]), vs, ve) written += 1 diff --git a/tests/oracle/test_p6_vcf.py b/tests/oracle/test_p6_vcf.py index 3acce8c..211cd83 100644 --- a/tests/oracle/test_p6_vcf.py +++ b/tests/oracle/test_p6_vcf.py @@ -148,3 +148,74 @@ def test_bedgraph_equals_gtc(tmp_path: Path) -> None: # OUT-050, OUT-052, EDGE- assert v == f"{float(g.logr_ratios[single[key]]):.6g}" checked += 1 assert checked > 500_000 + + +def test_info_and_extra_format_match_gtc2vcf(tmp_path: Path) -> None: # OUT-012, plan §8.3 + """INFO ASSAY_TYPE/GENTRAIN/CLUSTER_SEP and FORMAT X/Y/NORMX/NORMY equal gtc2vcf's; THETA/R too, except where the + normalized intensities are both 0 (SimurgArray: missing, plan §7.4; gtc2vcf: 0).""" + from tests.conftest import GSA_EGT + + plugin = BCFTOOLS.parent / "plugins" / "gtc2vcf.so" + need(GSA_BPM, GSA_CSV, GSA_EGT, FASTA, ORACLE_GTCS) + if not plugin.exists(): + pytest.skip("gtc2vcf plugin not built") + import os + + from typer.testing import CliRunner + + from simurg.cli.main import app + + gtcs = tmp_path / "g" + gtcs.mkdir() + src = sorted(ORACLE_GTCS.glob("*.gtc"))[0] + (gtcs / src.name).symlink_to(src) + out = tmp_path / "v" + res = CliRunner().invoke( + app, + ["gtc-to-vcf", "--bpm", str(GSA_BPM), "--csv", str(GSA_CSV), "--fasta", str(FASTA), "--gtcs", str(gtcs), + "--out", str(out), "--threads", "1", "--unsquash-duplicates", "--quiet", "--egt", str(GSA_EGT), + "--info-fields", "ASSAY_TYPE,GENTRAIN,CLUSTER_SEP", "--extra-format", "X,Y,NORMX,NORMY,THETA,R"], + ) # fmt: skip + assert res.exit_code == 0, res.output + env = dict(os.environ, BCFTOOLS_PLUGINS=str(plugin.parent)) + g2v = subprocess.run( + [str(BCFTOOLS), "+gtc2vcf", "--no-version", "-t", "GT,X,Y,NORMX,NORMY,THETA,R", "-b", str(GSA_BPM), "-c", + str(GSA_CSV), "-e", str(GSA_EGT), "-f", str(FASTA), "-g", str(gtcs), "-O", "v"], + capture_output=True, text=True, env=env, check=True, + ).stdout # fmt: skip + ref = {} + for line in g2v.splitlines(): + if line.startswith("#"): + continue + c = line.split("\t") + info = dict(kv.split("=", 1) for kv in c[7].split(";") if "=" in kv) + ref[c[2]] = (info, dict(zip(c[8].split(":"), c[9].split(":"), strict=True))) + + def close(a: str, b: str) -> bool: + na, nb = a in (".", "") or "nan" in a.lower(), b in (".", "") or "nan" in b.lower() + if na or nb: + return na and nb + return abs(float(a) - float(b)) <= 1e-5 * max(1.0, abs(float(b))) + + checked = zero_norm = 0 + with gzip.open(out / f"{src.stem}.snv.vcf.gz", "rt") as fh: + for line in fh: + if line.startswith("#"): + continue + c = line.rstrip("\n").split("\t") + if "," in c[2] or ";" in c[2] or c[2] not in ref: + continue + info = dict(kv.split("=", 1) for kv in c[7].split(";")) + fmt = dict(zip(c[8].split(":"), c[9].split(":"), strict=True)) + ri, rf = ref[c[2]] + assert info["ASSAY_TYPE"] == ri["ASSAY_TYPE"] + assert close(info["GENTRAIN"], ri["GenTrain_Score"]) and close(info["CLUSTER_SEP"], ri["Cluster_Sep"]) + for k in ("X", "Y", "NORMX", "NORMY"): + assert close(fmt[k], rf[k]), (c[2], k, fmt[k], rf[k]) + if fmt["NORMX"] == "0" and fmt["NORMY"] == "0": + zero_norm += 1 + assert fmt["THETA"] == "." and float(rf["THETA"]) == 0.0 + else: + assert close(fmt["THETA"], rf["THETA"]) and close(fmt["R"], rf["R"]), (c[2], fmt, rf) + checked += 1 + assert checked > 600_000 and zero_norm < 1000 diff --git a/tests/unit/test_vcf.py b/tests/unit/test_vcf.py index 1951ec6..395e115 100644 --- a/tests/unit/test_vcf.py +++ b/tests/unit/test_vcf.py @@ -232,3 +232,64 @@ def test_ploidy_rules() -> None: # OUT-018, EDGE-038 def test_plink_codes() -> None: # OUT-043 assert [bed_code(g, 1) for g in ("0/0", "0/1", "1/1", "./.", "1", "0")] == [3, 2, 0, 1, 0, 3] assert bed_code("0/2", 1) == 3 and bed_code("0/2", 2) == 2 # split multi-allelic record + + +# ------------------------------------------------------------------------------------------------ P6 options +def test_merge_snv_indel_edge_013() -> None: # OUT-021 + snv = ProbeSite(0, "snv", "1", 5, "A", ("A", "G"), "snv", 2, "1") + dele = ProbeSite(1, "del", "1", 5, "ACT", ("ACT", "A"), "indel", 1, "1") + sep = group_sites([snv, dele], {"1": 0}) + assert len(sep) == 2 # default: separate records + (s,) = group_sites([snv, dele], {"1": 0}, merge_snv_indel=True) + assert s.ref == "ACT" and s.alts == ["GCT", "A"] and s.id == "snv;del" + ext_snv, ext_del = s.probes + assert ext_del.alleles == ("ACT", "A") + assert ext_snv.alleles == ("ACT", "GCT") + # SNV probe homozygous A (sees the anchor base of the deletion too), deletion probe het → REF/DEL + assert gt_text(genotype(s, (AA, AB), 2), 2) == "0/2" + # SNV het A/G, deletion hom REF → REF/SNV + assert gt_text(genotype(s, (AB, AA), 2), 2) == "0/1" + # incompatible REFs (not prefixes of each other) stay separate + other = ProbeSite(2, "x", "1", 5, "GT", ("GT", "G"), "indel", 1, "1") + assert len(group_sites([snv, other], {"1": 0}, merge_snv_indel=True)) == 2 + + +def _aux_probe(alleles: tuple[str, str] = ("A", "C"), idx: int = 0, name: str = "mnv") -> ProbeSite: + return ProbeSite(idx, name, "1", 5, "A", alleles, "snv", 2, "1") + + +def test_auxiliary_loci_rules() -> None: # OUT-035, EDGE-015 + from simurg.harmonize.sites import AuxRecord, apply_auxiliary + + genome = "TTTTATGCC" # positions 5-6 = AT + f = fetcher(genome) + plus = lambda i: "+" # noqa: E731 + ok = AuxRecord("mnv", "1", 5, "AT", ["CG"]) + probes, errs = apply_auxiliary([_aux_probe()], [ok], f, lambda c: c, plus) + assert not errs and probes[0].ref == "AT" and probes[0].alleles == ("AT", "CG") and probes[0].reason == "aux_locus" + cases = { + "(a)": AuxRecord("mnv", "1", 5, "AT", ["CG", "GG"]), + "(b)": AuxRecord("mnv", "1", 5, "AT", ["C"]), + "(d)": AuxRecord("mnv", "1", 5, "AT", ["TG"]), + "(e)": AuxRecord("mnv", "1", 5, "AG", ["CG"]), + } + for rule, rec in cases.items(): + _, errs = apply_auxiliary([_aux_probe()], [rec], f, lambda c: c, plus) + assert len(errs) == 1 and rule in errs[0], (rule, errs) + dup = [_aux_probe(), _aux_probe(idx=1, name="twin")] + _, errs = apply_auxiliary(dup, [ok], f, lambda c: c, plus) + assert "(c)" in errs[0] + _, errs = apply_auxiliary([_aux_probe()], [AuxRecord("nope", "1", 5, "AT", ["CG"])], f, lambda c: c, plus) + assert "not a probe" in errs[0] + # minus-strand design: the interrogated base is the MNV's last base on the plus strand + probes, errs = apply_auxiliary([_aux_probe(("T", "G"))], [AuxRecord("mnv", "1", 5, "AT", ["CG"])], f, + lambda c: c, lambda i: "-") # fmt: skip + assert not errs and probes[0].alleles == ("AT", "CG") + + +def test_pgen_codes() -> None: # plan §8.4 + from simurg.pipeline.export import PGEN_FROM_BED, _pack + + # PLINK 1 (A1 = ALT): 0 hom ALT, 1 missing, 2 het, 3 hom REF → PLINK 2: 2, 3, 1, 0 + assert list(PGEN_FROM_BED) == [2, 3, 1, 0] + assert _pack([0, 1, 2, 3, 0], 2) == bytes([0b11100100, 0b00000000])