"""``simurg`` command-line interface (spec §2, plan §14). Available: ``version``, ``inspect`` (P1), ``ingest`` (P2), ``normalize`` (P3), ``genotype`` (P4). Later phases add the other commands of CLI-001. """ from __future__ import annotations import json import os import sys from collections.abc import Callable from pathlib import Path from typing import Annotated, Any, TypeVar import numpy as np import typer from simurg import PRODUCT_NAME, __version__ from simurg.core.errors import ExitCode, SimurgError from simurg.core.log import Reporter from simurg.core.provenance import git_revision app = typer.Typer( name="simurg", help=f"{PRODUCT_NAME}: analysis toolkit for Illumina Infinium genotyping BeadChip data.", no_args_is_help=True, add_completion=False, pretty_exceptions_enable=False, context_settings={"help_option_names": ["-h", "--help"]}, ) F = TypeVar("F", bound=Callable[..., Any]) def version_string() -> str: rev = git_revision() return f"{PRODUCT_NAME} {__version__}" + (f"+{rev}" if rev else "") def _version_cb(value: bool) -> None: if value: typer.echo(version_string()) raise typer.Exit(0) @app.callback() def main( version: Annotated[ bool, typer.Option("--version", callback=_version_cb, is_eager=True, help="Print the version and exit.") ] = False, ) -> None: """Top-level options.""" def run_guarded(reporter: Reporter, fn: Callable[[], int]) -> None: """Run a command body mapping errors to catalogue codes and exit codes (CLI-014).""" reporter.open() code = ExitCode.OK try: code = ExitCode(fn()) sys.stdout.flush() # surface a closed pipe here, not at interpreter exit except BrokenPipeError: # the consumer of stdout stopped reading (`| head`): silence further writes, exit like a SIGPIPE'd tool os.dup2(os.open(os.devnull, os.O_WRONLY), sys.stdout.fileno()) code = ExitCode.SIGPIPE except SimurgError as e: reporter.error(e.code, e.message + (f"\n hint: {e.hint}" if e.hint else ""), exc=e, file=e.path) code = e.exit_code except KeyboardInterrupt: reporter.error("E_INTERNAL", "interrupted") code = ExitCode.SIGINT except Exception as e: reporter.error( "E_INTERNAL", f"internal error: {type(e).__name__}: {e}. Rerun with --debug and report it.", exc=e ) code = ExitCode.INTERNAL finally: reporter.close() raise typer.Exit(int(code)) # -------------------------------------------------------------------------------------------------------- @app.command("version") def version_cmd() -> None: """Print the version.""" typer.echo(version_string()) # -------------------------------------------------------------------------------------------------------- def _jsonable(o: Any) -> Any: if isinstance(o, np.generic): return o.item() if isinstance(o, np.ndarray): return o.tolist() if isinstance(o, Path): return str(o) if isinstance(o, bytes): return o.hex() if isinstance(o, float) and o != o: # NaN return None return str(o) def inspect_file(path: Path) -> dict[str, Any]: from simurg.formats.detect import detect_kind kind = detect_kind(path) info: dict[str, Any] = {"file": str(path), "kind": kind} if kind == "idat": from simurg.formats.idat import read_idat md = read_idat(path).metadata() md.pop("path", None) info.update(md) elif kind == "bpm": from simurg.formats.bpm import read_bpm info.update(read_bpm(path).summary()) elif kind == "csv_manifest": from simurg.formats.csv_manifest import assay_format, read_csv_manifest m = read_csv_manifest(path, load_sequences=False) s = m.summary() s.pop("normalization_bins", None) s.pop("infinium_i_at", None) s.pop("infinium_i_gc", None) s["heading"] = m.heading s["assay_format"] = assay_format(m) info.update(s) elif kind == "egt": from simurg.formats.egt import read_egt info.update(read_egt(path).summary()) elif kind == "gtc": from simurg.formats.gtc import read_gtc, validate_gtc g = read_gtc(path, check_write_complete=False) info.update(g.summary()) info["problems"] = validate_gtc(g) elif kind == "samplesheet": from simurg.formats.samplesheet import parse_samplesheet sh = parse_samplesheet(path) info.update( { "rows": len(sh.rows), "columns": sh.columns, "key_columns": [sh.barcode_column, sh.position_column], "root_folder": str(sh.root_folder) if sh.root_folder else None, "header": sh.header, "has_sample_name": sh.has_sample_name, "warnings": [f"{c}: {m}" for c, m in sh.warnings], } ) else: raise SimurgError("E_FILE_UNREADABLE", f"{path}: unrecognised file type", path=str(path)) info.pop("path", None) return info def _print_human(info: dict[str, Any]) -> None: typer.echo(f"== {info['file']} ({info['kind']})") for k, v in info.items(): if k in ("file", "kind"): continue if isinstance(v, list) and len(v) > 12: v = [*v[:12], f"… ({len(v)} total)"] if isinstance(v, (dict, list)): v = json.dumps(v, default=_jsonable, ensure_ascii=False) typer.echo(f" {k}: {v}") @app.command("inspect") def inspect_cmd( files: Annotated[list[Path], typer.Argument(help="IDAT, BPM, CSV manifest, EGT, GTC or sample-sheet files.")], as_json: Annotated[bool, typer.Option("--json", help="Machine-readable JSON output.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, log_level: Annotated[str, typer.Option("--log-level", help="debug|info|warning|error")] = "info", quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, ) -> None: """Print the header / table of contents of array files (spec CLI-020).""" rep = Reporter( command="inspect", level="debug" if debug else log_level, json_stderr=log_json, debug=debug, quiet=quiet ) def body() -> int: results = [] failures = 0 for f in files: if not f.exists(): raise SimurgError("E_FILE_NOT_FOUND", f"{f}: no such file", path=str(f)) try: info = inspect_file(f) except SimurgError as e: rep.error(e.code, e.message, file=str(f)) failures += 1 continue results.append(info) if not as_json: _print_human(info) if as_json: out = results[0] if len(results) == 1 else results sys.stdout.write(json.dumps(out, default=_jsonable, ensure_ascii=False, indent=2) + "\n") if failures and not results: return ExitCode.USAGE return ExitCode.PARTIAL if failures else ExitCode.OK run_guarded(rep, body) # -------------------------------------------------------------------------------------------------------- @app.command("ingest") def ingest_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest, checked against the BPM.")] = None, threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, resume: Annotated[bool, typer.Option("--resume", help="Keep outputs whose inputs are unchanged.")] = False, layout: Annotated[str, typer.Option("--layout", help="nested|flat")] = "nested", on_duplicate_idat: Annotated[str, typer.Option("--on-duplicate-idat", help="error|first|newest")] = "error", log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, log_level: Annotated[str, typer.Option("--log-level", help="debug|info|warning|error")] = "info", quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, ) -> None: """IDAT pairs + BPM → canonical per-sample raw intensities (Parquet). Phase P2.""" from simurg.pipeline.ingest import IngestOptions, run_ingest if layout not in ("nested", "flat"): raise typer.BadParameter("--layout must be nested or flat") if on_duplicate_idat not in ("error", "first", "newest"): raise typer.BadParameter("--on-duplicate-idat must be error, first or newest") rep = Reporter( command="ingest", level="debug" if debug else log_level, json_stderr=log_json, debug=debug, quiet=quiet, log_file=out / "logs" / "ingest.log", ) opts = IngestOptions( out=out, bpm=bpm, idats=idats, sheet=sheet, csv=csv, threads=threads, force=force, resume=resume, layout=layout, on_duplicate_idat=on_duplicate_idat, argv=sys.argv[1:], ) run_guarded(rep, lambda: run_ingest(opts, rep)) # -------------------------------------------------------------------------------------------------------- @app.command("normalize") def normalize_cmd( indir: Annotated[Path, typer.Option("--in", help="Output folder of `simurg ingest`.")], norm_version: Annotated[str, typer.Option("--norm-version", help="1.2.0 (current) or 1.1.2 (older)")] = "1.2.0", threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, ) -> None: """Fit per-bin normalization transforms and write normalized X/Y, theta and R (phase P3).""" from simurg.pipeline.normalize import NormalizeOptions, run_normalize rep = Reporter( command="normalize", level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet ) opts = NormalizeOptions(indir=indir, version=norm_version, threads=threads, force=force, argv=sys.argv[1:]) run_guarded(rep, lambda: run_normalize(opts, rep)) # -------------------------------------------------------------------------------------------------------- @app.command("genotype") def genotype_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], egt: Annotated[Path, typer.Option("--egt", help="EGT cluster file (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest, checked against the BPM.")] = None, cutoff: Annotated[float, typer.Option("--cutoff", help="GenCall score below which a call is a no-call.")] = 0.15, preset: Annotated[ str, typer.Option("--preset", help="current (default) or autocall-1.6.2 (legacy vendor behaviour)") ] = "current", norm_version: Annotated[ str | None, typer.Option("--norm-version", help="1.2.0 or 1.1.2 (default: from --preset)") ] = None, sex_version: Annotated[ int | None, typer.Option("--sex-version", min=1, max=2, help="Sex estimate rule 1 or 2 (default: preset)") ] = None, gtc_version: Annotated[ int | None, typer.Option("--gtc-version", help="GTC format 3 or 5 (default: preset)") ] = None, autocall_date: Annotated[ str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).") ] = None, threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, controls: Annotated[ bool, typer.Option("--controls/--no-controls", help="Write controls_raw.tsv and controls_qc.tsv.") ] = True, cohort_qc: Annotated[ bool, typer.Option("--cohort-qc/--no-cohort-qc", help="Contamination and relatedness after genotyping.") ] = True, qc_config: Annotated[Path | None, typer.Option("--qc-config", help="QC threshold YAML (QC-020).")] = None, qc_profile: Annotated[str | None, typer.Option("--qc-profile", help="QC profile name or 'auto'.")] = None, threshold_call_rate: Annotated[float | None, typer.Option("--threshold-call-rate")] = None, threshold_lrr_sd: Annotated[float | None, typer.Option("--threshold-lrr-sd")] = None, threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None, bg_offset: Annotated[float | None, typer.Option("--bg-offset", help="Control background offset (3000).")] = None, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, layout: Annotated[str, typer.Option("--layout", help="nested|flat")] = "nested", on_duplicate_idat: Annotated[str, typer.Option("--on-duplicate-idat", help="error|first|newest")] = "error", log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, ) -> None: """IDAT pairs → GTC files and a sample summary (GenCall-equivalent calling; phase P4).""" from simurg.pipeline.genotype import GenotypeOptions, run_genotype if layout not in ("nested", "flat"): raise typer.BadParameter("--layout must be nested or flat") from simurg.genotype.sample import PRESETS if preset not in PRESETS: raise typer.BadParameter(f"--preset must be one of {sorted(PRESETS)}") if gtc_version not in (None, 3, 5): raise typer.BadParameter("--gtc-version must be 3 or 5") if norm_version not in (None, "1.2.0", "1.1.2"): raise typer.BadParameter("--norm-version must be 1.2.0 or 1.1.2") rep = Reporter( command="genotype", level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet, log_file=out / "logs" / "genotype.log", ) opts = GenotypeOptions( out=out, bpm=bpm, egt=egt, idats=idats, sheet=sheet, csv=csv, threads=threads, force=force, layout=layout, on_duplicate_idat=on_duplicate_idat, cutoff=cutoff, preset=preset, norm_version=norm_version, sex_version=sex_version, gtc_version=gtc_version, autocall_date=autocall_date, controls=controls, cohort_qc=cohort_qc, qc_config=qc_config, qc_profile=qc_profile, threshold_call_rate=threshold_call_rate, threshold_lrr_sd=threshold_lrr_sd, thresholds=list(threshold or []), bg_offset=bg_offset, argv=sys.argv[1:], ) run_guarded(rep, lambda: run_genotype(opts, rep)) # -------------------------------------------------------------------------------------------------------- @app.command("gtc-to-vcf") def gtc_to_vcf_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], fasta: Annotated[Path, typer.Option("--fasta", help="Reference FASTA (plain or bgzip; required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (SourceSeq; required for indels).")] = None, no_bgzip: Annotated[bool, typer.Option("--no-bgzip", help="Plain VCF, no index.")] = False, genome_cache: Annotated[bool, typer.Option("--genome-cache/--no-genome-cache", help="mmap genome cache.")] = True, exclude_loci: Annotated[ Path | None, typer.Option("--exclude-loci", "--filter-loci", help="Probe names to drop (one per line).") ] = None, include_loci: Annotated[Path | None, typer.Option("--include-loci", help="Only these probe names.")] = None, unsquash_duplicates: Annotated[bool, typer.Option("--unsquash-duplicates", help="One record per probe.")] = False, inf1_nocall_informative: Annotated[ bool, typer.Option("--inf1-nocall-informative", help="Read Infinium I no-calls as allele absence.") ] = False, drop_ambiguous: Annotated[bool, typer.Option("--drop-ambiguous", help="Drop low-confidence strand loci.")] = False, filter_gs: Annotated[float | None, typer.Option("--filter-gs", help="FILTER LowGS when GS < X.")] = None, contig_alias: Annotated[str, typer.Option("--contig-alias", help="strict|auto|FILE")] = "strict", ploidy_policy: Annotated[str, typer.Option("--ploidy-policy", help="sex-aware|diploid")] = "sex-aware", threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, merge: Annotated[ bool, typer.Option("--merge", help="One multi-sample BCF (+ .csi) instead of per-sample VCFs.") ] = False, auxiliary_loci: Annotated[ Path | None, typer.Option("--auxiliary-loci", help="VCF of true MNV definitions for manifest probes (OUT-035).") ] = None, merge_snv_indel: Annotated[ bool, typer.Option("--merge-snv-indel", help="One record for an SNV and an indel at the same position.") ] = False, info_fields: Annotated[ str, typer.Option( "--info-fields", help="Comma list: ASSAY_TYPE,NORM_BIN,GENTRAIN,CLUSTER_SEP,ILMN_STRAND,PROBE_IDS" ), ] = "", extra_format: Annotated[ str, typer.Option("--extra-format", help="Comma list: X,Y,NORMX,NORMY,THETA,R (single-probe records)") ] = "", egt: Annotated[Path | None, typer.Option("--egt", help="Cluster file (INFO GENTRAIN, CLUSTER_SEP).")] = None, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, ) -> None: """GTC files → one bgzipped, indexed SNV VCF per sample (plus-strand alleles, merged multi-probe sites).""" from simurg.pipeline.vcf import VcfOptions, run_gtc_to_vcf rep = Reporter( command="gtc-to-vcf", level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet, log_file=out / "logs" / "gtc-to-vcf.log", ) opts = VcfOptions( out=out, bpm=bpm, fasta=fasta, gtcs=gtcs, sheet=sheet, csv=csv, bgzip=not no_bgzip, genome_cache=genome_cache, exclude_loci=exclude_loci, include_loci=include_loci, unsquash_duplicates=unsquash_duplicates, inf1_nocall_informative=inf1_nocall_informative, drop_ambiguous=drop_ambiguous, filter_gs=filter_gs, contig_alias=contig_alias, ploidy_policy=ploidy_policy, threads=threads, merge=merge, merge_snv_indel=merge_snv_indel, auxiliary_loci=auxiliary_loci, force=force, info_fields=[x.strip().upper() for x in info_fields.split(",") if x.strip()], extra_format=[x.strip().upper() for x in extra_format.split(",") if x.strip()], egt=egt, argv=sys.argv[1:], ) # fmt: skip run_guarded(rep, lambda: run_gtc_to_vcf(opts, rep)) # -------------------------------------------------------------------------------------------------------- @app.command("compare") def compare_cmd( a: Annotated[Path, typer.Argument(help="GTC/VCF file or folder (A).")], b: Annotated[Path, typer.Argument(help="GTC/VCF file or folder (B), paired with A by relative path.")], out: Annotated[Path, typer.Option("--out", help="Output folder (summary.json, discrepancies, dashboard).")], bpm: Annotated[Path | None, typer.Option("--bpm", help="Manifest, for per-locus-class counts of GTCs.")] = None, tolerance: Annotated[ list[str] | None, typer.Option("--tolerance", help="FIELD=exact|abs=X|ulp=N (repeatable); VCF floats: float=abs=X"), ] = None, sample: Annotated[ list[str] | None, typer.Option("--sample", help="VCF sample(s) to compare (default: all shared).") ] = None, title: Annotated[str, typer.Option("--title", help="Dashboard title.")] = "", threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, fail_on_discrepancy: Annotated[ bool, typer.Option("--fail-on-discrepancy", help="Exit 3 when any DISCREPANCY is found.") ] = False, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, debug: Annotated[bool, typer.Option("--debug")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, ) -> None: """Field-level differential test of two GTCs/VCFs (or folders): EXACT/NEAR/DISCREPANCY/UNKNOWN (plan §16.3).""" from simurg.pipeline.compare import CompareOptions, run_compare rep = Reporter( command="compare", level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet ) def body() -> int: tol: dict[str, str] = {} for item in tolerance or []: key, sep, val = item.partition("=") if not sep: raise SimurgError("E_CONFIG_KEY", f"--tolerance expects FIELD=VALUE, got {item!r}") tol[key.strip()] = val.strip() opts = CompareOptions( a=a, b=b, out=out, bpm=bpm, tolerances=tol, samples=list(sample) if sample else None, threads=threads, fail_on_discrepancy=fail_on_discrepancy, title=title, force=force, ) return run_compare(opts, rep) run_guarded(rep, body) # -------------------------------------------------------------------------------------------------------- export_app = typer.Typer(help="Exports from GTCs or VCFs: bedgraph, finalreport, plink, locus-summary.") app.add_typer(export_app, name="export") def _plain_reporter(command: str, log_json: bool, debug: bool, quiet: bool) -> Reporter: return Reporter(command=command, level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet) @export_app.command("bedgraph") def export_bedgraph_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, no_gzip: Annotated[bool, typer.Option("--no-gzip", help="Plain .bedgraph files.")] = False, chrom_style: Annotated[str, typer.Option("--chrom-style", help="manifest|ucsc")] = "manifest", smoothed_lrr: Annotated[int, typer.Option("--smoothed-lrr", min=0, help="Extra moving-median LRR track, K.")] = 0, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, debug: Annotated[bool, typer.Option("--debug")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, ) -> None: """BAF and LRR bedGraph tracks per sample (OUT-050 … OUT-052).""" from simurg.pipeline.export import BedgraphOptions, run_export_bedgraph rep = _plain_reporter("export bedgraph", log_json, debug, quiet) opts = BedgraphOptions(out=out, bpm=bpm, gtcs=gtcs, sheet=sheet, gzip=not no_gzip, chrom_style=chrom_style, smoothed_lrr=smoothed_lrr, force=force, argv=sys.argv[1:]) # fmt: skip run_guarded(rep, lambda: run_export_bedgraph(opts, rep)) @export_app.command("finalreport") def export_finalreport_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (source/reference strands).")] = None, egt: Annotated[Path | None, typer.Option("--egt", help="Cluster file (GT Score column).")] = None, matrix: Annotated[bool, typer.Option("--matrix", help="Loci × samples genotype matrix.")] = False, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, debug: Annotated[bool, typer.Option("--debug")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, ) -> None: """Long-format final report (one row per locus × sample), or a genotype matrix (OUT-071).""" from simurg.pipeline.export import FinalReportOptions, run_export_finalreport rep = _plain_reporter("export finalreport", log_json, debug, quiet) opts = FinalReportOptions(out=out, bpm=bpm, gtcs=gtcs, sheet=sheet, csv=csv, egt=egt, matrix=matrix, force=force, argv=sys.argv[1:]) # fmt: skip run_guarded(rep, lambda: run_export_finalreport(opts, rep)) @export_app.command("plink") def export_plink_cmd( vcfs: Annotated[Path, typer.Option("--vcfs", help="Folder with SimurgArray per-sample *.snv.vcf.gz.")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (reported sex).")] = None, gtcs: Annotated[Path | None, typer.Option("--gtcs", help="GTC folder (estimated sex).")] = None, multiallelic: Annotated[str, typer.Option("--multiallelic", help="split|drop")] = "split", fmt: Annotated[str, typer.Option("--format", help="bed (PLINK 1), pgen (PLINK 2) or both")] = "bed", prefix: Annotated[str, typer.Option("--prefix", help="Output file prefix.")] = "cohort", force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, debug: Annotated[bool, typer.Option("--debug")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, ) -> None: """PLINK 1 .bed/.bim/.fam from SimurgArray VCFs (OUT-043).""" from simurg.pipeline.export import PlinkOptions, run_export_plink rep = _plain_reporter("export plink", log_json, debug, quiet) opts = PlinkOptions( out=out, vcfs=vcfs, sheet=sheet, gtcs=gtcs, multiallelic=multiallelic, format=fmt, prefix=prefix, force=force, argv=sys.argv[1:], ) run_guarded(rep, lambda: run_export_plink(opts, rep)) @export_app.command("locus-summary") def export_locus_summary_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], egt: Annotated[Path, typer.Option("--egt", help="EGT cluster file (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, threads: Annotated[int, typer.Option("--threads", min=0, max=1024)] = 0, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: Annotated[bool, typer.Option("--log-json", "--json-log")] = False, debug: Annotated[bool, typer.Option("--debug")] = False, quiet: Annotated[bool, typer.Option("--quiet", "--no-progress")] = False, ) -> None: """Per-locus summary statistics (OUT-072); same as ``simurg qc loci``.""" _qc_cohort("loci", out, bpm, gtcs, sheet, egt, threads, force, log_json, debug, quiet) # -------------------------------------------------------------------------------------------------------- qc_app = typer.Typer(help="Quality control: controls, samples, contamination, relatedness, loci, thresholds.") app.add_typer(qc_app, name="qc") def _reporter(command: str, out: Path | None, log_json: bool, debug: bool, quiet: bool) -> Reporter: return Reporter( command=command, level="debug" if debug else "info", json_stderr=log_json, debug=debug, quiet=quiet, log_file=(out / "logs" / f"{command.replace(' ', '_')}.log") if out is not None else None, ) LogJson = Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] Debug = Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] Quiet = Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] @qc_app.command("controls") def qc_controls_cmd( out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], array_type: Annotated[str, typer.Option("--array-type", help="genotyping|methylation (required).")], bpm: Annotated[Path | None, typer.Option("--bpm", help="BPM manifest (control definitions).")] = None, csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (control definitions).")] = None, idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, bg_offset: Annotated[float, typer.Option("--bg-offset", help="Background offset x (QC-003).")] = 3000.0, compat: Annotated[bool, typer.Option("--compat-column-names", help="Reference-style CSV names.")] = False, force: Annotated[bool, typer.Option("--force", help="Replace existing control tables.")] = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Control-probe QC from IDATs without genotyping (QC-012).""" from simurg.pipeline.qc import QcControlsOptions, run_qc_controls if array_type not in ("genotyping", "methylation"): raise typer.BadParameter("--array-type must be genotyping or methylation") rep = _reporter("qc controls", out, log_json, debug, quiet) opts = QcControlsOptions( out=out, bpm=bpm, csv=csv, idats=idats, sheet=sheet, array_type=array_type, bg_offset=bg_offset, compat=compat, force=force, argv=sys.argv[1:], ) run_guarded(rep, lambda: run_qc_controls(opts, rep)) # -------------------------------------------------------------------------------------------------------- cn_app = typer.Typer(help="Genome-wide copy number, LOH and mosaicism (CYTO-001).") app.add_typer(cn_app, name="cn") @cn_app.command("build-model") def cn_build_model_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], fasta: Annotated[Path, typer.Option("--fasta", help="Reference FASTA of the manifest's build (required).")], out: Annotated[Path, typer.Option("--out", help="Output model file (.simurg-cyto).")], egt: Annotated[Path | None, typer.Option("--egt", help="Cluster file: population B-allele frequencies.")] = None, reference_gtcs: Annotated[ Path | None, typer.Option("--reference-gtcs", help="GTCs of normal samples (>= 24) for the batch reference.") ] = None, max_reference: Annotated[int, typer.Option("--max-reference", min=1, help="Random subset size.")] = 100, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Cyto model: GC covariates, PFB and an optional batch reference per probe (CYTO-038, IN-051).""" from simurg.pipeline.cn import run_build_model rep = _reporter("cn build-model", out.parent, log_json, debug, quiet) run_guarded(rep, lambda: run_build_model(bpm, fasta, out, rep, egt, "auto", reference_gtcs, max_reference)) def _kb(v: int) -> int: if not 0 <= v <= 250_000: raise typer.BadParameter("sizes are 0-250000 kb") return v * 1000 @cn_app.command("call") def cn_call_cmd( bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], model: Annotated[Path, typer.Option("--model", help="Cyto model from cn build-model (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None, summary: Annotated[ Path | None, typer.Option("--summary", help="sample_summary.tsv (default: next to GTCs).") ] = None, fasta: Annotated[Path | None, typer.Option("--fasta", help="FASTA for REF padding bases.")] = None, min_del_probes: Annotated[int, typer.Option("--min-del-probes", min=0, max=250000)] = 10, min_del_size: Annotated[int, typer.Option("--min-del-size", help="kb")] = 25, min_dup_probes: Annotated[int, typer.Option("--min-dup-probes", min=0, max=250000)] = 10, min_dup_size: Annotated[int, typer.Option("--min-dup-size", help="kb")] = 50, min_loh_probes: Annotated[int, typer.Option("--min-loh-probes", min=0, max=250000)] = 500, min_loh_size: Annotated[int, typer.Option("--min-loh-size", help="kb")] = 3000, min_gainloh_probes: Annotated[int, typer.Option("--min-gainloh-probes", min=0, max=250000)] = 300, min_gainloh_size: Annotated[int, typer.Option("--min-gainloh-size", help="kb")] = 500, min_mdel_probes: Annotated[int, typer.Option("--min-mdel-probes", min=0, max=250000)] = 10, min_mdel_size: Annotated[int, typer.Option("--min-mdel-size", help="kb")] = 25, min_mdup_probes: Annotated[int, typer.Option("--min-mdup-probes", min=0, max=250000)] = 10, min_mdup_size: Annotated[int, typer.Option("--min-mdup-size", help="kb")] = 50, min_mloh_probes: Annotated[int, typer.Option("--min-mloh-probes", min=0, max=250000)] = 500, min_mloh_size: Annotated[int, typer.Option("--min-mloh-size", help="kb")] = 3000, smoothing: Annotated[int, typer.Option("--smoothing", min=0, max=1000)] = 5, max_mosaic_fraction: Annotated[float, typer.Option("--max-mosaic-fraction")] = 1.0, min_mosaic_fraction: Annotated[float, typer.Option("--min-mosaic-fraction")] = 0.10, filter_qual: Annotated[int, typer.Option("--filter-qual", min=0, max=60)] = 10, cn_cap: Annotated[int, typer.Option("--cn-cap", min=3, max=8)] = 4, gc_correct: Annotated[bool, typer.Option("--gc-correct/--no-gc-correct")] = True, unknown_sex_policy: Annotated[str, typer.Option("--unknown-sex-policy")] = "skip-sex-chromosomes", gt_mode: Annotated[str, typer.Option("--gt-mode", help="vcf|compat")] = "vcf", svlen_mode: Annotated[str, typer.Option("--svlen-mode", help="signed|positive")] = "signed", mask_bed: Annotated[Path | None, typer.Option("--mask-bed", help="Mask BED (default: built-in).")] = None, no_bgzip: Annotated[bool, typer.Option("--no-bgzip", help="Plain VCF (cn annotate then refuses it).")] = False, mosaic_prior: Annotated[ float, typer.Option("--mosaic-prior", help="Advanced: mosaic states' switch mass.") ] = 0.004, switch_distance: Annotated[ float, typer.Option("--switch-distance", help="Advanced: HMM switch distance D (bp).") ] = 5e8, probe_lrr_sd: Annotated[ bool, typer.Option("--probe-lrr-sd/--no-probe-lrr-sd", help="Advanced: per-probe LRR noise from the model.") ] = True, threads: ThreadsOpt = 0, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """GTCs -> per-sample CNV/LOH/mosaic VCF, segments and summary (CLI-028).""" from simurg.cn.events import Thresholds from simurg.pipeline.cn import CnCallOptions, run_cn_call if unknown_sex_policy not in ("skip-sex-chromosomes", "female", "male", "error"): raise typer.BadParameter("--unknown-sex-policy: skip-sex-chromosomes|female|male|error") if gt_mode not in ("vcf", "compat") or svlen_mode not in ("signed", "positive"): raise typer.BadParameter("--gt-mode vcf|compat, --svlen-mode signed|positive") if not 0 < max_mosaic_fraction <= 1 or not 0 <= min_mosaic_fraction < 1: raise typer.BadParameter("mosaic fractions: 0 < max <= 1, 0 <= min < 1") th = Thresholds( probes={"del": min_del_probes, "dup": min_dup_probes, "loh": min_loh_probes, "gainloh": min_gainloh_probes, "mdel": min_mdel_probes, "mdup": min_mdup_probes, "mloh": min_mloh_probes}, size_bp={"del": _kb(min_del_size), "dup": _kb(min_dup_size), "loh": _kb(min_loh_size), "gainloh": _kb(min_gainloh_size), "mdel": _kb(min_mdel_size), "mdup": _kb(min_mdup_size), "mloh": _kb(min_mloh_size)}, max_mosaic_fraction=max_mosaic_fraction, min_mosaic_fraction=min_mosaic_fraction, filter_qual=filter_qual, ) # fmt: skip rep = _reporter("cn call", out, log_json, debug, quiet) opts = CnCallOptions( out=out, bpm=bpm, model=model, gtcs=gtcs, summary=summary, fasta=fasta, thresholds=th, smoothing=smoothing, cn_cap=cn_cap, gc_correct=gc_correct, unknown_sex_policy=unknown_sex_policy, gt_mode=gt_mode, svlen_mode=svlen_mode, mask_bed=mask_bed, bgzip=not no_bgzip, mosaic_prior=mosaic_prior, switch_d=switch_distance, probe_lrr_sd=probe_lrr_sd, threads=threads, force=force, argv=sys.argv[1:], ) # fmt: skip run_guarded(rep, lambda: run_cn_call(opts, rep)) @cn_app.command("annotate") def cn_annotate_cmd( vcfs: Annotated[Path, typer.Option("--vcfs", help="Output folder of cn call (required).")], bundle: Annotated[ Path | None, typer.Option("--annotation-bundle", help="Bundle from scripts/build_annotation_bundle.sh.") ] = None, out: Annotated[Path | None, typer.Option("--out", help="Output folder (default: --vcfs).")] = None, min_qual: Annotated[int, typer.Option("--min-qual", min=0, max=60)] = 20, size_basis: Annotated[str, typer.Option("--size-basis", help="effective|raw")] = "effective", iscn_min_size: Annotated[int, typer.Option("--iscn-min-size", min=0, help="bp")] = 0, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """cn call output -> per-sample JSON report with ISCN 2020 strings, cytobands and genes (CLI-029).""" from simurg.pipeline.cn import CnAnnotateOptions, run_cn_annotate if size_basis not in ("effective", "raw"): raise typer.BadParameter("--size-basis effective|raw") rep = _reporter("cn annotate", out or vcfs, log_json, debug, quiet) opts = CnAnnotateOptions(vcfs=vcfs, bundle=bundle, out=out, min_qual=min_qual, size_basis=size_basis, iscn_min_size=iscn_min_size, force=force) # fmt: skip run_guarded(rep, lambda: run_cn_annotate(opts, rep)) # -------------------------------------------------------------------------------------------------------- pgx_app = typer.Typer(help="Pharmacogenomics: targeted copy number, star alleles, phenotypes (v2).") app.add_typer(pgx_app, name="pgx") pgx_cn_app = typer.Typer(help="Targeted PGx copy number (PCN-001).") pgx_app.add_typer(pgx_cn_app, name="copy-number") pgx_db_app = typer.Typer(help="PGx allele-definition database (PGX-001).") pgx_app.add_typer(pgx_db_app, name="db") @pgx_db_app.command("build") def pgx_db_build_cmd( sources: Annotated[Path, typer.Option("--sources", help="Folder from scripts/pgx/fetch_sources.sh (required).")], fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA for variant normalisation (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder for pgx_db--.sqlite.")] = Path("."), manifests: Annotated[ list[Path] | None, typer.Option("--manifests", help="BPM manifest(s) for the probe map (repeatable).") ] = None, csv: Annotated[list[Path] | None, typer.Option("--csv", help="CSV manifest(s), same order (indels).")] = None, mask: Annotated[list[Path] | None, typer.Option("--mask", help="Probe mask TSV(s) (PGX-010).")] = None, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Build the PGx DB from public sources (CLI-035, PGX-002).""" from simurg.pipeline.pgx_db import PgxDbOptions, run_db_build rep = _reporter("pgx db build", out, log_json, debug, quiet) o = PgxDbOptions(sources=sources, fasta=fasta, out=out, manifests=manifests or [], csvs=csv or [], masks=mask or [], force=force) # fmt: skip run_guarded(rep, lambda: run_db_build(o, rep)) pgx_star_app = typer.Typer(help="Star-allele calling (PGX-011).") pgx_app.add_typer(pgx_star_app, name="star-allele") @pgx_star_app.command("call") def pgx_star_call_cmd( vcfs: Annotated[list[Path], typer.Option("--vcfs", help="Folder(s) with SNV + PGx CN VCFs (repeatable).")], db: Annotated[Path, typer.Option("--db", help="PGx DB from pgx db build (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], population: Annotated[str, typer.Option("--population", help="Frequency prior population.")] = "global", max_candidates: Annotated[int, typer.Option("--max-candidates", min=1)] = 10, feasible_ratio: Annotated[float, typer.Option("--feasible-ratio", min=0.0, max=1.0)] = 0.5, low_gs: Annotated[float, typer.Option("--low-gs")] = 0.15, no_cn: Annotated[bool, typer.Option("--no-cn", help="Run without PGx CN VCFs.")] = False, allow_foreign_vcf: Annotated[bool, typer.Option("--allow-foreign-vcf")] = False, gene: Annotated[list[str] | None, typer.Option("--gene", help="Restrict to gene(s).")] = None, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Call star-allele diplotypes (CLI-033, PGX-012 ... PGX-021).""" from simurg.pipeline.pgx_star import StarOptions, run_star_call rep = _reporter("pgx star-allele call", out, log_json, debug, quiet) o = StarOptions(vcfs=vcfs, db=db, out=out, population=population, max_candidates=max_candidates, feasible_ratio=feasible_ratio, low_gs=low_gs, no_cn=no_cn, allow_foreign_vcf=allow_foreign_vcf, genes=gene or [], force=force) # fmt: skip run_guarded(rep, lambda: run_star_call(o, rep)) def _pgx_phenotype(calls: Path, db: Path, out: Path | None, guideline: str, force: bool, log_json: bool, debug: bool, quiet: bool, name: str) -> None: # fmt: skip from simurg.pipeline.pgx_pheno import PhenoOptions, run_phenotype if guideline not in ("cpic", "dpwg"): raise typer.BadParameter("--guideline cpic|dpwg") rep = _reporter(name, out or calls, log_json, debug, quiet) o = PhenoOptions(calls=calls, db=db, out=out, guideline=guideline, force=force) run_guarded(rep, lambda: run_phenotype(o, rep)) @pgx_app.command("phenotype") def pgx_phenotype_cmd( calls: Annotated[Path, typer.Option("--calls", help="Output folder of pgx star-allele call (required).")], db: Annotated[Path, typer.Option("--db", help="PGx DB (required).")], out: Annotated[Path | None, typer.Option("--out", help="Output folder (default: --calls).")] = None, guideline: Annotated[str, typer.Option("--guideline", help="cpic|dpwg")] = "cpic", force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Phenotypes from star-allele calls: per-sample PGx JSON and diplotype_annotations.tsv (CLI-034).""" _pgx_phenotype(calls, db, out, guideline, force, log_json, debug, quiet, "pgx phenotype") @pgx_app.command("annotate", hidden=True) def pgx_annotate_cmd( calls: Annotated[Path, typer.Option("--calls")], db: Annotated[Path, typer.Option("--db")], out: Annotated[Path | None, typer.Option("--out")] = None, guideline: Annotated[str, typer.Option("--guideline")] = "cpic", force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Alias of pgx phenotype (CLI-034).""" _pgx_phenotype(calls, db, out, guideline, force, log_json, debug, quiet, "pgx annotate") @pgx_cn_app.command("train") def pgx_cn_train_cmd( gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for training GTCs (required).")], bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder.")] = Path("."), mask: Annotated[Path | None, typer.Option("--mask", help="Probe mask TSV (IN-052).")] = None, regions: Annotated[Path | None, typer.Option("--regions", help="Region BED (default: built-in v1).")] = None, summary: Annotated[Path | None, typer.Option("--summary", help="sample_summary.tsv (assay control).")] = None, platform: Annotated[str, typer.Option("--platform", help="auto|LCG|EX")] = "auto", min_samples: Annotated[int, typer.Option("--min-samples", min=1)] = 96, min_passing: Annotated[int, typer.Option("--min-passing", min=1)] = 90, max_samples: Annotated[int, typer.Option("--max-samples", min=1)] = 3000, train_lrr_sd_max: Annotated[float, typer.Option("--train-lrr-sd-max")] = 0.20, train_call_rate_min: Annotated[float, typer.Option("--train-call-rate-min")] = 0.99, train_assay_control_min: Annotated[float, typer.Option("--train-assay-control-min")] = 1.0, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Train a PGx CN model (CLI-032, PCN-005).""" from simurg.pipeline.pgx_cn import PgxTrainOptions, run_train if platform not in ("auto", "LCG", "EX"): raise typer.BadParameter("--platform auto|LCG|EX") rep = _reporter("pgx copy-number train", out, log_json, debug, quiet) o = PgxTrainOptions(gtcs=gtcs, bpm=bpm, fasta=fasta, out=out, mask=mask, regions=regions, summary=summary, platform=platform, min_samples=min_samples, min_passing=min_passing, max_samples=max_samples, lrr_sd_max=train_lrr_sd_max, call_rate_min=train_call_rate_min, assay_control_min=train_assay_control_min, force=force) # fmt: skip run_guarded(rep, lambda: run_train(o, rep)) @pgx_cn_app.command("call") def pgx_cn_call_cmd( gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for GTC files (required).")], model: Annotated[Path, typer.Option("--model", help="PGx CN model (.simurg-pgxcn, required).")], fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA with alt contigs (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], min_batch: Annotated[int, typer.Option("--min-batch", min=1)] = 24, min_passing: Annotated[int, typer.Option("--min-passing", min=1)] = 22, recommended_passing: Annotated[int, typer.Option("--recommended-passing", min=1)] = 96, lrr_sd_max: Annotated[float, typer.Option("--lrr-sd-max")] = 0.20, cluster_mismatch: Annotated[str, typer.Option("--cluster-mismatch", help="error|warn")] = "error", no_bgzip: Annotated[bool, typer.Option("--no-bgzip")] = False, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Call targeted PGx copy number per region (CLI-031, PCN-009).""" from simurg.pipeline.pgx_cn import PgxCallOptions, run_call if cluster_mismatch not in ("error", "warn"): raise typer.BadParameter("--cluster-mismatch error|warn") if min_passing > min_batch: raise typer.BadParameter("--min-passing must be <= --min-batch") rep = _reporter("pgx copy-number call", out, log_json, debug, quiet) o = PgxCallOptions(gtcs=gtcs, model=model, fasta=fasta, out=out, min_batch=min_batch, min_passing=min_passing, recommended_passing=recommended_passing, lrr_sd_max=lrr_sd_max, cluster_mismatch=cluster_mismatch, bgzip=not no_bgzip, force=force) # fmt: skip run_guarded(rep, lambda: run_call(o, rep)) # -------------------------------------------------------------------------------------------------------- methyl_app = typer.Typer(help="Methylation arrays: control QC, detection, beta values (MET-001).") app.add_typer(methyl_app, name="methyl") @methyl_app.command("qc") def methyl_qc_cmd( manifest: Annotated[Path, typer.Option("--manifest", help="Methylation CSV manifest (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None, sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, config: Annotated[Path | None, typer.Option("--config", help="QC threshold YAML.")] = None, profile: Annotated[str, typer.Option("--profile", help="QC profile name, 'auto' or 'none'.")] = "auto", threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None, bg_offset: Annotated[float | None, typer.Option("--bg-offset", help="Background offset x (3000).")] = None, bs_set: Annotated[str, typer.Option("--bs-conv-red-bg-set", help="U3-U5 (default) or U4-U6.")] = "U3-U5", pval: Annotated[float, typer.Option("--pval-threshold", help="Detection p threshold.")] = 0.05, beta_offset: Annotated[float, typer.Option("--beta-offset", help="alpha in M/(M+U+alpha).")] = 0.0, dye_bias: Annotated[str, typer.Option("--dye-bias", help="linear|nonlinear")] = "linear", min_prop: Annotated[ float | None, typer.Option("--min-prop-passing", help="Fail samples below this detection rate.") ] = None, pca_full: Annotated[int, typer.Option("--pca-max-samples-full", min=2)] = 100, pca_probes: Annotated[int, typer.Option("--pca-probes", min=2)] = 10_000, threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Methylation control QC: 21 control metrics, detection p-values, noob/dye-bias corrected β (MET-001).""" from simurg.pipeline.methyl import MethylQcOptions, run_methyl_qc if dye_bias not in ("linear", "nonlinear"): raise typer.BadParameter("--dye-bias must be linear or nonlinear") if bs_set not in ("U3-U5", "U4-U6"): raise typer.BadParameter("--bs-conv-red-bg-set must be U3-U5 or U4-U6") if not 0 < pval <= 1: raise typer.BadParameter("--pval-threshold must be in (0, 1]") rep = _reporter("methyl qc", out, log_json, debug, quiet) opts = MethylQcOptions( out=out, manifest=manifest, idats=idats, sheet=sheet, config=config, profile=None if profile == "none" else profile, thresholds=threshold or [], bg_offset=bg_offset, bs_conv_red_bg_set=bs_set, pval_threshold=pval, beta_offset=beta_offset, dye_bias=dye_bias, min_prop_passing=min_prop, pca_max_samples_full=pca_full, pca_probes=pca_probes, threads=threads, force=force, argv=sys.argv[1:], ) # fmt: skip run_guarded(rep, lambda: run_methyl_qc(opts, rep)) @methyl_app.command("convert-manifest") def methyl_convert_cmd( manifest: Annotated[Path, typer.Option("--manifest", help="Methylation CSV manifest (required).")], out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Methylation manifest -> probes.tsv and controls.tsv (MET-003).""" from simurg.pipeline.methyl import run_convert_manifest rep = _reporter("methyl convert-manifest", out, log_json, debug, quiet) run_guarded(rep, lambda: run_convert_manifest(manifest, out, force, rep)) def _qc_cohort( command: str, out: Path, bpm: Path, gtcs: Path | None, sheet: Path | None, egt: Path | None, threads: int, force: bool, log_json: bool, debug: bool, quiet: bool, *, all_pairs: bool = False, pedigree_from_sheet: bool = True, median: str = "reference", config: Path | None = None, profile: str | None = None, threshold_call_rate: float | None = None, threshold_lrr_sd: float | None = None, threshold: list[str] | None = None, ) -> None: from simurg.pipeline.qc import QcCohortOptions, run_qc_cohort from simurg.qc.thresholds import parse_cli_pairs, resolve_thresholds if median not in ("reference", "exact"): raise typer.BadParameter("--median must be reference or exact") rep = _reporter(f"qc {command}", out, log_json, debug, quiet) def body() -> int: cli: dict[str, Any] = parse_cli_pairs(list(threshold or [])) if threshold_call_rate is not None: cli["callRate"] = threshold_call_rate if threshold_lrr_sd is not None: cli["logRDev"] = threshold_lrr_sd th = resolve_thresholds(config=config, profile=profile, cli=cli) opts = QcCohortOptions( command=command, out=out, bpm=bpm, gtcs=gtcs, sheet=sheet, egt=egt, threads=threads, force=force, all_pairs=all_pairs, pedigree_from_sheet=pedigree_from_sheet, median=median, thresholds=th, argv=sys.argv[1:], ) return run_qc_cohort(opts, rep) run_guarded(rep, body) GtcsOpt = Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] SheetOpt = Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] BpmOpt = Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")] OutOpt = Annotated[Path, typer.Option("--out", help="Output folder (required).")] ThreadsOpt = Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] ForceOpt = Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] MedianOpt = Annotated[str, typer.Option("--median", help="reference (default, bit-compatible) or exact.")] @qc_app.command("samples") def qc_samples_cmd( out: OutOpt, bpm: BpmOpt, gtcs: GtcsOpt = None, sheet: SheetOpt = None, config: Annotated[Path | None, typer.Option("--config", help="QC threshold YAML.")] = None, profile: Annotated[str | None, typer.Option("--profile", help="QC profile name or 'auto'.")] = None, threshold_call_rate: Annotated[float | None, typer.Option("--threshold-call-rate")] = None, threshold_lrr_sd: Annotated[float | None, typer.Option("--threshold-lrr-sd")] = None, threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None, all_pairs: Annotated[bool, typer.Option("--all-pairs", help="Write every sample pair.")] = False, median: MedianOpt = "reference", threads: ThreadsOpt = 0, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Sample QC from GTCs: het rate, contamination, relatedness, flags (sample_qc.tsv).""" _qc_cohort( "samples", out, bpm, gtcs, sheet, None, threads, force, log_json, debug, quiet, all_pairs=all_pairs, median=median, config=config, profile=profile, threshold_call_rate=threshold_call_rate, threshold_lrr_sd=threshold_lrr_sd, threshold=threshold, ) @qc_app.command("contamination") def qc_contamination_cmd( out: OutOpt, bpm: BpmOpt, gtcs: GtcsOpt = None, sheet: SheetOpt = None, median: MedianOpt = "reference", threads: ThreadsOpt = 0, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """BAF-regression contamination estimate per sample (contamination.tsv).""" _qc_cohort("contamination", out, bpm, gtcs, sheet, None, threads, force, log_json, debug, quiet, median=median) @qc_app.command("relatedness") def qc_relatedness_cmd( out: OutOpt, bpm: BpmOpt, gtcs: GtcsOpt = None, sheet: SheetOpt = None, pedigree_from_sheet: Annotated[ bool, typer.Option("--pedigree-from-sheet/--no-pedigree", help="Use Parent1/Parent2/Replicate.") ] = True, all_pairs: Annotated[bool, typer.Option("--all-pairs", help="Write every sample pair.")] = False, threads: ThreadsOpt = 0, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """KING-robust kinship, duplicates and Mendelian errors (relatedness.tsv, mendel.tsv).""" _qc_cohort( "relatedness", out, bpm, gtcs, sheet, None, threads, force, log_json, debug, quiet, all_pairs=all_pairs, pedigree_from_sheet=pedigree_from_sheet, ) @qc_app.command("loci") def qc_loci_cmd( out: OutOpt, bpm: BpmOpt, egt: Annotated[Path, typer.Option("--egt", help="EGT cluster file (required).")], gtcs: GtcsOpt = None, sheet: SheetOpt = None, threads: ThreadsOpt = 0, force: ForceOpt = False, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Per-locus summary statistics across samples (locus_summary.tsv, OUT-072).""" _qc_cohort("loci", out, bpm, gtcs, sheet, egt, threads, force, log_json, debug, quiet) @qc_app.command("report") def qc_report_cmd( data: Annotated[list[str], typer.Option("--data", help="Dataset folder(s): one, a comma list, or a parent.")], out: Annotated[Path, typer.Option("--out", help="Output folder.")] = Path("."), label: Annotated[ list[str] | None, typer.Option("--label", help="Dataset labels (list or discovery order).") ] = None, config: Annotated[Path | None, typer.Option("--config", help="QC threshold YAML (QC-020).")] = None, profile: Annotated[str | None, typer.Option("--profile", help="QC profile name or 'auto'.")] = None, threshold_call_rate: Annotated[float | None, typer.Option("--threshold-call-rate")] = None, threshold_lrr_sd: Annotated[float | None, typer.Option("--threshold-lrr-sd")] = None, threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None, output_format: Annotated[str, typer.Option("--output-format", help="xlsx|csv|both")] = "xlsx", report_name: Annotated[str | None, typer.Option("--report-name", help="Fixed output file stem.")] = None, drop_metadata: Annotated[ list[str] | None, typer.Option("--drop-metadata", help="Metadata columns to drop.") ] = None, keep_metadata: Annotated[list[str] | None, typer.Option("--keep-metadata", help="Metadata allowlist.")] = None, no_metadata: Annotated[bool, typer.Option("--no-metadata", help="Embed no metadata columns.")] = False, redact_paths: Annotated[bool, typer.Option("--redact-paths", help="Blank sample_folder.")] = False, split_by: Annotated[str, typer.Option("--split-by", help="none|dataset")] = "none", log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Self-contained HTML QC report and QC table from dataset folders (QC-040 … QC-063).""" from simurg.pipeline.qc_report import QcReportOptions, run_qc_report from simurg.qc.thresholds import parse_cli_pairs, resolve_thresholds rep = _reporter("qc report", None, log_json, debug, quiet) def body() -> int: cli: dict[str, Any] = parse_cli_pairs(list(threshold or [])) if threshold_call_rate is not None: cli["callRate"] = threshold_call_rate if threshold_lrr_sd is not None: cli["logRDev"] = threshold_lrr_sd th = resolve_thresholds(config=config, profile=profile, cli=cli) opts = QcReportOptions( data=list(data), out=out, labels=list(label) if label else None, output_format=output_format, report_name=report_name, drop_metadata=list(drop_metadata or []), keep_metadata=list(keep_metadata) if keep_metadata else None, no_metadata=no_metadata, redact_paths=redact_paths, split_by=split_by, thresholds=th, ) return run_qc_report(opts, rep) run_guarded(rep, body) @qc_app.command("suggest-thresholds") def qc_suggest_cmd( data: Annotated[list[Path], typer.Option("--data", help="Dataset folder(s) with controls_qc.tsv.")], out: Annotated[Path, typer.Option("--out", help="Output YAML file (required).")], k: Annotated[float, typer.Option("--k", min=0.0, help="Spread multiplier.")] = 4.0, min_samples: Annotated[int, typer.Option("--min-samples", min=1, help="Minimum samples per metric.")] = 30, config: Annotated[Path | None, typer.Option("--config", help="Current thresholds (YAML).")] = None, profile: Annotated[str | None, typer.Option("--profile", help="Current profile.")] = None, log_json: LogJson = False, debug: Debug = False, quiet: Quiet = False, ) -> None: """Data-driven QC threshold suggestion (QC-026).""" from simurg.pipeline.qc import run_suggest_thresholds from simurg.qc.thresholds import resolve_thresholds rep = _reporter("qc suggest-thresholds", None, log_json, debug, quiet) folders = [Path(p) for d in data for p in str(d).split(",") if p] def body() -> int: th = resolve_thresholds(config=config, profile=profile) return run_suggest_thresholds(folders, out, th, rep, k=k, min_samples=min_samples) run_guarded(rep, body) if __name__ == "__main__": # pragma: no cover app()