#!/usr/bin/env python3 """Run PLINK-derived kinship/PCA and EMMAX association atomically.""" from __future__ import annotations import argparse import csv import fcntl import gzip import hashlib import json import math import os import platform import re import shlex import shutil import subprocess import sys import time import uuid from dataclasses import dataclass from datetime import datetime, timezone from pathlib import Path from typing import Any, Iterable, Optional, Sequence import numpy as np from gfa_core.gwas_scientific_identity import build_identity_recipe, materialize_stage_identity, new_run_id STAGE = "gwas_emmax" METHOD_VERSION = "plink_emmax_pca.v1" MARKER_ID_SCHEME = "sha256_normalized_variant_tuple.v1.prefix24" PSD_RELATIVE_TOLERANCE = 0.01 RANK_RELATIVE_TOLERANCE = 1e-10 KNOWN_EMMAX_SHA256 = "ae2a0ffa285b846b785efe68f0cefeb9bcd4b2a288e53e8368fd6ef01dd8a76b" ESTIMABILITY_POLICY = "at_least_two_observed_diploid_genotype_states.v1" def sha256_file(path: Path) -> str: digest = hashlib.sha256() with path.open("rb") as handle: for chunk in iter(lambda: handle.read(1024 * 1024), b""): digest.update(chunk) return digest.hexdigest() def _write_json(path: Path, value: dict[str, Any]) -> None: path.write_text(json.dumps(value, ensure_ascii=False, indent=2) + "\n", encoding="utf-8") def _open_text(path: Path): return gzip.open(path, "rt", encoding="utf-8") if path.suffix == ".gz" else path.open("r", encoding="utf-8") @dataclass(frozen=True) class CanonicalVcf: samples: list[str] marker_map: list[dict[str, Any]] input_variant_count: int nonbiallelic_excluded_count: int half_call_count: int @dataclass(frozen=True) class Relationship: matrix: np.ndarray sample_ids: list[str] eigenvalues: np.ndarray minimum_eigenvalue: float maximum_eigenvalue: float tolerance: float marker_rank: int diagonal_min: float diagonal_max: float def _canonical_id(chrom: str, pos: int, ref: str, alt: str) -> str: payload = json.dumps([chrom, pos, ref, alt], ensure_ascii=False, separators=(",", ":")) return "v_" + hashlib.sha256(payload.encode()).hexdigest()[:24] def _canonicalize_vcf(source: Path, destination: Path) -> CanonicalVcf: samples: Optional[list[str]] = None marker_map: list[dict[str, Any]] = [] tuples: set[tuple[str, int, str, str]] = set() ids: dict[str, tuple[str, int, str, str]] = {} input_count = excluded = half_calls = 0 destination.parent.mkdir(parents=True, exist_ok=True) with _open_text(source) as incoming, destination.open("w", encoding="utf-8") as outgoing: for line_number, raw in enumerate(incoming, 1): if raw.startswith("##"): outgoing.write(raw) continue if raw.startswith("#CHROM\t"): if samples is not None: raise ValueError("sorted VCF contains multiple #CHROM headers") fields = raw.rstrip("\r\n").split("\t") if len(fields) < 10: raise ValueError("sorted VCF must contain genotype samples") samples = fields[9:] if any(not sample for sample in samples) or len(samples) != len(set(samples)): raise ValueError("sorted VCF sample IDs must be nonempty and unique") outgoing.write(raw) continue if not raw.strip() or raw.startswith("#"): continue if samples is None: raise ValueError("sorted VCF #CHROM header is missing") fields = raw.rstrip("\r\n").split("\t") if len(fields) != 9 + len(samples): raise ValueError(f"malformed sorted VCF variant at line {line_number}") input_count += 1 chrom, pos_text, original_id, ref, alt = fields[:5] try: pos = int(pos_text) except ValueError as exc: raise ValueError(f"VCF position must be an integer at line {line_number}") from exc chrom, ref, alt = chrom.strip(), ref.strip().upper(), alt.strip().upper() if not chrom or pos <= 0 or not ref or not alt or alt == ".": raise ValueError(f"VCF variant tuple is invalid at line {line_number}") normalized = (chrom, pos, ref, alt) if normalized in tuples: raise ValueError("duplicate normalized variant tuple") tuples.add(normalized) format_fields = fields[8].split(":") if "GT" in format_fields: gt_index = format_fields.index("GT") for sample_field in fields[9:]: values = sample_field.split(":") gt = values[gt_index].replace("|", "/") if gt_index < len(values) else "." alleles = gt.split("/") if len(alleles) >= 2 and any(value == "." for value in alleles) and not all(value == "." for value in alleles): half_calls += 1 if "," in alt: excluded += 1 continue canonical = _canonical_id(chrom, pos, ref, alt) if canonical in ids and ids[canonical] != normalized: raise ValueError("canonical marker hash-prefix collision") ids[canonical] = normalized marker_map.append({"canonical_id": canonical, "original_id": original_id, "chrom": chrom, "pos": pos, "ref": ref, "alt": alt}) fields[2], fields[3], fields[4] = canonical, ref, alt outgoing.write("\t".join(fields) + "\n") if samples is None: raise ValueError("sorted VCF #CHROM header is missing") if input_count == 0 or not marker_map: raise ValueError("sorted VCF contains no eligible biallelic variants") return CanonicalVcf(samples, marker_map, input_count, excluded, half_calls) def _read_trait(path: Path) -> dict[str, Any]: lines = [line for line in path.read_text(encoding="utf-8").splitlines() if line.strip()] if not lines or not lines[0].startswith("\t") or len(lines[0].split("\t")) != 2: raise ValueError("prepared trait table header is invalid") trait_name = lines[0].split("\t")[1] values: dict[str, float] = {} for line in lines[1:]: fields = line.split("\t") if len(fields) != 2 or not fields[0]: raise ValueError("prepared trait table row is invalid") if fields[0] in values: raise ValueError("prepared trait sample IDs must be unique") try: value = float(fields[1]) except ValueError as exc: raise ValueError("prepared trait values must be numeric") from exc if not math.isfinite(value): raise ValueError("prepared trait values must be finite") values[fields[0]] = value if not values: raise ValueError("prepared trait table contains no samples") if min(values.values()) == max(values.values()): raise ValueError("phenotype variance must be positive") return {"trait_name": trait_name, "values": values} def _validate_trait_subset(trait: dict[str, Any], genotype_samples: Sequence[str]) -> list[str]: missing = sorted(set(trait["values"]) - set(genotype_samples)) if missing: raise ValueError("phenotype samples absent from genotype: " + ",".join(missing)) return [sample for sample in genotype_samples if sample not in trait["values"]] def _ordered_phenotype(trait: dict[str, Any], sample_order: Sequence[str]) -> list[tuple[str, float]]: missing = [sample for sample in sample_order if sample not in trait["values"]] if missing or set(sample_order) != set(trait["values"]): raise ValueError("TFAM sample set does not match phenotype sample set") return [(sample, float(trait["values"][sample])) for sample in sample_order] def _read_two_id_file(path: Path, label: str) -> list[str]: rows: list[str] = [] for line in path.read_text(encoding="utf-8").splitlines(): if not line.strip(): continue fields = line.split() if len(fields) < 2 or fields[0] != fields[1]: raise ValueError(f"{label} must contain matching FID/IID") rows.append(fields[1]) if not rows or len(rows) != len(set(rows)): raise ValueError(f"{label} sample IDs must be nonempty and unique") return rows def _read_relationship(rel_path: Path, rel_id_path: Path, expected_order: Sequence[str]) -> Relationship: ids = _read_two_id_file(rel_id_path, "REL-ID") if ids != list(expected_order): raise ValueError("REL-ID order does not match TFAM order") lines = [line for line in rel_path.read_text(encoding="utf-8").splitlines() if line.strip()] n = len(ids) if len(lines) != n: raise ValueError("relationship lower triangle row count is invalid") matrix = np.zeros((n, n), dtype=float) for i, line in enumerate(lines): fields = line.split() if len(fields) != i + 1: raise ValueError("relationship lower triangle width is invalid") try: values = [float(value) for value in fields] except ValueError as exc: raise ValueError("relationship matrix contains nonnumeric values") from exc if not all(math.isfinite(value) for value in values): raise ValueError("relationship matrix contains nonfinite values") matrix[i, : i + 1] = values matrix[: i + 1, i] = values if not np.allclose(matrix, matrix.T, rtol=0, atol=1e-12): raise ValueError("relationship matrix is asymmetric") eigenvalues = np.linalg.eigvalsh(matrix) scale = max(float(np.max(np.abs(eigenvalues))), np.finfo(float).eps) tolerance = PSD_RELATIVE_TOLERANCE * scale if float(eigenvalues[0]) < -tolerance: raise ValueError("relationship matrix exceeds positive-semidefinite tolerance") rank_tolerance = RANK_RELATIVE_TOLERANCE * scale rank = int(np.count_nonzero(eigenvalues > rank_tolerance)) diagonal = np.diag(matrix) return Relationship(matrix, ids, eigenvalues, float(eigenvalues[0]), float(eigenvalues[-1]), tolerance, rank, float(diagonal.min()), float(diagonal.max())) def _realized_pc_count(*, requested: int, sample_count: int, marker_rank: int) -> int: if requested < 0 or sample_count < 3 or marker_rank < 0: raise ValueError("PCA dimensions are invalid") return min(requested, sample_count - 2, marker_rank) def _marker_matrix_rank(tped_path: Path, sample_count: int) -> tuple[int, np.ndarray]: """Return rank of the centered genotype matrix without materializing M x N.""" if sample_count < 1: raise ValueError("TPED sample count is invalid") gram = np.zeros((sample_count, sample_count), dtype=float) marker_count = 0 for line_number, line in enumerate(tped_path.read_text(encoding="utf-8").splitlines(), 1): if not line.strip(): continue fields = line.split() if len(fields) != 4 + 2 * sample_count: raise ValueError(f"TPED genotype width is invalid at line {line_number}") dosage = np.empty(sample_count, dtype=float) dosage.fill(np.nan) for index in range(sample_count): alleles = fields[4 + 2 * index: 6 + 2 * index] if "0" in alleles: continue if any(allele not in {"1", "2"} for allele in alleles): raise ValueError(f"TPED genotype allele is invalid at line {line_number}") dosage[index] = float(alleles.count("2")) observed = np.isfinite(dosage) if not observed.any(): continue mean = float(np.mean(dosage[observed])) dosage[~observed] = mean dosage -= mean gram += np.outer(dosage, dosage) marker_count += 1 if marker_count == 0: raise ValueError("TPED contains no usable genotype markers") eigenvalues = np.linalg.eigvalsh(gram) scale = max(float(np.max(np.abs(eigenvalues))), np.finfo(float).eps) rank = int(np.count_nonzero(eigenvalues > RANK_RELATIVE_TOLERANCE * scale)) return rank, eigenvalues def _filter_tped_estimable_markers( tped_path: Path, *, sample_count: int, marker_map: dict[str, dict[str, Any]] ) -> tuple[dict[str, dict[str, Any]], dict[str, Any]]: """Atomically retain TPED markers with at least two observed genotype states.""" if sample_count < 1: raise ValueError("TPED sample count is invalid") seen: set[str] = set() retained_ids: list[str] = [] retained_lines: list[str] = [] no_called = nonvariable = 0 for line_number, line in enumerate(tped_path.read_text(encoding="utf-8").splitlines(), 1): if not line.strip(): continue fields = line.split() if len(fields) != 4 + 2 * sample_count: raise ValueError(f"TPED genotype width is invalid at line {line_number}") marker_id = fields[1] if marker_id in seen: raise ValueError("TPED marker IDs must be unique") seen.add(marker_id) if marker_id not in marker_map: raise ValueError(f"TPED contains unexpected marker ID: {marker_id}") states: set[tuple[str, str]] = set() for index in range(sample_count): first, second = fields[4 + 2 * index: 6 + 2 * index] if (first == "0") != (second == "0"): raise ValueError(f"TPED contains a half-call at line {line_number}") if first == second == "0": continue if first not in {"1", "2"} or second not in {"1", "2"}: raise ValueError(f"TPED genotype allele is invalid at line {line_number}") states.add(tuple(sorted((first, second)))) if not states: no_called += 1 elif len(states) == 1: nonvariable += 1 else: retained_ids.append(marker_id) retained_lines.append(line) if not retained_ids: raise ValueError("TPED contains no estimable markers") temporary = tped_path.with_name(f".{tped_path.name}.{uuid.uuid4().hex}.estimability") try: temporary.write_text("\n".join(retained_lines) + "\n", encoding="utf-8") os.replace(temporary, tped_path) finally: temporary.unlink(missing_ok=True) candidate_count = len(seen) tested_count = len(retained_ids) if candidate_count != no_called + nonvariable + tested_count: raise AssertionError("TPED estimability accounting is inconsistent") return ( {marker_id: marker_map[marker_id] for marker_id in retained_ids}, { "emmax_candidate_marker_count": candidate_count, "no_called_genotype_excluded_count": no_called, "nonvariable_genotype_excluded_count": nonvariable, "emmax_tested_count": tested_count, "estimability_policy": ESTIMABILITY_POLICY, }, ) def _validate_design(eigenvectors: np.ndarray, *, realized_pcs: int) -> tuple[int, int]: if eigenvectors.ndim != 2 or eigenvectors.shape[1] < realized_pcs: raise ValueError("PCA eigenvector dimensions are invalid") design = np.column_stack([np.ones(eigenvectors.shape[0]), eigenvectors[:, :realized_pcs]]) rank = int(np.linalg.matrix_rank(design)) if rank != design.shape[1]: raise ValueError("intercept plus PCA covariates must have full rank") residual_df = eigenvectors.shape[0] - rank if residual_df <= 0: raise ValueError("model must retain positive residual degrees of freedom") return rank, residual_df def _run_required_tool(command: list[str], label: str, *, required_output: Optional[Path] = None, cwd: Optional[Path] = None, evidence_sink: Optional[list[str]] = None) -> str: completed = subprocess.run(command, cwd=cwd, text=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, check=False) output = completed.stdout or "" if evidence_sink is not None: evidence_sink.append(f"[tool:{label}]\n{output.rstrip()}\n") if completed.returncode != 0: raise RuntimeError(f"{label} failed with exit code {completed.returncode}") if required_output is not None and (not required_output.is_file() or required_output.stat().st_size == 0): raise RuntimeError(f"{label} completed without required output") return output def _parse_emmax_ps(source: Path, marker_map: dict[str, dict[str, Any]], full_path: Path, pvalue_path: Path) -> dict[str, Any]: full = ["SNPID\tChr\tPos\tbeta\tSE\tp"] pvalues = ["chr\tpos\tp"] seen: set[str] = set() source_field_count: Optional[int] = None for line_number, line in enumerate(source.read_text(encoding="utf-8").splitlines(), 1): if not line.strip(): continue fields = line.split() if len(fields) not in {3, 4}: raise ValueError("EMMAX .ps rows must contain exactly three or four fields") if source_field_count is None: source_field_count = len(fields) elif len(fields) != source_field_count: raise ValueError("EMMAX .ps must use one consistent three- or four-field schema") if source_field_count == 3: marker, beta_text, p_text = fields se_text: Optional[str] = None else: marker, beta_text, se_text, p_text = fields if marker in seen: raise ValueError("EMMAX .ps contains duplicate SNPID") if marker not in marker_map: raise ValueError(f"EMMAX .ps contains unexpected SNPID: {marker}") try: beta, pvalue = float(beta_text), float(p_text) se = float(se_text) if se_text is not None else None except ValueError as exc: raise ValueError("EMMAX beta/SE/p must be numeric") from exc if not math.isfinite(beta) or not math.isfinite(pvalue) or (se is not None and not math.isfinite(se)): raise ValueError("EMMAX beta/SE/p must be finite") if se is not None and se < 0: raise ValueError("EMMAX SE must be nonnegative") if not 0 < pvalue <= 1: raise ValueError("EMMAX p-value must satisfy 0 < p <= 1") seen.add(marker) record = marker_map[marker] rendered_se = se_text if se_text is not None else "NA" full.append(f"{marker}\t{record['chrom']}\t{record['pos']}\t{beta_text}\t{rendered_se}\t{p_text}") pvalues.append(f"{record['chrom']}\t{record['pos']}\t{p_text}") if not seen: raise ValueError("EMMAX .ps contains no association rows") missing = set(marker_map) - seen if missing: raise ValueError(f"EMMAX .ps is missing {len(missing)} mapped SNPIDs") full_path.write_text("\n".join(full) + "\n", encoding="utf-8") pvalue_path.write_text("\n".join(pvalues) + "\n", encoding="utf-8") assert source_field_count in {3, 4} return {"emmax_tested_count": len(seen), "mapped_count": len(seen), "duplicate_count": 0, "unmapped_count": 0, "unexpected_count": 0, "finite_beta_count": len(seen), "finite_se_count": len(seen) if source_field_count == 4 else 0, "missing_se_count": len(seen) if source_field_count == 3 else 0, "valid_p_count": len(seen), "source_ps_field_count": source_field_count, "source_ps_schema": "SNPID_beta_p" if source_field_count == 3 else "SNPID_beta_SE_p"} def _safe_prefix(prefix: str) -> str: return re.sub(r"[^A-Za-z0-9._-]+", "_", prefix).strip("._") or "cohort" def _runner_command(args: argparse.Namespace) -> list[str]: command = ["python", "-m", "cli.run_gwas_emmax", "--workspace", str(args.workspace.resolve()), "--job-id", args.job_id, "--prefix", args.prefix, "--vcf", str(args.vcf.resolve()), "--trait", str(args.trait.resolve()), "--vcf-manifest", str(args.vcf_manifest.resolve()), "--trait-manifest", str(args.trait_manifest.resolve()), "--genotype-scope-id", args.genotype_scope_id, "--phenotype-scope-id", args.phenotype_scope_id, "--plink", args.plink, "--emmax", args.emmax, "--rscript", args.rscript, "--requested-pcs", str(args.requested_pcs), "--threads", str(args.threads)] if args.log_file: command.extend(["--log-file", str(args.log_file.resolve())]) return command def _artifact(workspace: Path, source: Path, destination: Path, role: str, media_type: str) -> dict[str, Any]: size = source.stat().st_size return {"artifact": role, "role": role, "path": str(destination.relative_to(workspace)), "full_path": str(destination.resolve()), "type": media_type, "media_type": media_type, "size": size, "size_bytes": size, "sha256": sha256_file(source)} def _authoritative_input_record(path: Path, role: str, producer_manifest_path: Path, producer_stage: str, scope_field: str, expected_scope: str, require_analysis_scope: bool = False) -> dict[str, Any]: path = path.resolve() producer_manifest_path = producer_manifest_path.resolve() if not path.is_file() or not producer_manifest_path.is_file(): raise FileNotFoundError(f"authoritative {role} or producer Manifest does not exist") try: producer = json.loads(producer_manifest_path.read_text(encoding="utf-8")) except json.JSONDecodeError as exc: raise ValueError(f"{producer_stage} producer Manifest is invalid JSON") from exc if producer.get("manifest_version") != "1.0": raise ValueError(f"{producer_stage} producer Manifest version must be 1.0") if producer.get("stage") != producer_stage or producer.get("status") != "success": raise ValueError(f"{role} producer Manifest stage/status is invalid") run_id = str(producer.get("run_id") or "") if not re.fullmatch(r"[A-Za-z0-9][A-Za-z0-9._:-]{0,255}", run_id): raise ValueError(f"{producer_stage} producer Manifest run_id is empty or unsafe") analysis_scope = str(producer.get("analysis_scope_id") or "") if ((require_analysis_scope and not analysis_scope) or (analysis_scope and not re.fullmatch(r"gwas-analysis:[0-9a-f]{64}", analysis_scope))): raise ValueError(f"{producer_stage} producer Manifest analysis_scope_id is invalid") if producer.get(scope_field) != expected_scope: raise ValueError(f"{role} producer Manifest {scope_field} does not match") artifact = (producer.get("artifacts") or {}).get(role) if not isinstance(artifact, dict) or artifact.get("role") != role or artifact.get("artifact") != role: raise ValueError(f"{role} producer Manifest artifact role is invalid") if Path(str(artifact.get("full_path") or "")).resolve() != path: raise ValueError(f"{role} producer Manifest artifact path does not match") size = path.stat().st_size digest = sha256_file(path) if artifact.get("size_bytes", artifact.get("size")) != size: raise ValueError(f"{role} producer Manifest artifact size does not match") if artifact.get("sha256") != digest: raise ValueError(f"{role} producer Manifest artifact SHA256 does not match") record = {"source": "manifest", "artifact": role, "role": role, "path": str(path), "full_path": str(path), "size": size, "size_bytes": size, "sha256": digest, "producer_stage": producer_stage, "producer_run_id": run_id, "producer_manifest_path": str(producer_manifest_path), "producer_manifest_sha256": sha256_file(producer_manifest_path), "producer_scope_field": scope_field, "producer_scope_id": expected_scope} if analysis_scope: record["producer_analysis_scope_id"] = analysis_scope return record def _revalidate_input_records(records: Sequence[dict[str, Any]]) -> None: for record in records: path = Path(str(record["full_path"])).resolve() if (not path.is_file() or path.stat().st_size != record["size_bytes"] or sha256_file(path) != record["sha256"]): raise RuntimeError(f"authoritative {record['artifact']} changed during execution") manifest_path = Path(str(record["producer_manifest_path"])).resolve() if (not manifest_path.is_file() or sha256_file(manifest_path) != record["producer_manifest_sha256"]): raise RuntimeError(f"authoritative {record['artifact']} producer Manifest changed during execution") def _warning_lines(outputs: Iterable[str]) -> list[str]: return [line.strip() for output in outputs for line in output.splitlines() if "warning" in line.lower()] def _write_reproducibility(workspace: Path, identity12: str, run_id: str, command: list[str], manifest_path: Path, manifest: dict[str, Any], environment: str, artifact_sources: Optional[dict[str, Path]] = None, public_manifest_path: Optional[Path] = None) -> Path: bundle = workspace / "reproducibility" / STAGE / f"{_safe_prefix(manifest['prefix'])}.{identity12}" / run_id bundle.mkdir(parents=True, exist_ok=False) canonical = shlex.join(command) (bundle / "commands.sh").write_text("#!/usr/bin/env bash\nset -euo pipefail\n" + canonical + "\n", encoding="utf-8") (bundle / "commands.sh").chmod(0o755) (bundle / "environment.txt").write_text(environment, encoding="utf-8") recorded_manifest = public_manifest_path or manifest_path _write_json(bundle / "progress.json", {"stage": STAGE, "status": manifest["status"], "run_id": run_id, "manifest_path": str(recorded_manifest.relative_to(workspace))}) _write_json(bundle / "provenance.json", {"stage": STAGE, "run_id": run_id, "runner_command": canonical, "software": manifest.get("software", {}), "inputs": manifest.get("inputs", []), "artifacts": manifest.get("artifacts", {})}) checksums = [(manifest_path, str(recorded_manifest.relative_to(workspace)))] for role, artifact in manifest.get("artifacts", {}).items(): source = artifact_sources.get(role) if artifact_sources else workspace / artifact["path"] if source and source.is_file(): checksums.append((source, artifact["path"])) (bundle / "checksums.sha256").write_text("".join(f"{sha256_file(path)} {label}\n" for path, label in sorted(checksums, key=lambda item: item[1])), encoding="utf-8") return bundle def _publish(sources: dict[str, Path], destinations: dict[str, Path], manifest_source: Path, manifest_destination: Path, lock_path: Path, started_at: str) -> bool: lock_path.parent.mkdir(parents=True, exist_ok=True) with lock_path.open("a+") as lock: fcntl.flock(lock.fileno(), fcntl.LOCK_EX) if manifest_destination.is_file(): try: previous = json.loads(manifest_destination.read_text(encoding="utf-8")) except (OSError, json.JSONDecodeError): previous = {} if previous.get("status") == "success" and str(previous.get("start_time", "")) > started_at: return False token = uuid.uuid4().hex staged: list[tuple[Path, Path]] = [] for role, source in sources.items(): destination = destinations[role] destination.parent.mkdir(parents=True, exist_ok=True) temporary = destination.with_name(f".{destination.name}.{token}.publishing") shutil.copy2(source, temporary) staged.append((temporary, destination)) temporary_manifest = manifest_destination.with_name(f".{manifest_destination.name}.{token}.publishing") shutil.copy2(manifest_source, temporary_manifest) backups: dict[Path, Optional[Path]] = {} try: for _temporary, destination in [*staged, (temporary_manifest, manifest_destination)]: if destination.is_file(): backup = destination.with_name(f".{destination.name}.{token}.rollback") shutil.copy2(destination, backup) backups[destination] = backup else: backups[destination] = None for temporary, destination in staged: os.replace(temporary, destination) os.replace(temporary_manifest, manifest_destination) except Exception: for destination, backup in backups.items(): if backup is None: destination.unlink(missing_ok=True) elif backup.is_file(): os.replace(backup, destination) for temporary, _destination in staged: temporary.unlink(missing_ok=True) temporary_manifest.unlink(missing_ok=True) for backup in backups.values(): if backup is not None: backup.unlink(missing_ok=True) raise for backup in backups.values(): if backup is not None: backup.unlink(missing_ok=True) return True def _publish_failed_manifest(path: Path, manifest: dict[str, Any], lock_path: Path, success_manifest_path: Optional[Path] = None) -> bool: lock_path.parent.mkdir(parents=True, exist_ok=True) with lock_path.open("a+") as lock: fcntl.flock(lock.fileno(), fcntl.LOCK_EX) started_at = str(manifest.get("start_time") or "") for candidate in (success_manifest_path, path): if candidate is None or not candidate.is_file(): continue try: existing = json.loads(candidate.read_text(encoding="utf-8")) except (OSError, json.JSONDecodeError): continue if str(existing.get("start_time") or "") > started_at: return False temporary = path.with_name(f".{path.name}.{uuid.uuid4().hex}.publishing") _write_json(temporary, manifest) os.replace(temporary, path) return True def _clear_failed_manifest_if_winner(failed_path: Path, lock_path: Path, success_manifest_path: Path, run_id: str) -> bool: lock_path.parent.mkdir(parents=True, exist_ok=True) with lock_path.open("a+") as lock: fcntl.flock(lock.fileno(), fcntl.LOCK_EX) try: success = json.loads(success_manifest_path.read_text(encoding="utf-8")) except (OSError, json.JSONDecodeError): return False if success.get("status") != "success" or success.get("run_id") != run_id: return False if not failed_path.is_file(): return True try: failed = json.loads(failed_path.read_text(encoding="utf-8")) except (OSError, json.JSONDecodeError): return False if str(failed.get("start_time") or "") > str(success.get("start_time") or ""): return False failed_path.unlink() return True def _publish_with_bundle(sources: dict[str, Path], destinations: dict[str, Path], manifest_source: Path, manifest_destination: Path, lock_path: Path, started_at: str, bundle: Path) -> bool: try: published = _publish(sources, destinations, manifest_source, manifest_destination, lock_path, started_at) except Exception: shutil.rmtree(bundle, ignore_errors=True) raise if not published: # Its checksum labels target the public paths, but those paths belong to # the newer winner. Keeping this bundle would make replay verification # falsely compare the older attempt hash with the newer Manifest. shutil.rmtree(bundle, ignore_errors=True) return published def _tool_version(command: list[str], label: str) -> str: completed = subprocess.run(command, text=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, check=False) text = (completed.stdout or "").strip().splitlines() if not text: raise RuntimeError(f"{label} version could not be proven") return text[0].strip() def _emmax_provenance(path: Path) -> dict[str, str]: executable_hash = sha256_file(path) if executable_hash == KNOWN_EMMAX_SHA256: return {"release": "20120210", "sha256": executable_hash, "package": "upstream-binary"} owner = subprocess.run(["dpkg-query", "-S", str(path)], text=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, check=False) if owner.returncode != 0 or ":" not in (owner.stdout or ""): raise RuntimeError("EMMAX executable provenance could not be proven") package = (owner.stdout or "").split(":", 1)[0].strip() if package != "emmax": raise RuntimeError("EMMAX executable is not owned by the official emmax package") version = subprocess.run(["dpkg-query", "-W", "-f=${Package}\t${Version}\\n", package], text=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, check=False) fields = (version.stdout or "").strip().split("\t") if version.returncode != 0 or len(fields) != 2 or fields[0] != package or not fields[1]: raise RuntimeError("EMMAX package version could not be proven") return {"release": fields[1], "sha256": executable_hash, "package": package} def _write_matrix_via_r(rscript: Path, relationship: Relationship, attempt: Path, evidence_sink: Optional[list[str]] = None) -> Path: source = attempt / "write_kinship.R" input_path = attempt / "kinship.python.tsv" output_path = attempt / "pop.kinship.txt" np.savetxt(input_path, relationship.matrix, delimiter="\t", fmt="%.10g") source.write_text("args <- commandArgs(trailingOnly=TRUE)\nx <- as.matrix(read.table(args[1], header=FALSE, check.names=FALSE))\nwrite.table(x, args[2], row.names=FALSE, col.names=FALSE, quote=FALSE, sep=' ')\n", encoding="utf-8") _run_required_tool([str(rscript), str(source), str(input_path), str(output_path)], "R kinship materialization", required_output=output_path, cwd=attempt, evidence_sink=evidence_sink) rendered = np.loadtxt(output_path, dtype=float, ndmin=2) if rendered.shape != relationship.matrix.shape or not np.allclose(rendered, relationship.matrix, rtol=1e-8, atol=1e-10): raise ValueError("R materialized relationship matrix does not preserve PLINK values") return output_path def _read_marker_ids(path: Path, column: int, label: str) -> list[str]: ids = [line.split()[column] for line in path.read_text(encoding="utf-8").splitlines() if line.strip()] if not ids or len(ids) != len(set(ids)): raise ValueError(f"{label} marker IDs must be nonempty and unique") return ids def _require_exact_marker_set(actual: Sequence[str], expected: set[str], label: str) -> None: if set(actual) != expected: raise ValueError(f"{label} is not authorized") def _read_eigenvectors(path: Path, expected_order: Sequence[str], realized_pcs: int) -> np.ndarray: ids: list[str] = [] rows: list[list[float]] = [] for line in path.read_text(encoding="utf-8").splitlines(): fields = line.split() if not fields: continue if len(fields) != 2 + realized_pcs or fields[0] != fields[1]: raise ValueError("PLINK PCA eigenvector row width or FID/IID is invalid") ids.append(fields[1]) try: rows.append([float(value) for value in fields[2:]]) except ValueError as exc: raise ValueError("PLINK PCA eigenvectors must be numeric") from exc if ids != list(expected_order): raise ValueError("PLINK PCA sample order does not match TFAM order") matrix = np.asarray(rows, dtype=float) if not np.isfinite(matrix).all(): raise ValueError("PLINK PCA eigenvectors must be finite") return matrix def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--workspace", required=True, type=Path) parser.add_argument("--job-id", required=True) parser.add_argument("--prefix", required=True) parser.add_argument("--vcf", required=True, type=Path) parser.add_argument("--trait", required=True, type=Path) parser.add_argument("--vcf-manifest", required=True, type=Path) parser.add_argument("--trait-manifest", required=True, type=Path) parser.add_argument("--genotype-scope-id", required=True) parser.add_argument("--phenotype-scope-id", required=True) parser.add_argument("--plink", required=True) parser.add_argument("--emmax", required=True) parser.add_argument("--rscript", required=True) parser.add_argument("--requested-pcs", type=int, default=3) parser.add_argument("--threads", type=int, default=1) parser.add_argument("--log-file", type=Path) return parser def main(argv: Optional[list[str]] = None) -> int: args = build_parser().parse_args(argv) workspace, vcf, trait_path = args.workspace.resolve(), args.vcf.resolve(), args.trait.resolve() final_dir = workspace / "final" final_dir.mkdir(parents=True, exist_ok=True) run_id = new_run_id(STAGE) attempt = workspace / ".attempts" / run_id attempt.mkdir(parents=True, exist_ok=False) started = datetime.now(timezone.utc) clock = time.monotonic() safe_prefix = _safe_prefix(args.prefix) runner = _runner_command(args) parameters = {"requested_pcs": args.requested_pcs, "pc_selection_rule": "min(requested_pcs,N-2,marker_matrix_rank)", "multiallelic_policy": "exclude", "vcf_half_call_policy": "missing", "estimability_policy": ESTIMABILITY_POLICY, "marker_id_scheme": MARKER_ID_SCHEME, "emmax_precision_digits": 10, "psd_relative_tolerance": PSD_RELATIVE_TOLERANCE, "rank_relative_tolerance": RANK_RELATIVE_TOLERANCE} manifest: dict[str, Any] = {"manifest_version": "1.0", "stage": STAGE, "prefix": args.prefix, "job_id": args.job_id, "run_id": run_id, "status": "failed", "timestamp": started.isoformat(), "start_time": started.isoformat(), "end_time": started.isoformat(), "elapsed_sec": 0.0, "cmd": runner, "workspace_root": str(workspace), "workspace": str(final_dir), "input_dir": str(workspace / "input"), "parameters": parameters, "metrics": {}, "artifacts": {}} environment = "" identity12 = "failed" logs: list[str] = [] log_file = args.log_file.resolve() if args.log_file else workspace / "logs" / f"{run_id}.log" log_file.parent.mkdir(parents=True, exist_ok=True) try: if args.threads < 1 or args.threads > 64 or args.requested_pcs < 0: raise ValueError("threads must be 1..64 and requested PCs must be nonnegative") for name, value, prefix in (("genotype_scope_id", args.genotype_scope_id, "genotype"), ("phenotype_scope_id", args.phenotype_scope_id, "phenotype")): if not re.fullmatch(prefix + r":[0-9a-f]{64}", value): raise ValueError(f"{name} is invalid") if not vcf.is_file() or not trait_path.is_file(): raise FileNotFoundError("authoritative sorted_vcf and trait_table are required") tools = {name: Path(value).resolve() for name, value in (("plink", args.plink), ("emmax", args.emmax), ("rscript", args.rscript))} for name, path in tools.items(): if not path.is_file() or not os.access(path, os.X_OK): raise FileNotFoundError(f"{name} executable was not found or is not executable") plink_version = _tool_version([str(tools["plink"]), "--version"], "PLINK") r_version = _tool_version([str(tools["rscript"]), "--version"], "Rscript") emmax_provenance = _emmax_provenance(tools["emmax"]) emmax_hash = emmax_provenance["sha256"] emmax_release = emmax_provenance["release"] environment = "\n".join(["runner=cli.run_gwas_emmax.py", f"runtime_environment={os.environ.get('CONDA_DEFAULT_ENV', 'unknown')}", f"python_resolved={Path(sys.executable).resolve()}", f"python_version={platform.python_version()}", f"plink_resolved={tools['plink']}", f"plink_version={plink_version}", f"emmax_resolved={tools['emmax']}", f"emmax_release={emmax_release}", f"emmax_package={emmax_provenance['package']}", f"emmax_sha256={emmax_hash}", "emmax_ps_supported_schemas=SNPID_beta_p|SNPID_beta_SE_p", "emmax_ps_official_schema=SNPID_beta_p", "emmax_ps_missing_se_policy=NA", f"rscript_resolved={tools['rscript']}", f"r_version={r_version}", f"platform={platform.platform()}"]) + "\n" inputs = [ _authoritative_input_record(vcf, "sorted_vcf", args.vcf_manifest, "gwas_fmt", "genotype_scope_id", args.genotype_scope_id), _authoritative_input_record(trait_path, "trait_table", args.trait_manifest, "gwas_trait", "phenotype_scope_id", args.phenotype_scope_id, require_analysis_scope=True), ] manifest["inputs"] = inputs scopes = {"genotype_scope_id": args.genotype_scope_id, "phenotype_scope_id": args.phenotype_scope_id} recipe = build_identity_recipe(stage=STAGE, required_artifact_roles=["sorted_vcf", "trait_table"], parameters=parameters, method_version=METHOD_VERSION, contract_version="v1", tool_versions={"plink": plink_version, "emmax_release": emmax_release, "emmax_sha256": emmax_hash, "R": r_version, "python": platform.python_version()}) identity = materialize_stage_identity(recipe, {item["artifact"]: item for item in inputs}, scopes) manifest.update(identity) identity12 = identity["analysis_scope_id"].split(":", 1)[1][:12] canonical = _canonicalize_vcf(vcf, attempt / "canonical.vcf") trait = _read_trait(trait_path) genotype_only = _validate_trait_subset(trait, canonical.samples) keep = attempt / "trait.keep" keep.write_text("".join(f"{sample}\t{sample}\n" for sample in trait["values"]), encoding="utf-8") bed = attempt / "genotype" plink_outputs: list[str] = [] plink_outputs.append(_run_required_tool([str(tools["plink"]), "--vcf", str(attempt / "canonical.vcf"), "--keep", str(keep), "--vcf-half-call", "missing", "--biallelic-only", "strict", "--allow-extra-chr", "--double-id", "--make-bed", "--threads", str(args.threads), "--out", str(bed)], "PLINK VCF conversion", required_output=bed.with_suffix(".bed"), cwd=attempt, evidence_sink=logs)) fam_ids = _read_two_id_file(bed.with_suffix(".fam"), "FAM") if set(fam_ids) != set(trait["values"]): raise ValueError("PLINK retained FAM sample set does not match phenotype") marker_map = {row["canonical_id"]: row for row in canonical.marker_map} bim_ids = _read_marker_ids(bed.with_suffix(".bim"), 1, "BIM") _require_exact_marker_set(bim_ids, set(marker_map), "PLINK conversion marker loss") if canonical.input_variant_count != canonical.nonbiallelic_excluded_count + len(bim_ids): raise ValueError("PLINK variant loss accounting is inconsistent") rel = attempt / "relationship" plink_outputs.append(_run_required_tool([str(tools["plink"]), "--bfile", str(bed), "--make-rel", "--allow-extra-chr", "--out", str(rel)], "PLINK relationship", required_output=rel.with_suffix(".rel"), cwd=attempt, evidence_sink=logs)) relationship = _read_relationship(rel.with_suffix(".rel"), rel.with_suffix(".rel.id"), fam_ids) kinship = _write_matrix_via_r(tools["rscript"], relationship, attempt, logs) trans = attempt / "emmax_genotype" plink_outputs.append(_run_required_tool([str(tools["plink"]), "--bfile", str(bed), "--recode", "12", "transpose", "--output-missing-genotype", "0", "--allow-extra-chr", "--out", str(trans)], "PLINK transpose", required_output=trans.with_suffix(".tped"), cwd=attempt, evidence_sink=logs)) tfam_ids = _read_two_id_file(trans.with_suffix(".tfam"), "TFAM") if tfam_ids != fam_ids: raise ValueError("TFAM sample order does not match FAM order") tped_ids = _read_marker_ids(trans.with_suffix(".tped"), 1, "TPED") _require_exact_marker_set(tped_ids, set(marker_map), "TPED marker mismatch") tested_marker_map, estimability_metrics = _filter_tped_estimable_markers( trans.with_suffix(".tped"), sample_count=len(tfam_ids), marker_map=marker_map ) environment += "\n".join( f"{key}={value}" for key, value in estimability_metrics.items() ) + "\n" marker_matrix_rank, marker_matrix_eigenvalues = _marker_matrix_rank( trans.with_suffix(".tped"), len(tfam_ids) ) realized = _realized_pc_count(requested=args.requested_pcs, sample_count=len(fam_ids), marker_rank=marker_matrix_rank) pca = attempt / "pca" if realized: plink_outputs.append(_run_required_tool([str(tools["plink"]), "--bfile", str(bed), "--pca", str(realized), "--allow-extra-chr", "--out", str(pca)], "PLINK PCA", required_output=pca.with_suffix(".eigenvec"), cwd=attempt, evidence_sink=logs)) eigenvectors = _read_eigenvectors(pca.with_suffix(".eigenvec"), tfam_ids, realized) else: eigenvectors = np.empty((len(tfam_ids), 0)) design_rank, residual_df = _validate_design(eigenvectors, realized_pcs=realized) phenotype = attempt / "phenotype.txt" phenotype.write_text("".join(f"{sample}\t{sample}\t{value:.10g}\n" for sample, value in _ordered_phenotype(trait, tfam_ids)), encoding="utf-8") covariates = attempt / "covariates.txt" covariates.write_text("".join("\t".join([sample, sample, "1"] + [f"{value:.10g}" for value in eigenvectors[i]]) + "\n" for i, sample in enumerate(tfam_ids)), encoding="utf-8") output = attempt / "emmax_output" _run_required_tool([str(tools["emmax"]), "-v", "-d", "10", "-t", str(trans), "-p", str(phenotype), "-k", str(kinship), "-c", str(covariates), "-o", str(output)], "EMMAX", required_output=output.with_suffix(".ps"), cwd=attempt, evidence_sink=logs) full_source, p_source = attempt / "emmax_full.tsv", attempt / "emmax_pvalues.tsv" ps_metrics = _parse_emmax_ps( output.with_suffix(".ps"), tested_marker_map, full_source, p_source ) if ps_metrics["emmax_tested_count"] != estimability_metrics["emmax_tested_count"]: raise ValueError("EMMAX tested-marker accounting is inconsistent") environment += f"emmax_ps_source_schema={ps_metrics['source_ps_schema']}\n" metrics = {"genotype_sample_count": len(canonical.samples), "phenotype_sample_count": len(trait["values"]), "retained_sample_count": len(tfam_ids), "genotype_only_sample_ids": genotype_only, "input_variant_count": canonical.input_variant_count, "nonbiallelic_excluded_count": canonical.nonbiallelic_excluded_count, "plink_retained_count": len(bim_ids), "tped_marker_count": len(tped_ids), "half_call_count": canonical.half_call_count, "vcf_half_call_policy": "missing", "requested_pcs": args.requested_pcs, "realized_pcs": realized, "marker_matrix_rank": marker_matrix_rank, "marker_matrix_eigenvalues": [float(value) for value in marker_matrix_eigenvalues], "marker_matrix_minimum_eigenvalue": float(marker_matrix_eigenvalues[0]), "marker_matrix_maximum_eigenvalue": float(marker_matrix_eigenvalues[-1]), "design_rank": design_rank, "residual_degrees_of_freedom": residual_df, "relationship_rows": len(tfam_ids), "relationship_columns": len(tfam_ids), "relationship_diagonal_min": relationship.diagonal_min, "relationship_diagonal_max": relationship.diagonal_max, "relationship_minimum_eigenvalue": relationship.minimum_eigenvalue, "relationship_maximum_eigenvalue": relationship.maximum_eigenvalue, "relationship_eigenvalues": [float(value) for value in relationship.eigenvalues], "relationship_psd_tolerance": relationship.tolerance, "plink_warnings": _warning_lines(plink_outputs), **estimability_metrics, **ps_metrics} stem = f"{safe_prefix}.{STAGE}.{identity12}" destinations = {"emmax_full": final_dir / f"{stem}.full.tsv", "emmax_pvalues": final_dir / f"{stem}.pvalues.tsv", "report": final_dir / f"{stem}.report.md", "environment": final_dir / f"{stem}.environment.txt"} report_source, environment_source = attempt / "report.md", attempt / "environment.txt" report_source.write_text("\n".join(["# EMMAX association report", "", f"- Samples retained: {len(tfam_ids)}", f"- Input variants: {canonical.input_variant_count}", f"- Multiallelic excluded: {canonical.nonbiallelic_excluded_count}", f"- Candidate markers: {estimability_metrics['emmax_candidate_marker_count']}", f"- No-called genotype excluded: {estimability_metrics['no_called_genotype_excluded_count']}", f"- Nonvariable genotype excluded: {estimability_metrics['nonvariable_genotype_excluded_count']}", f"- Estimability policy: {ESTIMABILITY_POLICY}", f"- EMMAX tested markers: {ps_metrics['emmax_tested_count']}", f"- Requested/realized PCs: {args.requested_pcs}/{realized}", f"- Source .ps schema: {ps_metrics['source_ps_schema']} ({ps_metrics['source_ps_field_count']} fields)", f"- Standard errors: {ps_metrics['finite_se_count']} finite; {ps_metrics['missing_se_count']} missing (NA)", "- Official package .ps schema: SNPID_beta_p (3 fields; SE=NA)", f"- Kinship diagonal range: {relationship.diagonal_min:.10g}..{relationship.diagonal_max:.10g}", f"- Kinship eigenvalue range: {relationship.minimum_eigenvalue:.10g}..{relationship.maximum_eigenvalue:.10g}", ""]) , encoding="utf-8") environment_source.write_text(environment, encoding="utf-8") sources = {"emmax_full": full_source, "emmax_pvalues": p_source, "report": report_source, "environment": environment_source} _revalidate_input_records(inputs) ended = datetime.now(timezone.utc) manifest.update({"status": "success", "timestamp": ended.isoformat(), "end_time": ended.isoformat(), "elapsed_sec": round(time.monotonic() - clock, 6), "inputs": inputs, "software": {"plink": plink_version, "emmax_release": emmax_release, "emmax_package": emmax_provenance["package"], "emmax_sha256": emmax_hash, "R": r_version, "python": platform.python_version()}, "metrics": metrics, "artifacts": {role: _artifact(workspace, source, destinations[role], role, "text/markdown" if role == "report" else "text/plain" if role == "environment" else "text/tab-separated-values") for role, source in sources.items()}}) manifest_source = attempt / "manifest.json" _write_json(manifest_source, manifest) manifest_destination = final_dir / f"{stem}.manifest.json" bundle = _write_reproducibility(workspace, identity12, run_id, runner, manifest_source, manifest, environment, sources, manifest_destination) published = _publish_with_bundle( sources, destinations, manifest_source, manifest_destination, workspace / ".locks" / f"{stem}.lock", started.isoformat(), bundle, ) failed = final_dir / f"{safe_prefix}.{STAGE}.failed.manifest.json" if published: _clear_failed_manifest_if_winner( failed, workspace / ".locks" / f"{safe_prefix}.{STAGE}.failed.lock", manifest_destination, run_id, ) disposition = "published" if published else "superseded by a newer success" log_file.write_text("".join(logs) + f"\n[success] {STAGE} completed ({disposition})\n", encoding="utf-8") return 0 except Exception as exc: ended = datetime.now(timezone.utc) root_cause = re.sub(r"\s+", " ", str(exc)).strip()[:1000] or exc.__class__.__name__ manifest.update({"status": "failed", "timestamp": ended.isoformat(), "end_time": ended.isoformat(), "elapsed_sec": round(time.monotonic() - clock, 6), "artifacts": {}, "error": root_cause}) failure_path = final_dir / f"{safe_prefix}.{STAGE}.failed.manifest.json" success_path = None if identity12 != "failed": success_path = final_dir / f"{safe_prefix}.{STAGE}.{identity12}.manifest.json" failure_written = _publish_failed_manifest( failure_path, manifest, workspace / ".locks" / f"{safe_prefix}.{STAGE}.failed.lock", success_path, ) try: if failure_written: _write_reproducibility(workspace, identity12, run_id, runner, failure_path, manifest, environment) except Exception: pass log_file.write_text("".join(logs) + f"\n[error] {root_cause}\n", encoding="utf-8") print(f"ERROR: {root_cause}", file=sys.stderr) return 1 if __name__ == "__main__": raise SystemExit(main())