#!/usr/bin/env python3 """Sort an external VCF and convert it to TASSEL HapMap format.""" from __future__ import annotations import argparse import gzip import json import os import platform import re import shlex import shutil import subprocess import sys import time from datetime import datetime, timezone from pathlib import Path from typing import TextIO from gfa_core.genotype_context import genotype_scope_id_from_sha256, sha256_file STAGE = "gwas_fmt" SOURCE_LABEL = "18_gwas_tassel_convert" def _open_text(path: Path) -> TextIO: if path.suffix == ".gz": return gzip.open(path, "rt", encoding="utf-8") return path.open("r", encoding="utf-8") def _vcf_shape(path: Path) -> tuple[list[str], int]: samples: list[str] | None = None variants = 0 with _open_text(path) as handle: for raw in handle: if raw.startswith("#CHROM\t"): samples = raw.rstrip("\r\n").split("\t")[9:] elif raw and not raw.startswith("#"): variants += 1 if samples is None: raise ValueError(f"VCF header #CHROM was not found: {path}") return samples, variants def _vcf_genotypes(path: Path) -> tuple[list[str], dict[tuple[str, str, str, str], tuple[str, ...]]]: samples: list[str] | None = None records: dict[tuple[str, str, str, str], tuple[str, ...]] = {} with _open_text(path) as handle: for line_number, raw in enumerate(handle, start=1): if raw.startswith("#CHROM\t"): columns = raw.rstrip("\r\n").split("\t") samples = columns[9:] continue if not raw.strip() or raw.startswith("#"): continue fields = raw.rstrip("\r\n").split("\t") if len(fields) < 9: raise ValueError(f"Malformed VCF record at line {line_number}") key = (fields[0], fields[1], fields[3], fields[4]) if key in records: raise ValueError(f"Duplicate VCF variant key {key} in {path}") format_fields = fields[8].split(":") try: gt_index = format_fields.index("GT") except ValueError as exc: raise ValueError(f"VCF record has no GT field at {fields[0]}:{fields[1]}") from exc calls = [] for sample_field in fields[9:]: values = sample_field.split(":") gt = values[gt_index] if gt_index < len(values) else "." calls.append(_canonical_gt(gt)) records[key] = tuple(calls) if samples is None: raise ValueError(f"VCF header #CHROM was not found: {path}") if not samples: raise ValueError(f"VCF must contain at least one sample: {path}") if not records: raise ValueError(f"VCF must contain at least one variant: {path}") return samples, records def _canonical_gt(value: str) -> str: raw = str(value or ".").replace("|", "/") tokens = raw.split("/") if not tokens or any(token in {"", "."} for token in tokens): return "." try: return "/".join(str(index) for index in sorted(int(token) for token in tokens)) except ValueError as exc: raise ValueError(f"Invalid VCF GT value: {value}") from exc _IUPAC = { frozenset(("A", "G")): "R", frozenset(("C", "T")): "Y", frozenset(("G", "C")): "S", frozenset(("A", "T")): "W", frozenset(("G", "T")): "K", frozenset(("A", "C")): "M", } def _sample_index_map(input_samples: list[str], converted_samples: list[str]) -> list[int]: """Map TASSEL's deterministic ``sample_sample`` aliases back to VCF IDs.""" if len(input_samples) != len(converted_samples): raise RuntimeError("TASSEL conversion changed sample count") indexes: list[int] = [] used: set[int] = set() for sample in input_samples: matches = [ index for index, converted in enumerate(converted_samples) if converted == sample or converted == f"{sample}_{sample}" ] matches = [index for index in matches if index not in used] if len(matches) != 1: raise RuntimeError(f"TASSEL conversion changed or ambiguously renamed sample: {sample}") indexes.append(matches[0]) used.add(matches[0]) return indexes def _gt_to_hapmap(gt: str, ref: str, alt: str) -> str: if gt == ".": return "N" alleles = [ref, *alt.split(",")] try: bases = [alleles[int(token)].upper() for token in gt.split("/")] except (IndexError, ValueError): return "N" if any(len(base) != 1 or base not in {"A", "C", "G", "T"} for base in bases): return "N" unique = frozenset(bases) if len(unique) == 1: return next(iter(unique)) return _IUPAC.get(unique, "N") def _strict_hapmap_marker(key: tuple[str, str, str, str], genotypes: tuple[str, ...]) -> bool: """Return whether a VCF marker is reversibly representable in HapMap.""" _chrom, _pos, ref, alt = key ref = ref.upper() alt = alt.upper() if len(ref) != 1 or ref not in {"A", "C", "G", "T"}: return False if "," in alt or len(alt) != 1 or alt not in {"A", "C", "G", "T"} or alt == ref: return False for gt in genotypes: if gt == ".": continue tokens = gt.split("/") if len(tokens) != 2 or any(token not in {"0", "1"} for token in tokens): return False return True def _marker_audit(records: dict[tuple[str, str, str, str], tuple[str, ...]]) -> dict[str, object]: strict = {key for key, genotypes in records.items() if _strict_hapmap_marker(key, genotypes)} unsupported = set(records) - strict return { "strict_keys": strict, "unsupported_keys": unsupported, "strict_marker_count": len(strict), "unsupported_marker_count": len(unsupported), "lossy_conversion": bool(unsupported), } def _hapmap_genotypes(path: Path) -> tuple[list[str], dict[tuple[str, str], tuple[str, tuple[str, ...]]]]: samples, _count = _hapmap_shape(path) records: dict[tuple[str, str], tuple[str, tuple[str, ...]]] = {} with path.open("r", encoding="utf-8") as handle: next(handle) for line_number, raw in enumerate(handle, start=2): if not raw.strip(): continue fields = raw.rstrip("\r\n").split("\t") if len(fields) < 11 + len(samples): raise ValueError(f"Malformed HapMap record at line {line_number}") key = (fields[2], fields[3]) if key in records: raise ValueError(f"Duplicate HapMap locus {key} in {path}") records[key] = (fields[1], tuple(value.upper() for value in fields[11:])) return samples, records def _prove_genotype_conservation(input_vcf: Path, sorted_vcf: Path, hapmap: Path) -> dict[str, object]: input_samples, input_records = _vcf_genotypes(input_vcf) sorted_samples, sorted_records = _vcf_genotypes(sorted_vcf) sorted_indexes = _sample_index_map(input_samples, sorted_samples) normalized_sorted: dict[tuple[str, str], tuple[tuple[str, str, str, str], tuple[str, ...]]] = {} for sorted_key, calls in sorted_records.items(): locus = (sorted_key[0].casefold(), sorted_key[1]) if locus in normalized_sorted: raise RuntimeError(f"TASSEL sorted VCF produced ambiguous normalized locus {sorted_key[0]}:{sorted_key[1]}") normalized_sorted[locus] = (sorted_key, calls) if {(key[0].casefold(), key[1]) for key in input_records} != set(normalized_sorted): raise RuntimeError("TASSEL sorted VCF changed chrom,pos locus identity") marker_audit = _marker_audit(input_records) strict_keys = marker_audit["strict_keys"] for key in strict_keys: expected = input_records[key] sorted_key, sorted_calls = normalized_sorted[(key[0].casefold(), key[1])] if sorted_key[2:] != key[2:]: raise RuntimeError(f"TASSEL sorted VCF changed strict SNP allele identity at {key[0]}:{key[1]}") observed_sorted = tuple(sorted_calls[index] for index in sorted_indexes) if observed_sorted != expected: raise RuntimeError(f"TASSEL sorted VCF changed genotype calls at {key[0]}:{key[1]}") hapmap_samples, hapmap_records = _hapmap_genotypes(hapmap) hapmap_indexes = _sample_index_map(input_samples, hapmap_samples) normalized_input_loci = {(key[0].casefold(), key[1]) for key in input_records} normalized_hapmap: dict[tuple[str, str], tuple[str, str, tuple[str, ...]]] = {} for (chrom, pos), (alleles, calls) in hapmap_records.items(): normalized_key = (chrom.casefold(), pos) if normalized_key in normalized_hapmap: raise RuntimeError(f"TASSEL HapMap produced ambiguous normalized locus {chrom}:{pos}") normalized_hapmap[normalized_key] = (chrom, alleles, calls) if normalized_input_loci != set(normalized_hapmap): raise RuntimeError("TASSEL HapMap changed chrom,pos locus identity") compared = 0 normalized_contigs = 0 for (chrom, pos, ref, alt), genotypes in input_records.items(): observed_chrom, allele_text, raw_observed = normalized_hapmap[(chrom.casefold(), pos)] if observed_chrom != chrom: normalized_contigs += 1 if (chrom, pos, ref, alt) not in strict_keys: continue observed = tuple(raw_observed[index] for index in hapmap_indexes) expected_alleles = {ref.upper(), *(item.upper() for item in alt.split(","))} observed_alleles = {item.upper() for item in allele_text.replace(",", "/").split("/") if item} if observed_alleles != expected_alleles: raise RuntimeError(f"TASSEL HapMap changed allele identity at {chrom}:{pos}") expected = tuple(_gt_to_hapmap(gt, ref, alt) for gt in genotypes) if observed != expected: raise RuntimeError(f"TASSEL HapMap genotype discordance at {chrom}:{pos}") compared += len(expected) return { "variants_compared": len(strict_keys), "genotypes_compared": compared, "strict_marker_count": marker_audit["strict_marker_count"], "unsupported_marker_count": marker_audit["unsupported_marker_count"], "normalized_contig_count": normalized_contigs, "strict_genotype_conservation": True, "lossy_conversion": marker_audit["lossy_conversion"], } def _hapmap_shape(path: Path) -> tuple[list[str], int]: with path.open("r", encoding="utf-8") as handle: header = handle.readline().rstrip("\r\n").split("\t") if len(header) < 11 or header[:4] != ["rs#", "alleles", "chrom", "pos"]: raise ValueError(f"Invalid TASSEL HapMap header: {path}") samples = header[11:] variants = sum(1 for raw in handle if raw.strip()) return samples, variants def _artifact(workspace: Path, path: Path, role: str, media_type: str) -> dict[str, object]: size = path.stat().st_size return { "artifact": role, "role": role, "path": str(path.relative_to(workspace)), "full_path": str(path.resolve()), "type": media_type, "media_type": media_type, "size": size, "size_bytes": size, "sha256": sha256_file(path), } def _run(command: list[str]) -> str: try: completed = subprocess.run( command, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True, check=False, ) except OSError: raise RuntimeError("TASSEL execution could not be started") from None output = completed.stdout or "" if output: print(output, end="" if output.endswith("\n") else "\n", flush=True) if completed.returncode != 0: raise RuntimeError(f"TASSEL execution failed with exit code {completed.returncode}") return output def _public_failure_error(exc: Exception) -> str: text = str(exc).strip() malformed = re.match(r"^(Malformed (?:VCF|HapMap) record)(?: at line (\d+))?", text) if malformed: return f"{malformed.group(1)} at line {malformed.group(2)}" if malformed.group(2) else malformed.group(1) for private_prefix, public_error in ( ("Invalid VCF GT value", "Invalid VCF GT value"), ("VCF record has no GT field", "VCF record has no GT field"), ("TASSEL conversion changed or ambiguously renamed sample", "TASSEL conversion changed or ambiguously renamed sample"), ("Duplicate VCF variant key", "Duplicate VCF variant key"), ("Duplicate HapMap locus", "Duplicate HapMap locus"), ("VCF header #CHROM was not found", "VCF header #CHROM was not found"), ("VCF must contain at least one sample", "VCF must contain at least one sample"), ("VCF must contain at least one variant", "VCF must contain at least one variant"), ("Invalid TASSEL HapMap header", "Invalid TASSEL HapMap header"), ("External VCF does not exist", "External VCF does not exist"), ): if text.startswith(private_prefix): return public_error safe_exact = { "TASSEL conversion changed sample count", "TASSEL sorted VCF changed chrom,pos locus identity", "TASSEL HapMap changed chrom,pos locus identity", "Mask policy cannot prove conversion: strict_marker_count is zero", "TASSEL pipeline executable was not found or is not executable", "TASSEL execution could not be started", "TASSEL completed without producing the sorted VCF and HapMap outputs", "TASSEL version could not be proven from executable or run log", "Mask conversion is unproven: strict_marker_count and genotypes_compared must both be positive", } if text in safe_exact: return text if re.fullmatch(r"TASSEL execution failed with exit code \d+", text): return text if re.fullmatch(r"\d+ unsupported HapMap marker\(s\) found under fail policy", text): return text if re.fullmatch( r"TASSEL conversion changed sample or variant cardinality: input=\d+/\d+, sorted=\d+/\d+, hapmap=\d+/\d+", text, ): return text for private_prefix in ( "TASSEL sorted VCF produced ambiguous normalized locus", "TASSEL sorted VCF changed strict SNP allele identity", "TASSEL sorted VCF changed genotype calls", "TASSEL HapMap produced ambiguous normalized locus", "TASSEL HapMap changed allele identity", "TASSEL HapMap genotype discordance", ): if text.startswith(private_prefix): return private_prefix actionable = re.search(r"(?i)(?:line|exit code)\s+\d+", text) if actionable: return f"{exc.__class__.__name__}: {actionable.group(0)}" return "GWAS TASSEL conversion failed" def _resolved_command(command: str) -> str: if command == "run_pipeline.pl": command = str(os.environ.get("TASSEL_PIPELINE") or command) candidate = Path(command).expanduser() if candidate.is_absolute() or "/" in command: return str(candidate.resolve()) return str(Path(shutil.which(command) or command)) def _version(command: list[str]) -> str: try: completed = subprocess.run(command, text=True, capture_output=True, check=False, timeout=10) text = (completed.stdout or completed.stderr or "").strip().splitlines() return text[0][:300] if text else "unknown" except Exception: return "unknown" def _parse_tassel_version(text: str) -> str: for pattern in ( r"(?i)TASSEL(?:\s+version)?\s*[:=]?\s*(\d+\.\d+(?:\.\d+)?)", r"(?i)version\s*[:=]?\s*(5\.\d+(?:\.\d+)?)", ): match = re.search(pattern, str(text or "")) if match: return match.group(1) return "unknown" def _tassel_version(resolved_tassel: str, run_output: str = "") -> str: parsed = _parse_tassel_version(run_output) if parsed != "unknown": return parsed try: completed = subprocess.run( [resolved_tassel, "-version"], text=True, capture_output=True, check=False, timeout=10, ) return _parse_tassel_version(f"{completed.stdout or ''}\n{completed.stderr or ''}") except Exception: return "unknown" def _write_json(path: Path, data: dict[str, object]) -> None: path.write_text(json.dumps(data, ensure_ascii=False, indent=2) + "\n", encoding="utf-8") def _safe_identity_token(value: str) -> str: token = re.sub(r"[^A-Za-z0-9._-]+", "_", str(value)).strip("._-") return token or "cohort" def _write_reproducibility_bundle( *, workspace: Path, prefix: str, input_sha256: str, manifest_path: Path, manifest: dict[str, object], runner_command: list[str], environment_text: str, log_path: Path | None, ) -> Path: bundle_dir = workspace / "reproducibility" / STAGE / f"{_safe_identity_token(prefix)}.{input_sha256[:12]}" if bundle_dir.is_symlink() or bundle_dir.is_file(): bundle_dir.unlink() elif bundle_dir.is_dir(): shutil.rmtree(bundle_dir) bundle_dir.mkdir(parents=True, exist_ok=True) canonical_runner = shlex.join(runner_command) (bundle_dir / "commands.sh").write_text( "#!/usr/bin/env bash\nset -euo pipefail\n" + canonical_runner + "\n", encoding="utf-8", ) (bundle_dir / "environment.txt").write_text(environment_text, encoding="utf-8") manifest_relative = str(manifest_path.relative_to(workspace)) log_relative = _workspace_relative_or_absolute(workspace, log_path) environment_relative = str((bundle_dir / "environment.txt").relative_to(workspace)) status = str(manifest.get("status") or "failed") progress = { "stage": STAGE, "status": status, "job_id": manifest.get("job_id", ""), "prefix": prefix, "genotype_scope_id": manifest.get("genotype_scope_id", ""), "input_sha256": input_sha256, "started_at": manifest.get("start_time", ""), "completed_at": manifest.get("end_time", ""), "manifest_path": manifest_relative, "log_path": log_relative, "message": manifest.get("error") if status == "failed" else f"{STAGE} completed successfully", } provenance = { "stage": STAGE, "status": status, "job_id": manifest.get("job_id", ""), "prefix": prefix, "genotype_scope_id": manifest.get("genotype_scope_id", ""), "input_sha256": input_sha256, "runner_command": canonical_runner, "manifest_path": manifest_relative, "log_path": log_relative, "environment_path": environment_relative, "software": manifest.get("software", {}), "artifacts": manifest.get("artifacts", {}), } _write_json(bundle_dir / "progress.json", progress) _write_json(bundle_dir / "provenance.json", provenance) checksum_targets = [manifest_path] artifacts = manifest.get("artifacts", {}) if isinstance(artifacts, dict): for artifact in artifacts.values(): if not isinstance(artifact, dict) or "path" not in artifact: continue candidate = workspace / str(artifact["path"]) if candidate.is_file(): checksum_targets.append(candidate) checksum_lines = [ f"{sha256_file(path)} {path.relative_to(workspace)}" for path in sorted(checksum_targets) ] (bundle_dir / "checksums.sha256").write_text("\n".join(checksum_lines) + "\n", encoding="utf-8") return bundle_dir def _workspace_relative_or_absolute(workspace: Path, path: Path | None) -> str: if path is None: return "" resolved = path.expanduser().resolve() try: return str(resolved.relative_to(workspace)) except ValueError: return str(resolved) def _runner_command(args: argparse.Namespace, workspace: Path, input_vcf: Path) -> list[str]: return [ "python", "-m", "apps.cli.run_gwas_tassel", "--workspace", str(workspace), "--job-id", args.job_id, "--prefix", args.prefix, "--vcf", str(input_vcf), "--xmx", args.xmx, "--unsupported-marker-policy", args.unsupported_marker_policy, ] def _environment_text( *, resolved_tassel: str, tassel_version: str, resolved_java: str, java_version: str, ) -> str: return "\n".join([ "runner=apps.cli.run_gwas_tassel.py", f"runtime_environment={os.environ.get('CONDA_DEFAULT_ENV', 'unknown')}", f"python_resolved={Path(sys.executable).resolve()}", f"python_version={platform.python_version()}", "external_tool=TASSEL", f"external_launcher={resolved_tassel}", f"tassel_version={tassel_version}", f"java_resolved={resolved_java}", f"java_version={java_version}", f"platform={platform.platform()}", ]) + "\n" 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", "-p", required=True) parser.add_argument("--vcf", required=True, type=Path) parser.add_argument("--tassel-pipeline", default="run_pipeline.pl") parser.add_argument("--xmx", default="10g") parser.add_argument("--unsupported-marker-policy", choices=("mask", "fail"), default="mask") parser.add_argument("--log-file", type=Path) return parser def main(argv: list[str] | None = None) -> int: args = build_parser().parse_args(argv) workspace = args.workspace.resolve() input_vcf = args.vcf.resolve() final_dir = workspace / "final" stage_dir = workspace / STAGE final_dir.mkdir(parents=True, exist_ok=True) stage_dir.mkdir(parents=True, exist_ok=True) manifest_path = final_dir / f"{args.prefix}.{STAGE}.manifest.json" sorted_work = stage_dir / "all.sorted.vcf.gz" hapmap_base = stage_dir / "hapmap" hapmap_work = stage_dir / "hapmap.hmp.txt" sorted_final = final_dir / f"{args.prefix}.{STAGE}.sorted.vcf.gz" hapmap_final = final_dir / f"{args.prefix}.{STAGE}.hapmap.hmp.txt" environment_path = final_dir / f"{args.prefix}.{STAGE}.environment.txt" report_path = final_dir / f"{args.prefix}.{STAGE}.report.md" formal_paths = (sorted_final, hapmap_final, environment_path, report_path) for formal_path in formal_paths: formal_path.unlink(missing_ok=True) (final_dir / f"{args.prefix}.{STAGE}.commands.sh").unlink(missing_ok=True) (final_dir / f"{args.prefix}.{STAGE}.log.txt").unlink(missing_ok=True) digest = "" genotype_scope_id = "" started_at = datetime.now(timezone.utc) started_clock = time.monotonic() runner_command = _runner_command(args, workspace, input_vcf) environment_text = _environment_text( resolved_tassel="unknown", tassel_version="unknown", resolved_java="unknown", java_version="unknown", ) manifest: dict[str, object] = { "manifest_version": "1.0", "stage": STAGE, "prefix": args.prefix, "job_id": args.job_id, "run_id": f"{STAGE}_{started_at.strftime('%Y%m%d_%H%M%S')}", "status": "failed", "timestamp": started_at.isoformat(), "start_time": started_at.isoformat(), "end_time": started_at.isoformat(), "elapsed_sec": 0.0, "cmd": runner_command, "tool_versions": {"tassel": "unknown", "java": "unknown", "python": platform.python_version()}, "workspace_root": str(workspace), "workspace": str(final_dir), "input_dir": str(input_vcf.parent), "parameters": { "source_label": SOURCE_LABEL, "tassel_pipeline": args.tassel_pipeline, "xmx": args.xmx, "unsupported_marker_policy": args.unsupported_marker_policy, }, "provenance": { "runner_command": shlex.join(runner_command), "log_path": _workspace_relative_or_absolute(workspace, args.log_file), }, "metrics": {}, "artifacts": {}, } try: if not input_vcf.is_file(): raise FileNotFoundError(f"External VCF does not exist: {input_vcf}") digest = sha256_file(input_vcf) genotype_scope_id = genotype_scope_id_from_sha256(digest) manifest["genotype_scope_id"] = genotype_scope_id manifest["input_sha256"] = digest manifest["inputs"] = [{ "role": "external_vcf", "path": str(input_vcf), "size_bytes": input_vcf.stat().st_size, "sha256": digest, }] # Reject semantically invalid or empty cohorts before invoking TASSEL. _input_samples_for_audit, input_records_for_audit = _vcf_genotypes(input_vcf) marker_audit = _marker_audit(input_records_for_audit) manifest["metrics"] = { "strict_marker_count": marker_audit["strict_marker_count"], "unsupported_marker_count": marker_audit["unsupported_marker_count"], "lossy_conversion": marker_audit["lossy_conversion"], } warnings: list[str] = [] if marker_audit["unsupported_marker_count"]: warnings.append( f"{marker_audit['unsupported_marker_count']} marker(s) are not reversibly representable in HapMap; " "strict genotype conservation applies only to the reversible diploid biallelic SNP subset." ) manifest["warnings"] = warnings if args.unsupported_marker_policy == "fail" and marker_audit["unsupported_marker_count"]: raise ValueError( f"{marker_audit['unsupported_marker_count']} unsupported HapMap marker(s) found under fail policy" ) if args.unsupported_marker_policy == "mask" and int(marker_audit["strict_marker_count"]) <= 0: raise RuntimeError("Mask policy cannot prove conversion: strict_marker_count is zero") resolved_tassel = _resolved_command(args.tassel_pipeline) if not Path(resolved_tassel).is_file() or not os.access(resolved_tassel, os.X_OK): raise FileNotFoundError("TASSEL pipeline executable was not found or is not executable") resolved_java = _resolved_command("java") software = { "tassel": {"requested": args.tassel_pipeline, "resolved_path": resolved_tassel, "version": "unknown"}, "java": {"resolved_path": resolved_java, "version": _version([resolved_java, "-version"])}, "python": {"resolved_path": sys.executable, "version": platform.python_version()}, } manifest["software"] = software manifest["tool_versions"] = { "tassel": "unknown", "java": str(software["java"]["version"]), "python": platform.python_version(), } environment_text = _environment_text( resolved_tassel=resolved_tassel, tassel_version="unknown", resolved_java=resolved_java, java_version=str(software["java"]["version"]), ) sort_command = [ args.tassel_pipeline, f"-Xmx{args.xmx}", "-SortGenotypeFilePlugin", "-inputFile", str(input_vcf), "-outputFile", str(sorted_work), "-fileType", "VCF", ] export_command = [ args.tassel_pipeline, f"-Xmx{args.xmx}", "-importGuess", str(sorted_work), "-export", str(hapmap_base), "-exportType", "Hapmap", ] sorted_work.unlink(missing_ok=True) hapmap_work.unlink(missing_ok=True) tassel_output = _run(sort_command) + _run(export_command) if not sorted_work.is_file() or not hapmap_work.is_file(): raise RuntimeError("TASSEL completed without producing the sorted VCF and HapMap outputs") input_samples, input_variants = _vcf_shape(input_vcf) sorted_samples, sorted_variants = _vcf_shape(sorted_work) hapmap_samples, hapmap_variants = _hapmap_shape(hapmap_work) _sample_index_map(input_samples, sorted_samples) _sample_index_map(input_samples, hapmap_samples) sample_conservation = len(input_samples) == len(sorted_samples) == len(hapmap_samples) variant_conservation = input_variants == sorted_variants == hapmap_variants if not sample_conservation or not variant_conservation: raise RuntimeError( "TASSEL conversion changed sample or variant cardinality: " f"input={len(input_samples)}/{input_variants}, " f"sorted={len(sorted_samples)}/{sorted_variants}, " f"hapmap={len(hapmap_samples)}/{hapmap_variants}" ) software["tassel"]["version"] = _tassel_version(resolved_tassel, tassel_output) if software["tassel"]["version"] == "unknown": raise RuntimeError("TASSEL version could not be proven from executable or run log") manifest["tool_versions"]["tassel"] = str(software["tassel"]["version"]) environment_text = _environment_text( resolved_tassel=resolved_tassel, tassel_version=str(software["tassel"]["version"]), resolved_java=resolved_java, java_version=str(software["java"]["version"]), ) environment_path.write_text(environment_text, encoding="utf-8") genotype_metrics = _prove_genotype_conservation(input_vcf, sorted_work, hapmap_work) if int(genotype_metrics["strict_marker_count"]) <= 0 or int(genotype_metrics["genotypes_compared"]) <= 0: raise RuntimeError( "Mask conversion is unproven: strict_marker_count and genotypes_compared must both be positive" ) shutil.copy2(sorted_work, sorted_final) shutil.copy2(hapmap_work, hapmap_final) report_path.write_text( f"# GWAS format conversion\n\n" f"- Input SHA256: `{digest}`\n" f"- Genotype scope: `{genotype_scope_id}`\n" f"- Samples conserved: {len(input_samples)}\n" f"- Variants conserved: {input_variants}\n" f"- Unsupported marker policy: `{args.unsupported_marker_policy}`\n" f"- Strict reversible marker subset: {genotype_metrics['strict_marker_count']}\n" f"- Unsupported/lossy markers: {genotype_metrics['unsupported_marker_count']}\n" f"- Strict genotype conservation: {genotype_metrics['strict_genotype_conservation']}\n", encoding="utf-8", ) manifest["metrics"] = { "input": {"samples": len(input_samples), "variants": input_variants}, "sorted_vcf": {"samples": len(sorted_samples), "variants": sorted_variants}, "hapmap": {"samples": len(hapmap_samples), "variants": hapmap_variants}, "sample_conservation": sample_conservation, "variant_conservation": variant_conservation, "genotype_conservation": not bool(genotype_metrics["lossy_conversion"]), **genotype_metrics, } manifest["artifacts"] = { "sorted_vcf": _artifact(workspace, sorted_final, "sorted_vcf", "application/gzip"), "tassel_hapmap": _artifact(workspace, hapmap_final, "tassel_hapmap", "text/tab-separated-values"), "environment": _artifact(workspace, environment_path, "environment", "text/plain"), "report": _artifact(workspace, report_path, "report", "text/markdown"), } manifest["status"] = "success" completed_at = datetime.now(timezone.utc).isoformat() manifest["timestamp"] = completed_at manifest["end_time"] = completed_at manifest["elapsed_sec"] = round(time.monotonic() - started_clock, 6) _write_json(manifest_path, manifest) _write_reproducibility_bundle( workspace=workspace, prefix=args.prefix, input_sha256=digest, manifest_path=manifest_path, manifest=manifest, runner_command=runner_command, environment_text=environment_text, log_path=args.log_file, ) return 0 except Exception as exc: reproducibility_publication_failed = manifest.get("status") == "success" failed_software = manifest.get("software") if isinstance(failed_software, dict) and isinstance(failed_software.get("tassel"), dict): failed_tassel = failed_software["tassel"] resolved_failed_tassel = str(failed_tassel.get("resolved_path") or "") if resolved_failed_tassel: failed_tassel["version"] = _tassel_version(resolved_failed_tassel) tool_versions = manifest.get("tool_versions") if isinstance(tool_versions, dict): tool_versions["tassel"] = str(failed_tassel["version"]) failed_java = failed_software.get("java", {}) environment_text = _environment_text( resolved_tassel=resolved_failed_tassel, tassel_version=str(failed_tassel["version"]), resolved_java=str(failed_java.get("resolved_path", "unknown")), java_version=str(failed_java.get("version", "unknown")), ) generic_reproducibility_error = "Stage-local reproducibility evidence could not be written" manifest["status"] = "failed" manifest["artifacts"] = {} manifest["error"] = generic_reproducibility_error if reproducibility_publication_failed else _public_failure_error(exc) if reproducibility_publication_failed: manifest["reproducibility_error"] = generic_reproducibility_error completed_at = datetime.now(timezone.utc).isoformat() manifest["timestamp"] = completed_at manifest["end_time"] = completed_at manifest["elapsed_sec"] = round(time.monotonic() - started_clock, 6) if genotype_scope_id: manifest["genotype_scope_id"] = genotype_scope_id for formal_path in formal_paths: formal_path.unlink(missing_ok=True) _write_json(manifest_path, manifest) if digest: try: _write_reproducibility_bundle( workspace=workspace, prefix=args.prefix, input_sha256=digest, manifest_path=manifest_path, manifest=manifest, runner_command=runner_command, environment_text=environment_text, log_path=args.log_file, ) except Exception: manifest["reproducibility_error"] = generic_reproducibility_error _write_json(manifest_path, manifest) print(f"ERROR: {manifest['error']}", file=sys.stderr) return 1 if __name__ == "__main__": raise SystemExit(main())