kemalandClaude Opus 5.5
fbc4080a50
- pipeline/qc_report: dataset discovery/labels/types (QC-041..043), robust TSV/CSV reading (QC-061), join + unmatched warning (QC-040), privacy options (QC-046), xlsx (4 sheets, conditional fills) / csv table (QC-045), byte-deterministic HTML with gzip+base64 embedded data (QC-044, QC-047) - qc/report: template, CSS and JS app — Control Dashboard, Automated QC (table with sort/filter/search/ paging/column menu/export, metric scatters with threshold lines), Sample Heatmaps (chip/plate, continuous/binary), Trend Analysis (grouped summary, Tukey box plots, counts, sex, PCA), thresholds modal (apply / download / upload YAML), linked selection across plots and table (QC-048..060) - bundled plotly.js 2.35.2 (MIT); openpyxl dependency (MIT) - tests: discovery/labels/errors/determinism/xlsx + headless-browser self-test - known_unknowns K-020: remaining report deviations Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
689 lines
28 KiB
Python
689 lines
28 KiB
Python
"""``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 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())
|
|
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))
|
|
|
|
|
|
# --------------------------------------------------------------------------------------------------------
|
|
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))
|
|
|
|
|
|
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()
|