kemalandClaude Opus 5.5
2e8478cc6a
- 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 Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
821 lines
37 KiB
Python
821 lines
37 KiB
Python
"""``simurg gtc-to-vcf``: GTC + manifest + FASTA → one bgzipped, indexed SNV VCF per sample (spec CLI-022,
|
||
OUT-010 … OUT-042, IN-060 … IN-064; plan §8).
|
||
|
||
The manifest is mapped to plus-strand sites once (``simurg.harmonize``); each sample is then written by a worker
|
||
process with vectorised genotype lookup for single-probe sites and the OUT-031 algorithm for merged sites.
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import hashlib
|
||
import math
|
||
import multiprocessing as mp
|
||
import os
|
||
from collections import Counter
|
||
from concurrent.futures import ProcessPoolExecutor
|
||
from dataclasses import dataclass, field
|
||
from pathlib import Path
|
||
from typing import Any
|
||
|
||
import numpy as np
|
||
import numpy.typing as npt
|
||
|
||
from simurg import PRODUCT_NAME, __version__
|
||
from simurg.core.atomic import atomic_open, check_outputs_absent
|
||
from simurg.core.errors import ExitCode, SimurgError
|
||
from simurg.core.log import Reporter
|
||
from simurg.core.provenance import RunManifest, sha256_file
|
||
from simurg.formats.bgzf import BgzfWriter, TabixBuilder, write_tbi
|
||
from simurg.harmonize.combine import genotype as combine_genotype
|
||
from simurg.harmonize.combine import gt_text
|
||
from simurg.harmonize.sites import AuditRow, ManifestView, Site, group_sites, map_locus
|
||
|
||
PAR = { # 1-based inclusive pseudo-autosomal intervals on X (IN-065, human)
|
||
"GRCh38": ((10_001, 2_781_479), (155_701_383, 156_030_895)),
|
||
"GRCh37": ((60_001, 2_699_520), (154_931_044, 155_260_560)),
|
||
}
|
||
CLASS_AUTO, CLASS_X, CLASS_Y, CLASS_MT, CLASS_PAR = 0, 1, 2, 3, 4
|
||
|
||
|
||
@dataclass
|
||
class VcfOptions:
|
||
out: Path
|
||
bpm: Path
|
||
fasta: Path
|
||
gtcs: Path | None = None
|
||
sheet: Path | None = None
|
||
csv: Path | None = None
|
||
bgzip: bool = True
|
||
genome_cache: bool = True
|
||
exclude_loci: Path | None = None
|
||
include_loci: Path | None = None
|
||
unsquash_duplicates: bool = False
|
||
inf1_nocall_informative: bool = False
|
||
drop_ambiguous: bool = False
|
||
filter_gs: float | None = None
|
||
contig_alias: str = "strict"
|
||
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)
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ site table
|
||
@dataclass
|
||
class SiteTable:
|
||
"""Everything a worker needs, precomputed once (picklable)."""
|
||
|
||
prefix: list[bytes] # "CHROM\tPOS\tID\tREF\tALT\t.\t" per site
|
||
single_rows: npt.NDArray[np.int64] # site indices of single-probe sites
|
||
single_locus: npt.NDArray[np.int64] # their manifest locus index
|
||
single_gt: npt.NDArray[np.object_] # (n_single, 2 ploidy, 4 calls) GT text
|
||
single_flip: npt.NDArray[np.bool_] # BAF must be 1 − BAF
|
||
merged_sites: list[Site]
|
||
merged_rows: npt.NDArray[np.int64]
|
||
site_class: npt.NDArray[np.int8]
|
||
contigs: list[str]
|
||
header: str
|
||
site_beg: npt.NDArray[np.int64]
|
||
site_end: npt.NDArray[np.int64]
|
||
site_contig: list[str]
|
||
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"):
|
||
return CLASS_MT
|
||
if c == "Y":
|
||
return CLASS_Y
|
||
if c == "X":
|
||
if par_flag or any(a <= pos <= b for a, b in PAR.get(build or "GRCh38", ())):
|
||
return CLASS_PAR
|
||
return CLASS_X
|
||
return CLASS_AUTO
|
||
|
||
|
||
def ploidy_for(cls: npt.NDArray[np.int8], sex: str, policy: str) -> npt.NDArray[np.int8]:
|
||
"""OUT-018 per site; 0 marks 'haploid and always missing' (Y of female samples)."""
|
||
p = np.full(len(cls), 2, dtype=np.int8)
|
||
p[cls == CLASS_MT] = 1
|
||
if policy == "diploid":
|
||
return p
|
||
p[cls == CLASS_Y] = 0 if sex == "F" else 1
|
||
if sex == "M":
|
||
p[cls == CLASS_X] = 1
|
||
return p
|
||
|
||
|
||
def read_locus_list(path: Path) -> list[str]:
|
||
out = []
|
||
for line in path.read_text(encoding="utf-8").splitlines():
|
||
line = line.split("#", 1)[0].strip()
|
||
if line:
|
||
out.append(line)
|
||
return out
|
||
|
||
|
||
def build_sites(
|
||
opts: VcfOptions, rep: Reporter, rm: RunManifest
|
||
) -> tuple[SiteTable, list[AuditRow], Counter[str], Any]:
|
||
from simurg.formats.bpm import read_bpm
|
||
from simurg.formats.fasta import ContigMap, Genome
|
||
|
||
bpm = read_bpm(opts.bpm)
|
||
csvm = None
|
||
if opts.csv is not None:
|
||
from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest
|
||
|
||
csvm = read_csv_manifest(opts.csv)
|
||
problems = compare_manifests(bpm, csvm)
|
||
if problems:
|
||
raise SimurgError("E_BPM_CSV_MISMATCH", "BPM and CSV manifests differ:\n " + "\n ".join(problems))
|
||
genome = Genome(opts.fasta, use_cache=opts.genome_cache)
|
||
for n in genome.notes:
|
||
rep.info(n)
|
||
builds = Counter(bpm.genome_builds)
|
||
mbuild = builds.most_common(1)[0][0] if builds else ""
|
||
want = {"37": "GRCh37", "38": "GRCh38"}.get(mbuild.strip(), mbuild.strip())
|
||
if genome.build and want and want != genome.build:
|
||
raise SimurgError(
|
||
"E_BUILD_MISMATCH", f"manifest genome build {want} but FASTA looks like {genome.build} ({opts.fasta})"
|
||
)
|
||
build = genome.build or want
|
||
src = csvm or bpm
|
||
mv = ManifestView(
|
||
bpm.names,
|
||
bpm.snps,
|
||
bpm.chroms,
|
||
bpm.map_infos.tolist(),
|
||
src.ilmn_strands,
|
||
src.source_strands,
|
||
src.ref_strands if any(src.ref_strands) else bpm.ref_strands,
|
||
csvm.source_seqs if csvm else [""] * bpm.num_loci,
|
||
bpm.assay_types.tolist(),
|
||
bpm.intensity_only.tolist(),
|
||
)
|
||
if csvm is None and any(("D" in s) for s in bpm.snps):
|
||
rep.error("W_SKIP_INDEL_CONTEXT", "no CSV manifest: indels are skipped (give --csv with SourceSeq)")
|
||
names = set(bpm.names)
|
||
excl: set[str] = set()
|
||
incl: set[str] | None = None
|
||
for path, kind in ((opts.exclude_loci, "exclude"), (opts.include_loci, "include")):
|
||
if path is None:
|
||
continue
|
||
lst = read_locus_list(path)
|
||
unknown = [x for x in lst if x not in names]
|
||
if unknown:
|
||
rep.warn("W_UNKNOWN_LOCI", f"{len(unknown)} name(s) in --{kind}-loci are not in the manifest")
|
||
if kind == "exclude":
|
||
excl = set(lst)
|
||
else:
|
||
incl = set(lst)
|
||
cmap = ContigMap([e.name for e in genome.entries], opts.contig_alias)
|
||
probes = []
|
||
audit: list[AuditRow] = []
|
||
skips: Counter[str] = Counter()
|
||
for i in range(bpm.num_loci):
|
||
nm = bpm.names[i]
|
||
if nm in excl or (incl is not None and nm not in incl):
|
||
a = AuditRow(nm, bpm.snps[i], "", "", "", "", "", False, False, "W_SKIP_EXCLUDED", build, "")
|
||
p = None
|
||
else:
|
||
p, a = map_locus(i, mv, genome.fetch, cmap.get, build)
|
||
if p is not None and opts.drop_ambiguous and p.confidence == "low":
|
||
a.reason, p = "W_SKIP_AMBIGUOUS", None
|
||
audit.append(a)
|
||
if p is None or a.reason.startswith("W_"):
|
||
skips[a.reason] += 1
|
||
if p is not None:
|
||
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, 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
|
||
prefix, cls = [], np.zeros(len(sites), np.int8)
|
||
single_rows, single_locus, single_flip, gts, single_ab = [], [], [], [], []
|
||
merged_sites, merged_rows = [], []
|
||
beg = np.zeros(len(sites), np.int64)
|
||
end = np.zeros(len(sites), np.int64)
|
||
for k, s in enumerate(sites):
|
||
prefix.append(f"{s.contig}\t{s.pos}\t{s.id}\t{s.ref}\t{','.join(s.alts) if s.alts else '.'}\t.\t".encode())
|
||
cls[k] = _site_class(s.contig, s.pos, any(p.par for p in s.probes), build)
|
||
beg[k], end[k] = s.pos - 1, s.pos - 1 + len(s.ref)
|
||
if len(s.probes) == 1:
|
||
from simurg.harmonize.combine import single
|
||
|
||
single_rows.append(k)
|
||
single_locus.append(s.probes[0].index)
|
||
al_a, al_b = s.probes[0].alleles
|
||
single_flip.append(al_a != s.ref and al_b == s.ref)
|
||
single_ab.append((s.allele_index(al_a), s.allele_index(al_b)))
|
||
gts.append([[gt_text(single(s, c, pl), pl) for c in range(4)] for pl in (1, 2)])
|
||
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,
|
||
single_rows=np.array(single_rows, np.int64),
|
||
single_locus=np.array(single_locus, np.int64),
|
||
single_gt=np.array(gts, dtype=object).reshape(len(single_rows), 2, 4),
|
||
single_flip=np.array(single_flip, np.bool_),
|
||
merged_sites=merged_sites,
|
||
merged_rows=np.array(merged_rows, np.int64),
|
||
site_class=cls,
|
||
contigs=contigs,
|
||
header=header,
|
||
site_beg=beg,
|
||
site_end=end,
|
||
site_contig=[s.contig for s in 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,
|
||
)
|
||
return table, audit, skips, genome
|
||
|
||
|
||
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)."""
|
||
|
||
def lst(p: Path | None) -> str:
|
||
return "none" if p is None else f"{p.name}:{hashlib.sha256(p.read_bytes()).hexdigest()[:16]}"
|
||
|
||
lines = [
|
||
"##fileformat=VCFv4.3",
|
||
f"##source={PRODUCT_NAME} {__version__}",
|
||
f"##reference=file://{Path(opts.fasta).resolve()}",
|
||
f"##simurg_genome_build={build}",
|
||
f"##simurg_manifest={opts.bpm.name};sha256={_sha(opts.bpm)}",
|
||
]
|
||
if opts.csv:
|
||
lines.append(f"##simurg_csv_manifest={opts.csv.name};sha256={_sha(opts.csv)}")
|
||
lines.append("##simurg_cluster_file=__CLUSTER__")
|
||
lines.append(
|
||
"##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={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=<ID={c},length={genome.contig_length(c)}>")
|
||
lines.append('##FILTER=<ID=PASS,Description="All filters passed">')
|
||
for key in opts.info_fields:
|
||
kind, desc = INFO_DEFS[key]
|
||
lines.append(f'##INFO=<ID={key},Number=.,Type={kind},Description="{desc}">')
|
||
if opts.filter_gs is not None:
|
||
lines.append(f'##FILTER=<ID=LowGS,Description="GenCall score below {opts.filter_gs}">')
|
||
lines += [
|
||
'##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">',
|
||
'##FORMAT=<ID=GS,Number=1,Type=Float,Description="GenCall score. For records merged from several assays, '
|
||
'the minimum GenCall score over the contributing assays">',
|
||
'##FORMAT=<ID=BAF,Number=1,Type=Float,Description="B allele frequency, oriented so that 1 means homozygous '
|
||
'ALT (REF=A, ALT=B relative to the reference genome)">',
|
||
'##FORMAT=<ID=LRR,Number=1,Type=Float,Description="Log R ratio">',
|
||
]
|
||
for key in opts.extra_format:
|
||
kind, desc = EXTRA_DEFS[key]
|
||
lines.append(f'##FORMAT=<ID={key},Number=1,Type={kind},Description="{desc}">')
|
||
return "\n".join(lines) + "\n"
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ per sample
|
||
_W: dict[str, Any] = {}
|
||
|
||
|
||
def _init(table: SiteTable) -> None:
|
||
_W["t"] = table
|
||
|
||
|
||
def _f4(v: float) -> str:
|
||
return "." if math.isnan(v) else f"{v:.4f}"
|
||
|
||
|
||
def _g8(v: float) -> str:
|
||
return "." if math.isnan(v) else f"{v:.8g}"
|
||
|
||
|
||
def write_sample(task: tuple[str, str, str, str, bool]) -> dict[str, Any]:
|
||
"""Worker: one GTC → one VCF (+ index). Returns counts for the run summary."""
|
||
from simurg.formats.gtc import read_gtc
|
||
|
||
sample_id, sample_name, gtc_path, out_path, bgzip = task
|
||
t: SiteTable = _W["t"]
|
||
res: dict[str, Any] = {"sample_id": sample_id, "records": 0, "reasons": Counter(), "merge_failures": []}
|
||
try:
|
||
g = read_gtc(gtc_path)
|
||
n = len(t.prefix)
|
||
sex = g.gender_char
|
||
ploidy = ploidy_for(t.site_class, sex, t.ploidy_policy)
|
||
gt = np.empty(n, dtype=object)
|
||
gs = np.full(n, np.nan, np.float64)
|
||
baf = np.full(n, np.nan, np.float64)
|
||
lrr = np.full(n, np.nan, np.float64)
|
||
has_baf = len(g.b_allele_freqs) == g.num_snps
|
||
# single-probe sites (vectorised)
|
||
rows, loc = t.single_rows, t.single_locus
|
||
calls = g.genotypes[loc].astype(np.int64)
|
||
pl = ploidy[rows]
|
||
pidx = np.where(pl == 2, 1, 0)
|
||
gt[rows] = t.single_gt[np.arange(len(rows)), pidx, calls]
|
||
gt[rows[pl == 0]] = "."
|
||
gs[rows] = g.genotype_scores[loc].astype(np.float64)
|
||
if has_baf:
|
||
b = g.b_allele_freqs[loc].astype(np.float32)
|
||
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
|
||
# merged sites
|
||
for s, k in zip(t.merged_sites, t.merged_rows.tolist(), strict=True):
|
||
c = tuple(int(g.genotypes[p.index]) for p in s.probes)
|
||
p_k = int(ploidy[k])
|
||
gs[k] = float(min(float(g.genotype_scores[p.index]) for p in s.probes))
|
||
if p_k == 0:
|
||
gt[k] = "."
|
||
continue
|
||
o = combine_genotype(s, c, p_k, t.inf1_nc_informative)
|
||
gt[k] = gt_text(o, p_k)
|
||
if o.reason:
|
||
res["reasons"][o.reason] += 1
|
||
if o.reason.startswith("W_MERGE") and o.reason != "W_MERGE_ALL_NC":
|
||
res["merge_failures"].append((s.id, c, o.reason))
|
||
merged_mask = np.zeros(n, np.bool_)
|
||
merged_mask[t.merged_rows] = True
|
||
header = t.header.replace("__CLUSTER__", g.cluster_file or "") + (
|
||
f"##SAMPLE=<ID={sample_name},SentrixID={sample_id}>\n"
|
||
f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}\n"
|
||
)
|
||
low = t.filter_gs
|
||
gsl, bafl, lrrl = gs.tolist(), baf.tolist(), lrr.tolist()
|
||
tbx = TabixBuilder()
|
||
out = Path(out_path)
|
||
with atomic_open(out, "wb") as raw:
|
||
if bgzip:
|
||
w = BgzfWriter(raw) # type: ignore[arg-type]
|
||
w.write_record(header.encode())
|
||
else:
|
||
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{info}\tGT:GS\t{gt[k]}:{_f4(gsl[k])}\n"
|
||
else:
|
||
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)
|
||
tbx.add(t.site_contig[k], int(t.site_beg[k]), int(t.site_end[k]), vs, ve)
|
||
else:
|
||
raw.write(line) # type: ignore[arg-type]
|
||
if bgzip:
|
||
w.close()
|
||
if bgzip:
|
||
write_tbi(Path(str(out) + ".tbi"), tbx)
|
||
res["records"] = n
|
||
res["sex"] = sex
|
||
except SimurgError as e:
|
||
res["error"] = f"{e.code}: {e.message}"
|
||
except Exception as e:
|
||
res["error"] = f"E_INTERNAL: {type(e).__name__}: {e}"
|
||
return res
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ command
|
||
def run_gtc_to_vcf(opts: VcfOptions, rep: Reporter) -> int:
|
||
from simurg.formats.discovery import discover_gtcs
|
||
from simurg.formats.samplesheet import output_stem, parse_samplesheet
|
||
from simurg.pipeline.inputs import check_paths
|
||
|
||
rm = RunManifest(command="gtc-to-vcf", argv=opts.argv)
|
||
rep.run_id = rm.run_id
|
||
check_paths(
|
||
[
|
||
(opts.bpm, "BPM"),
|
||
(opts.fasta, "FASTA"),
|
||
(opts.csv, "CSV manifest"),
|
||
(opts.gtcs, "GTC folder"),
|
||
(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:
|
||
raise SimurgError("E_NO_INPUT", "give --gtcs or a sample sheet with a RootFolder")
|
||
found = discover_gtcs(root)
|
||
samples: list[tuple[str, str, str]] = [] # (sample_id, sample_name, path)
|
||
if sheet is not None:
|
||
for r in sheet.rows:
|
||
if r.sample_id in found:
|
||
samples.append((r.sample_id, r.sample_name or r.sample_id, str(found[r.sample_id])))
|
||
else:
|
||
rep.warn("W_QC_UNMATCHED_SAMPLES", f"{r.sample_id}: no GTC under {root}", sample_id=r.sample_id)
|
||
else:
|
||
samples = [(k, k, str(v)) for k, v in sorted(found.items())]
|
||
suffix = ".snv.vcf.gz" if opts.bgzip else ".snv.vcf"
|
||
outs = (
|
||
[opts.out / "cohort.snv.bcf"]
|
||
if opts.merge
|
||
else [opts.out / f"{output_stem(name)}{suffix}" for _, name, _ in samples]
|
||
)
|
||
check_outputs_absent([*outs, opts.out / "audit" / "loci.parquet"], force=opts.force)
|
||
if not opts.bgzip:
|
||
rep.warn("W_NO_BGZIP", "writing plain VCF without index; pgx and cn commands need bgzipped, indexed VCFs")
|
||
rep.info(f"mapping manifest {opts.bpm.name} to {opts.fasta.name}")
|
||
table, audit, skips, _genome = build_sites(opts, rep, rm)
|
||
for role, p in (("bpm", opts.bpm), ("csv_manifest", opts.csv), ("fasta", opts.fasta), ("sample_sheet", opts.sheet)):
|
||
if p is not None:
|
||
rm.add_input(role, p)
|
||
opts.out.mkdir(parents=True, exist_ok=True)
|
||
write_audit(opts.out / "audit" / "loci.parquet", audit)
|
||
if opts.merge:
|
||
n = write_cohort_bcf(table, samples, opts.out / "cohort.snv.bcf", opts, rep)
|
||
rm.summary = {"samples": len(samples), "records": n, "skipped_loci_by_code": dict(skips)}
|
||
rm.add_output(opts.out / "cohort.snv.bcf", "bcf/2.2")
|
||
rm.status = "completed"
|
||
rm.write(opts.out / "run_manifest.json")
|
||
rep.info(f"done: cohort BCF with {len(samples)} sample(s), {n} records")
|
||
return ExitCode.OK
|
||
workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(samples)))
|
||
tasks = [(sid, name, path, str(o), opts.bgzip) for (sid, name, path), o in zip(samples, outs, strict=True)]
|
||
rep.info(f"{len(table.prefix)} VCF records; writing {len(tasks)} sample(s) with {workers} worker(s)")
|
||
if workers == 1:
|
||
_init(table)
|
||
results = [write_sample(t) for t in tasks]
|
||
else:
|
||
with ProcessPoolExecutor(
|
||
workers, mp_context=mp.get_context("spawn"), initializer=_init, initargs=(table,)
|
||
) as ex:
|
||
results = list(ex.map(write_sample, tasks))
|
||
reasons: Counter[str] = Counter()
|
||
failures = []
|
||
for res in results:
|
||
if "error" in res:
|
||
rep.error(res["error"].split(":", 1)[0], f"{res['sample_id']}: {res['error']}", sample_id=res["sample_id"])
|
||
continue
|
||
reasons.update(res["reasons"])
|
||
failures += [(res["sample_id"], *f) for f in res["merge_failures"]]
|
||
for code, n in sorted(reasons.items()):
|
||
if code.startswith("W_"):
|
||
rep.warn(code, f"{n} genotype(s): {code}")
|
||
write_merge_audit(opts.out / "audit" / "merges.parquet", failures)
|
||
ok = sum("error" not in r for r in results)
|
||
rm.summary = {"samples": len(results), "succeeded": ok, "records_per_sample": len(table.prefix),
|
||
"skipped_loci_by_code": dict(skips), "genotype_notes": dict(reasons)} # fmt: skip
|
||
for o in outs:
|
||
if o.exists():
|
||
rm.add_output(o, "vcf/4.3")
|
||
rm.status = "completed"
|
||
rm.write(opts.out / "run_manifest.json")
|
||
rep.info(f"done: {ok}/{len(results)} VCF(s) in {opts.out}")
|
||
if ok == len(results):
|
||
return ExitCode.OK
|
||
return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED
|
||
|
||
|
||
def write_audit(path: Path, audit: list[AuditRow]) -> None:
|
||
import pyarrow as pa
|
||
|
||
from simurg.core.schemas import write_parquet
|
||
|
||
cols = list(AuditRow.__dataclass_fields__)
|
||
table = pa.table({c: [getattr(a, c) for a in audit] for c in cols})
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
write_parquet(table, path)
|
||
|
||
|
||
def write_merge_audit(path: Path, rows: list[tuple[str, str, tuple[int, ...], str]]) -> None:
|
||
import pyarrow as pa
|
||
|
||
from simurg.core.schemas import write_parquet
|
||
|
||
table = pa.table(
|
||
{
|
||
"sample_id": [r[0] for r in rows],
|
||
"site_id": [r[1] for r in rows],
|
||
"calls": ["/".join(str(c) for c in r[2]) for r in rows],
|
||
"reason": [r[3] for r in rows],
|
||
}
|
||
)
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
write_parquet(table, path)
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ cohort BCF
|
||
def _f32_payload(v: npt.NDArray[np.float32]) -> bytes:
|
||
u = np.asarray(v, dtype="<f4").view("<u4").copy()
|
||
u[np.isnan(np.asarray(v, dtype=np.float32))] = 0x7F800001
|
||
return u.tobytes()
|
||
|
||
|
||
def write_cohort_bcf(
|
||
t: SiteTable, samples: list[tuple[str, str, str]], out: Path, opts: VcfOptions, rep: Reporter
|
||
) -> int:
|
||
"""OUT-040: one multi-sample BCF (+ .csi), same records and FORMAT rules as the per-sample VCFs."""
|
||
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
|
||
|
||
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],
|
||
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:
|
||
if s.error:
|
||
raise SimurgError("E_GTC_UNSUITABLE", f"{s.sample_id}: {s.error}")
|
||
names = [name for _, name, _ in samples]
|
||
cluster = read_cluster_name(samples[0][2])
|
||
header = t.header.replace("__CLUSTER__", cluster)
|
||
header += "".join(f"##SAMPLE=<ID={name},SentrixID={sid}>\n" for sid, name, _ in samples)
|
||
header += "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" + "\t".join(names) + "\n"
|
||
strings, contigs = header_dictionaries(header)
|
||
k_gt, k_gs, k_baf, k_lrr = strings["GT"], strings["GS"], strings["BAF"], strings["LRR"]
|
||
f_pass, f_low = strings["PASS"], strings.get("LowGS", -1)
|
||
n = len(names)
|
||
sexes = [s.sex_est for s in infos]
|
||
ploidy = np.stack([ploidy_for(t.site_class, sx, t.ploidy_policy) for sx in sexes], axis=1) # sites × samples
|
||
single_pos = np.full(len(t.prefix), -1, np.int64)
|
||
single_pos[t.single_rows] = np.arange(len(t.single_rows))
|
||
merged_pos = {int(k): j for j, k in enumerate(t.merged_rows.tolist())}
|
||
csi = CsiBuilder(len(contigs))
|
||
fdesc = _type_byte(1, FLOAT)
|
||
written = 0
|
||
with atomic_open(out, "wb") as raw:
|
||
w = BcfWriter(raw, header)
|
||
chunk = 20000
|
||
for start in range(0, len(t.prefix), chunk):
|
||
rows = range(start, min(len(t.prefix), start + chunk))
|
||
need = {int(t.single_locus[single_pos[k]]) for k in rows if single_pos[k] >= 0}
|
||
for k in rows:
|
||
if k in merged_pos:
|
||
need.update(p.index for p in t.merged_sites[merged_pos[k]].probes)
|
||
loci = sorted(need)
|
||
col = {loc: j for j, loc in enumerate(loci)}
|
||
sel = np.array(loci, np.int64)
|
||
gt = np.asarray(store.gt[:, sel]) if len(sel) else np.zeros((n, 0), np.uint8)
|
||
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:
|
||
j = int(single_pos[k])
|
||
c = col[int(t.single_locus[j])]
|
||
calls = gt[:, c].astype(np.int64)
|
||
ia, ib = (int(x) for x in t.single_ab[j]) # type: ignore[index]
|
||
a1 = np.where(calls == 3, ib, ia)
|
||
a2 = np.where(calls == 1, ia, ib)
|
||
lo, hi = np.minimum(a1, a2), np.maximum(a1, a2)
|
||
width = max(1, int(pl.max()))
|
||
cells = np.full((n, width), -127, np.int8)
|
||
miss = (calls == 0) | (pl == 0) | ((pl == 1) & (calls == 2))
|
||
hap = pl <= 1
|
||
cells[:, 0] = np.where(miss, 0, ((np.where(hap & (calls == 3), ib, lo) + 1) << 1)).astype(
|
||
np.int8
|
||
)
|
||
if width == 2:
|
||
second = np.where(miss, 0, (hi + 1) << 1).astype(np.int8)
|
||
cells[:, 1] = np.where(hap, np.int8(-127), second)
|
||
gs = sc[:, c]
|
||
b = bf[:, c]
|
||
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("<i4").tobytes()
|
||
fmt.append((strings[key], _type_byte(1, INT32), payload))
|
||
else:
|
||
fmt.append((strings[key], fdesc, _f32_payload(xsrc[key].astype(np.float32))))
|
||
gs_min = gs
|
||
else:
|
||
site = t.merged_sites[merged_pos[k]]
|
||
cs = [col[p.index] for p in site.probes]
|
||
g_all = gt[:, cs].astype(np.int64)
|
||
gs_min = sc[:, cs].min(axis=1)
|
||
width = max(1, int(pl.max()))
|
||
cells = np.full((n, width), -127, np.int8)
|
||
for si in range(n):
|
||
p_k = int(pl[si])
|
||
if p_k == 0:
|
||
cells[si, 0] = 0
|
||
continue
|
||
o = combine_genotype(site, tuple(int(x) for x in g_all[si]), p_k, t.inf1_nc_informative)
|
||
if o.alleles is None:
|
||
cells[si, :p_k] = 0
|
||
else:
|
||
for q, a in enumerate(o.alleles):
|
||
cells[si, q] = (a + 1) << 1
|
||
fmt = [(k_gt, _type_byte(width, INT8), cells.tobytes()), (k_gs, fdesc, _f32_payload(gs_min))]
|
||
low = opts.filter_gs is not None and bool(np.nanmin(gs_min) < opts.filter_gs) if n else False
|
||
info = []
|
||
for key, kind, vals in t.site_info_typed[k] if t.site_info_typed else []:
|
||
if kind == "Integer":
|
||
enc = typed_ints([int(v) for v in vals])
|
||
elif kind == "Float":
|
||
enc = typed_floats([float(v) for v in vals])
|
||
else:
|
||
enc = typed_string(",".join(str(v) for v in vals))
|
||
info.append((strings[key], enc))
|
||
rec = encode_record(contigs[t.site_contig[k]], int(t.site_beg[k]),
|
||
int(t.site_end[k] - t.site_beg[k]), t.site_id[k], t.site_alleles[k],
|
||
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
|
||
w.close()
|
||
Path(str(out) + ".csi").write_bytes(csi.encode())
|
||
finally:
|
||
store.close()
|
||
return written
|
||
|
||
|
||
def read_cluster_name(gtc_path: str) -> str:
|
||
from simurg.formats.gtc import read_gtc
|
||
|
||
return read_gtc(gtc_path, check_write_complete=False).cluster_file or ""
|