P6 complete: optional VCF fields, SNV+indel merging, auxiliary MNV loci, PLINK 2 .pgen
- 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
This commit is contained in:
kemal committed
1 parent
5b85e4e35f
commit
4d190d559e
13 files changed
+611
-50
No files matched your search
@@ -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("<II", 0, n) + b"\x40") # variant count patched below
|
||||
pvar.write(pvar_header + "#CHROM\tPOS\tID\tREF\tALT\n")
|
||||
for recs in zip(*streams, strict=True):
|
||||
r0 = recs[0]
|
||||
if any(r[:5] != r0[:5] for r in recs[1:]):
|
||||
@@ -405,18 +434,49 @@ def run_export_plink(opts: PlinkOptions, rep: Reporter) -> 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("<I", n_var))
|
||||
sex_code = {"M": 1, "F": 2}
|
||||
if want_bed:
|
||||
fam = "".join(f"{s}\t{s}\t0\t0\t{sex_code.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:
|
||||
psam = "#IID\tSEX\n" + "".join(f"{s}\t{ {'M': '1', 'F': '2'}.get(sex.get(s, 'U'), 'NA') }\n" for s in names)
|
||||
with atomic_open(base.with_suffix(".psam"), "w", "utf-8") as psam_raw:
|
||||
cast(IO[str], psam_raw).write(psam)
|
||||
rm.status = "completed"
|
||||
rm.summary = {"samples": n, "variants": n_var, "multiallelic_dropped": dropped}
|
||||
rm.summary = {"samples": n, "variants": n_var, "multiallelic_dropped": dropped, "format": opts.format}
|
||||
rm.write(opts.out / "run_manifest.export_plink.json")
|
||||
rep.info(f"PLINK fileset {base}.bed/.bim/.fam: {n} samples, {n_var} variants")
|
||||
rep.info(f"PLINK fileset {base} ({opts.format}): {n} samples, {n_var} variants")
|
||||
return ExitCode.OK
|
||||
|
||||
|
||||
class _Null:
|
||||
"""Context manager standing in for an output that is not requested."""
|
||||
|
||||
def __enter__(self) -> 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
|
||||
Reference in new issue
Block a user