"""Synthetic software verification of the anchored depth caller. This experiment deliberately separates: 1. normal-structure fitting; 2. null/threshold calibration; 3. development comparison of a small, predeclared component grid; and 4. one locked same-assay verification stream; 5. within-assay low-depth/site/mosaic stress; and 6. transport to a deliberately changed assay with aligned target IDs. The simulator is count-level and is not a substitute for physical reference materials, analytical validation, or clinical validation. Event truth is always evaluated against its complete target list; it is never clipped to the caller's callable mask. Callability is reported as a separate endpoint. """ from __future__ import annotations import argparse import hashlib import json import math import platform import time from dataclasses import dataclass from datetime import datetime, timezone from pathlib import Path from typing import Sequence import numpy as np import sklearn from lgasieve.clinical_depth import AnchoredDepthModel, DepthTarget from lgasieve.panel import Panel from lgasieve.provenance import atomic_write_json SCHEMA = "lgasieve.synthetic-software-verification.v2" BASE_SEED = 20260731 COMPONENT_GRID = (0, 4, 8, 12) ALPHA_GRID = (0.01, 0.02, 0.05) @dataclass(frozen=True) class Event: gene: str kind: str indices: tuple[int, ...] amplitude: str factor: float event_class: str @dataclass class Assay: targets: tuple[DepthTarget, ...] efficiency: np.ndarray gc: np.ndarray latent: np.ndarray instability: np.ndarray reportable_by_gene: dict[str, tuple[int, ...]] unstable_reportable: tuple[int, ...] def sha256_file(path: Path) -> str: digest = hashlib.sha256() with Path(path).open("rb") as handle: for block in iter(lambda: handle.read(1024 * 1024), b""): digest.update(block) return digest.hexdigest() def sha256_json(value: object) -> str: encoded = json.dumps( value, sort_keys=True, separators=(",", ":"), ensure_ascii=False ).encode("utf-8") return hashlib.sha256(encoded).hexdigest() def build_assay( panel: Panel, *, n_controls: int = 500, seed: int = BASE_SEED, ) -> Assay: rng = np.random.default_rng(seed) controls = [ DepthTarget( target_id=f"CTRL_{index + 1:04d}", gene="CONTROL", exon=str(index + 1), order=index + 1, reportable=False, ) for index in range(n_controls) ] reportable: list[DepthTarget] = [] panel_indices: list[int] = [] for gene in sorted(panel.order): for panel_index in sorted( panel.order[gene], key=lambda index: panel.targets[index].exon ): target = panel.targets[panel_index] reportable.append( DepthTarget( target_id=target.id, gene=target.gene, exon=target.legacy_exon, order=int(target.exon), reportable=True, ) ) panel_indices.append(panel_index) targets = tuple(controls + reportable) n_targets = len(targets) control_gc = np.clip(rng.beta(5.5, 6.0, n_controls), 0.20, 0.78) reportable_gc = np.asarray( [panel.targets[index].gc for index in panel_indices], dtype=float ) gc = np.concatenate((control_gc, reportable_gc)) fixed = rng.normal(0.0, 0.22, n_targets) fixed += -3.2 * (gc - 0.47) ** 2 efficiency = np.exp(fixed) efficiency /= np.median(efficiency) latent = rng.normal(0.0, 0.075, (8, n_targets)) # Smooth loadings over adjacent BRCA bins make the simulation harder than # independent target noise and reward proper factor correction. for component in range(latent.shape[0]): tail = latent[component, n_controls:] for _ in range(3): tail[1:-1] = 0.25 * tail[:-2] + 0.50 * tail[1:-1] + 0.25 * tail[2:] latent[component, n_controls:] = tail instability = np.full(n_targets, 0.035) unstable_control = rng.choice(n_controls, size=12, replace=False) instability[unstable_control] = rng.uniform(0.08, 0.13, unstable_control.size) reportable_pool = np.arange(n_controls, n_targets) unstable_reportable = tuple( sorted( int(value) for value in rng.choice(reportable_pool, size=4, replace=False) ) ) instability[list(unstable_reportable)] = rng.uniform(0.11, 0.17, 4) reportable_by_gene: dict[str, tuple[int, ...]] = {} for gene in sorted(panel.order): reportable_by_gene[gene] = tuple( index for index, target in enumerate(targets) if target.reportable and target.gene == gene ) return Assay( targets=targets, efficiency=efficiency, gc=gc, latent=latent, instability=instability, reportable_by_gene=reportable_by_gene, unstable_reportable=unstable_reportable, ) def build_transport_shift_assay( source: Assay, *, seed: int, ) -> tuple[Assay, dict]: """Create a changed assay/chemistry while preserving target identity. The shift is deliberately target-specific rather than a scalar site offset: target efficiencies, latent-factor loadings, and target instability all change. GC and target coordinates remain fixed because the same target IDs are being transported. This is a software stress construction, not a model of a named commercial assay. """ rng = np.random.default_rng(seed) n_targets = len(source.targets) n_components = source.latent.shape[0] reportable = np.asarray( [ index for index, target in enumerate(source.targets) if target.reportable ], dtype=int, ) controls = np.asarray( [ index for index, target in enumerate(source.targets) if not target.reportable ], dtype=int, ) gc_delta = source.gc - float(np.median(source.gc)) log_efficiency_delta = ( rng.normal(0.0, 0.16, n_targets) + 0.55 * gc_delta + 0.06 * np.sin(np.arange(n_targets) / 17.0) ) # Make reportable-bin transport nontrivial without changing target names. log_efficiency_delta[reportable] += rng.normal( 0.0, 0.11, reportable.size ) efficiency = source.efficiency * np.exp(log_efficiency_delta) efficiency /= np.median(efficiency) mixing = np.eye(n_components) + rng.normal( 0.0, 0.16, (n_components, n_components) ) latent = mixing @ source.latent latent += rng.normal(0.0, 0.025, source.latent.shape) latent[:, reportable] += rng.normal( 0.0, 0.025, (n_components, reportable.size) ) instability = source.instability * rng.lognormal( mean=math.log(1.30), sigma=0.22, size=n_targets ) shifted_controls = rng.choice( controls, size=min(24, controls.size), replace=False ) shifted_reportable = np.asarray( rng.choice( reportable, size=min(8, reportable.size), replace=False ), dtype=int, ) instability[shifted_controls] = np.maximum( instability[shifted_controls], rng.uniform(0.10, 0.18, shifted_controls.size), ) instability[shifted_reportable] = np.maximum( instability[shifted_reportable], rng.uniform(0.15, 0.24, shifted_reportable.size), ) shifted = Assay( targets=source.targets, efficiency=efficiency, gc=source.gc.copy(), latent=latent, instability=instability, reportable_by_gene={ gene: tuple(indices) for gene, indices in source.reportable_by_gene.items() }, unstable_reportable=tuple( sorted(int(index) for index in shifted_reportable) ), ) source_ids = [target.target_id for target in source.targets] shifted_ids = [target.target_id for target in shifted.targets] source_log_efficiency = np.log(source.efficiency) shifted_log_efficiency = np.log(shifted.efficiency) metadata = { "construction": ( "target-specific efficiency, latent-loading, and instability " "shift; same target IDs and GC values" ), "seed": seed, "target_ids_aligned": source_ids == shifted_ids, "target_ids_sha256": sha256_json(source_ids), "efficiency_log2_ratio_rms": float( np.sqrt( np.mean( np.log2(shifted.efficiency / source.efficiency) ** 2 ) ) ), "log_efficiency_correlation": float( np.corrcoef( source_log_efficiency, shifted_log_efficiency )[0, 1] ), "latent_loading_rms_delta": float( np.sqrt(np.mean((shifted.latent - source.latent) ** 2)) ), "latent_loading_correlation": float( np.corrcoef( source.latent.ravel(), shifted.latent.ravel() )[0, 1] ), "source_median_instability": float( np.median(source.instability) ), "shifted_median_instability": float( np.median(shifted.instability) ), "forced_unstable_reportable_target_ids": [ shifted.targets[index].target_id for index in shifted.unstable_reportable ], } return shifted, metadata def event_design( assay: Assay, *, n: int, seed: int, include_mosaic: bool = False, ) -> list[Event]: rng = np.random.default_rng(seed) widths = (1, 2, 3, 5) amplitudes = ( ( ("germline", {"DEL": 0.50, "DUP": 1.50}), ("mosaic_50pct", {"DEL": 0.75, "DUP": 1.25}), ) if include_mosaic else (("germline", {"DEL": 0.50, "DUP": 1.50}),) ) placements: list[tuple[str, tuple[int, ...], str]] = [] for gene, indices in assay.reportable_by_gene.items(): for width in widths: for start in range(len(indices) - width + 1): placements.append( ( gene, tuple(indices[start : start + width]), "single" if width == 1 else f"width_{width}", ) ) placements.append((gene, tuple(indices), "whole_gene")) candidates = [ Event( gene=gene, kind=kind, indices=indices, amplitude=amplitude, factor=float(factors[kind]), event_class=event_class, ) for gene, indices, event_class in placements for kind in ("DEL", "DUP") for amplitude, factors in amplitudes ] order = rng.integers(0, len(candidates), size=n) return [candidates[int(index)] for index in order] def simulate_counts( assay: Assay, n_samples: int, *, seed: int, events: Sequence[Event | None] | None = None, site_shift: float = 0.0, median_depth: float = 360.0, low_depth_fraction: float = 0.0, artifact_sample_fraction: float = 0.015, ) -> tuple[np.ndarray, list[dict]]: rng = np.random.default_rng(seed) n_targets = len(assay.targets) if events is None: event_rows: list[Event | None] = [None] * n_samples else: event_rows = list(events) if len(event_rows) != n_samples: raise ValueError("one event row is required per simulated sample") depth = np.exp(rng.normal(math.log(median_depth), 0.28, n_samples)) if low_depth_fraction: low = rng.random(n_samples) < low_depth_fraction depth[low] *= rng.uniform(0.18, 0.45, int(low.sum())) factor_score = rng.normal(0.0, 1.0, (n_samples, assay.latent.shape[0])) if site_shift: factor_score[:, :3] += site_shift * np.asarray([0.8, -0.6, 0.5]) gc_linear = rng.normal(0.0 + 0.10 * site_shift, 0.42, n_samples) gc_quad = rng.normal(0.0 - 0.08 * site_shift, 0.30, n_samples) gc_delta = assay.gc - 0.47 log_rate = ( np.log(depth)[:, None] + np.log(assay.efficiency)[None, :] + factor_score @ assay.latent + gc_linear[:, None] * gc_delta[None, :] + gc_quad[:, None] * 5.0 * gc_delta[None, :] ** 2 + rng.normal(0.0, assay.instability, (n_samples, n_targets)) ) rate = np.exp(log_rate) truth: list[dict] = [] for sample_index, event in enumerate(event_rows): if event is not None: rate[sample_index, list(event.indices)] *= event.factor truth.append( { "gene": event.gene, "kind": event.kind, "indices": list(event.indices), "target_ids": [ assay.targets[index].target_id for index in event.indices ], "amplitude": event.amplitude, "factor": event.factor, "event_class": event.event_class, } ) else: truth.append({}) # Probe-site variants and sporadic capture failures are deliberately not # labelled as CNVs. A pure depth screen may flag them; in a clinical # workflow the reflex assay resolves them. artifact = rng.random(n_samples) < artifact_sample_fraction artifact_rows = np.flatnonzero(artifact) reportable = np.asarray( [ index for index, target in enumerate(assay.targets) if target.reportable ], dtype=int, ) for sample_index in artifact_rows: target = int(rng.choice(reportable)) factor = float(rng.choice((0.42, 0.55, 1.75))) rate[sample_index, target] *= factor truth[sample_index]["technical_artifact"] = { "target_id": assay.targets[target].target_id, "factor": factor, } return rng.poisson(rate).astype(np.int64), truth def simulate_normal_site_mixture( assay: Assay, n_samples: int, *, seed: int, site_shifts: Sequence[float], ) -> np.ndarray: """Generate a balanced QC-clean normal baseline across intended sites.""" if not site_shifts: raise ValueError("at least one intended site shift is required") base = n_samples // len(site_shifts) remainder = n_samples % len(site_shifts) chunks = [] for index, shift in enumerate(site_shifts): size = base + int(index < remainder) counts, _ = simulate_counts( assay, size, seed=seed + index, site_shift=float(shift), artifact_sample_fraction=0.0, ) chunks.append(counts) return np.vstack(chunks) def simulate_orthogonal_log2( assay: Assay, truth: Sequence[dict], *, seed: int, sigma: float = 0.055, independent_artifact_fraction: float = 0.002, ) -> np.ndarray: """Simulate an independently probed dosage confirmation assay. Biological copy-number events are shared with the primary assay. Primary capture artifacts are deliberately *not* copied; rare independent orthogonal artifacts are generated from a separate random stream. This is a workflow stress model, not a claim about MLPA or aCGH performance. """ rng = np.random.default_rng(seed) values = rng.normal(0.0, sigma, (len(truth), len(assay.targets))) for sample_index, event in enumerate(truth): if event.get("gene"): values[sample_index, event["indices"]] += math.log2( float(event["factor"]) ) reportable = np.asarray( [ index for index, target in enumerate(assay.targets) if target.reportable ], dtype=int, ) independent_artifact = ( rng.random(len(truth)) < independent_artifact_fraction ) for sample_index in np.flatnonzero(independent_artifact): target = int(rng.choice(reportable)) values[sample_index, target] += float( rng.choice((-0.55, 0.38)) ) return values def _matches_80pct_reciprocal(call, event: dict) -> bool: """Match against the complete event truth with 80% reciprocal overlap.""" if call.gene != event["gene"] or call.kind != event["kind"]: return False called = set(call.target_ids) actual = set(event["target_ids"]) overlap = called & actual return bool( called and actual and overlap and len(overlap) / len(called) >= 0.80 and len(overlap) / len(actual) >= 0.80 ) def _has_exact_endpoints(call, event: dict) -> bool: return bool( call.gene == event["gene"] and call.kind == event["kind"] and call.target_ids and event["target_ids"] and call.target_ids[0] == event["target_ids"][0] and call.target_ids[-1] == event["target_ids"][-1] ) def _has_exact_target_set(call, event: dict) -> bool: return bool( call.gene == event["gene"] and call.kind == event["kind"] and set(call.target_ids) == set(event["target_ids"]) ) def _independent_dosage_logic_retains( call, values: np.ndarray, target_index: dict[str, int], ) -> bool: observed = np.asarray( [values[target_index[target_id]] for target_id in call.target_ids], dtype=float, ) if observed.size == 0: return False required = max(1, int(math.ceil(observed.size * 0.5))) if call.kind == "DEL": return int(np.sum(observed < -0.25)) >= required return int(np.sum(observed > 0.20)) >= required def _wilson95(successes: int, total: int) -> list[float] | None: if total <= 0: return None z = 1.959963984540054 proportion = successes / total denominator = 1.0 + z * z / total centre = (proportion + z * z / (2.0 * total)) / denominator half_width = ( z * math.sqrt( proportion * (1.0 - proportion) / total + z * z / (4.0 * total * total) ) / denominator ) return [max(0.0, centre - half_width), min(1.0, centre + half_width)] def _endpoint_record(successes: int, total: int) -> dict: return { "successes": int(successes), "total": int(total), "rate": successes / total if total else None, "conditional_wilson_95": _wilson95(successes, total), } def _endpoint_by_population( intention_to_test_rows: Sequence[dict], completed_rows: Sequence[dict], field: str, ) -> dict: return { "intention_to_test": _endpoint_record( sum(bool(row[field]) for row in intention_to_test_rows), len(intention_to_test_rows), ), "completed": _endpoint_record( sum(bool(row[field]) for row in completed_rows), len(completed_rows), ), } def evaluate( model: AnchoredDepthModel, counts: np.ndarray, truth: Sequence[dict], *, prefix: str, orthogonal_log2: np.ndarray | None = None, ) -> dict: """Evaluate full synthetic truth without altering it for callability. ``completed`` positive-event denominators require both sample-QC completion and every truth target to be callable. Intention-to-test denominators keep every generated event. This makes target exclusion visible rather than converting an uncallable truth interval into a shorter, easier event. """ counts = np.asarray(counts) if counts.ndim != 2 or counts.shape[1] != len(model.targets): raise ValueError( "counts must be a two-dimensional sample-by-target matrix" ) if len(truth) != counts.shape[0]: raise ValueError("truth must contain one row per count row") callable_target_ids = { target.target_id for target, callable_target in zip( model.targets, model.callable_mask ) if target.reportable and callable_target } target_index = { target.target_id: index for index, target in enumerate(model.targets) } if orthogonal_log2 is not None: orthogonal_log2 = np.asarray(orthogonal_log2, dtype=float) if orthogonal_log2.shape != counts.shape: raise ValueError( "orthogonal_log2 must have the same dimensions as counts" ) rows = [] for index, (sample_counts, event) in enumerate(zip(counts, truth)): result = model.screen(sample_counts, sample=f"{prefix}_{index:05d}") dosage_retained_calls = ( tuple( call for call in result.calls if _independent_dosage_logic_retains( call, orthogonal_log2[index], target_index ) ) if orthogonal_log2 is not None else () ) if event.get("gene"): truth_target_ids = tuple(event["target_ids"]) truth_target_set = set(truth_target_ids) any_reflex_trigger = bool(result.calls) matched_80pct_reciprocal_overlap = any( _matches_80pct_reciprocal(call, event) for call in result.calls ) exact_endpoints = any( _has_exact_endpoints(call, event) for call in result.calls ) exact_target_set = any( _has_exact_target_set(call, event) for call in result.calls ) dosage_any_reflex_trigger = bool(dosage_retained_calls) dosage_matched_80pct_reciprocal_overlap = ( any( _matches_80pct_reciprocal(call, event) for call in dosage_retained_calls ) if orthogonal_log2 is not None else False ) dosage_exact_endpoints = ( any( _has_exact_endpoints(call, event) for call in dosage_retained_calls ) if orthogonal_log2 is not None else False ) dosage_exact_target_set = ( any( _has_exact_target_set(call, event) for call in dosage_retained_calls ) if orthogonal_log2 is not None else False ) n_truth_targets_callable = len( truth_target_set & callable_target_ids ) truth_callable_fraction = ( n_truth_targets_callable / len(truth_target_set) ) all_truth_targets_callable = ( n_truth_targets_callable == len(truth_target_set) ) any_truth_target_callable = n_truth_targets_callable > 0 else: any_reflex_trigger = bool(result.calls) matched_80pct_reciprocal_overlap = False exact_endpoints = False exact_target_set = False dosage_any_reflex_trigger = bool(dosage_retained_calls) dosage_matched_80pct_reciprocal_overlap = False dosage_exact_endpoints = False dosage_exact_target_set = False n_truth_targets_callable = 0 truth_callable_fraction = None all_truth_targets_callable = True any_truth_target_callable = True sample_qc_completed = result.status != "NO_CALL_QC" event_endpoint_completed = bool( sample_qc_completed and all_truth_targets_callable ) rows.append( { "event": event, "status": result.status, "sample_qc_completed": sample_qc_completed, "event_endpoint_completed": event_endpoint_completed, "any_reflex_trigger": any_reflex_trigger, "no_primary_referral": not any_reflex_trigger, "matched_80pct_reciprocal_overlap": matched_80pct_reciprocal_overlap, "exact_endpoints": exact_endpoints, "exact_target_set": exact_target_set, "dosage_any_reflex_trigger": dosage_any_reflex_trigger, "no_referral_after_dosage_logic": not dosage_any_reflex_trigger, "dosage_matched_80pct_reciprocal_overlap": dosage_matched_80pct_reciprocal_overlap, "dosage_exact_endpoints": dosage_exact_endpoints, "dosage_exact_target_set": dosage_exact_target_set, "n_calls": len(result.calls), "n_dosage_retained_calls": len(dosage_retained_calls), "n_truth_targets_callable": n_truth_targets_callable, "truth_callable_fraction": truth_callable_fraction, "all_truth_targets_callable": all_truth_targets_callable, "any_truth_target_callable": any_truth_target_callable, "median_control_depth": result.median_control_depth, "control_noise": result.control_noise, } ) positives = [row for row in rows if row["event"].get("gene")] negatives = [row for row in rows if not row["event"].get("gene")] completed_positive = [ row for row in positives if row["event_endpoint_completed"] ] completed_negative = [ row for row in negatives if row["sample_qc_completed"] ] clean_completed_negative = [ row for row in completed_negative if not row["event"].get("technical_artifact") ] artifact_completed_negative = [ row for row in completed_negative if row["event"].get("technical_artifact") ] summary = { "evaluation_label": "synthetic software verification", "truth_policy": ( "full event target IDs are retained and evaluated without " "intersection with the caller callable mask" ), "population_definitions": { "intention_to_test_positive": ( "every generated biological event, including sample-QC " "no-calls and events touching uncallable targets" ), "completed_positive": ( "sample passed caller QC and every full-truth event target " "was callable; truth is still not clipped" ), "intention_to_test_negative": ( "every generated copy-number-negative sample" ), "completed_negative": "copy-number-negative sample passed caller QC", }, "sample_counts": { "all": len(rows), "positive_events_intention_to_test": len(positives), "positive_events_completed": len(completed_positive), "positive_sample_qc_no_calls": sum( not row["sample_qc_completed"] for row in positives ), "positive_events_not_fully_callable": sum( not row["all_truth_targets_callable"] for row in positives ), "negative_samples_intention_to_test": len(negatives), "negative_samples_completed": len(completed_negative), "negative_sample_qc_no_calls": ( len(negatives) - len(completed_negative) ), }, "caller_callable_reportable_targets": { "n": len(callable_target_ids), "target_ids": sorted(callable_target_ids), }, "positive_truth_target_callability": { "all_truth_targets_callable": _endpoint_record( sum( row["all_truth_targets_callable"] for row in positives ), len(positives), ), "at_least_one_truth_target_callable": _endpoint_record( sum( row["any_truth_target_callable"] for row in positives ), len(positives), ), "mean_truth_target_callable_fraction": ( float( np.mean( [ row["truth_callable_fraction"] for row in positives ] ) ) if positives else None ), }, "positive_event_endpoints": { "any_reflex_trigger": _endpoint_by_population( positives, completed_positive, "any_reflex_trigger" ), "matched_80pct_reciprocal_overlap": _endpoint_by_population( positives, completed_positive, "matched_80pct_reciprocal_overlap", ), "exact_endpoints": _endpoint_by_population( positives, completed_positive, "exact_endpoints" ), "exact_target_set": _endpoint_by_population( positives, completed_positive, "exact_target_set" ), }, "negative_sample_referral_endpoints": { "any_primary_reflex_trigger": _endpoint_by_population( negatives, completed_negative, "any_reflex_trigger" ), "no_referral_rate_intention_to_test": _endpoint_record( sum(row["no_primary_referral"] for row in negatives), len(negatives), ), "no_referral_rate_completed": _endpoint_record( sum( row["no_primary_referral"] for row in completed_negative ), len(completed_negative), ), "primary_reflex_calls_in_completed_samples": int( sum(row["n_calls"] for row in completed_negative) ), "artifact_negative_subset": { "samples": len(artifact_completed_negative), "any_primary_reflex_trigger": _endpoint_record( sum( row["any_reflex_trigger"] for row in artifact_completed_negative ), len(artifact_completed_negative), ), }, "artifact_free_negative_subset": { "samples": len(clean_completed_negative), "any_primary_reflex_trigger": _endpoint_record( sum( row["any_reflex_trigger"] for row in clean_completed_negative ), len(clean_completed_negative), ), }, }, } strata_rows: dict[str, list[dict]] = {} for row in positives: event = row["event"] depth_bin = ( "lt150" if row["median_control_depth"] < 150 else "150to249" if row["median_control_depth"] < 250 else "ge250" ) keys = ( f"kind:{event['kind']}", f"class:{event['event_class']}", f"amplitude:{event['amplitude']}", f"depth:{depth_bin}", ) for key in keys: strata_rows.setdefault(key, []).append(row) summary["positive_event_strata"] = {} for key, bucket_rows in strata_rows.items(): completed_bucket = [ row for row in bucket_rows if row["event_endpoint_completed"] ] summary["positive_event_strata"][key] = { "sample_counts": { "intention_to_test": len(bucket_rows), "completed": len(completed_bucket), }, "any_reflex_trigger": _endpoint_by_population( bucket_rows, completed_bucket, "any_reflex_trigger" ), "matched_80pct_reciprocal_overlap": _endpoint_by_population( bucket_rows, completed_bucket, "matched_80pct_reciprocal_overlap", ), "exact_endpoints": _endpoint_by_population( bucket_rows, completed_bucket, "exact_endpoints" ), "exact_target_set": _endpoint_by_population( bucket_rows, completed_bucket, "exact_target_set" ), } if orthogonal_log2 is not None: summary["independent_dosage_channel_workflow_logic_only"] = { "scope": ( "software workflow logic using a second synthetic channel; " "not performance of MLPA, aCGH, ddPCR, or any clinical assay" ), "positive_event_endpoints_after_logic": { "any_reflex_trigger_retained": _endpoint_by_population( positives, completed_positive, "dosage_any_reflex_trigger", ), "matched_80pct_reciprocal_overlap_retained": _endpoint_by_population( positives, completed_positive, "dosage_matched_80pct_reciprocal_overlap", ), "exact_endpoints_retained": _endpoint_by_population( positives, completed_positive, "dosage_exact_endpoints", ), "exact_target_set_retained": _endpoint_by_population( positives, completed_positive, "dosage_exact_target_set", ), }, "negative_sample_referrals_retained_after_logic": _endpoint_by_population( negatives, completed_negative, "dosage_any_reflex_trigger", ), "negative_sample_no_referral_rate_after_logic_intention_to_test": _endpoint_record( sum( row["no_referral_after_dosage_logic"] for row in negatives ), len(negatives), ), "negative_sample_no_referral_rate_after_logic_completed": _endpoint_record( sum( row["no_referral_after_dosage_logic"] for row in completed_negative ), len(completed_negative), ), "retained_calls_in_completed_negative_samples": int( sum( row["n_dosage_retained_calls"] for row in completed_negative ) ), } return summary def _development_utility( result: dict, callable_reportable_targets: int, total_reportable_targets: int, ) -> tuple: """Engineering-only selection: gates first, then fewer referrals.""" overlap_rate = ( result["positive_event_endpoints"][ "matched_80pct_reciprocal_overlap" ]["completed"]["rate"] or 0.0 ) negative_referral_rate = ( result["negative_sample_referral_endpoints"][ "any_primary_reflex_trigger" ]["completed"]["rate"] or 0.0 ) callable_fraction = ( callable_reportable_targets / total_reportable_targets if total_reportable_targets else 0.0 ) positive_completion = ( result["sample_counts"]["positive_events_completed"] / result["sample_counts"]["positive_events_intention_to_test"] if result["sample_counts"]["positive_events_intention_to_test"] else 0.0 ) gates = ( int(overlap_rate >= 0.98) + int(negative_referral_rate <= 0.025) + int(callable_fraction >= 0.90) + int(positive_completion >= 0.85) ) return ( gates, callable_fraction, positive_completion, overlap_rate, -negative_referral_rate, ) def _headline_findings(result: dict) -> dict: positive = result["positive_event_endpoints"] negative = result["negative_sample_referral_endpoints"] counts = result["sample_counts"] return { "positive_events_intention_to_test": counts["positive_events_intention_to_test"], "positive_events_completed": counts["positive_events_completed"], "positive_sample_qc_no_calls": counts["positive_sample_qc_no_calls"], "positive_events_not_fully_callable": counts["positive_events_not_fully_callable"], "positive_completion_rate": ( counts["positive_events_completed"] / counts["positive_events_intention_to_test"] if counts["positive_events_intention_to_test"] else None ), "any_reflex_trigger_rate_intention_to_test": positive["any_reflex_trigger"]["intention_to_test"]["rate"], "matched_80pct_reciprocal_overlap_rate_intention_to_test": positive["matched_80pct_reciprocal_overlap"][ "intention_to_test" ]["rate"], "matched_80pct_reciprocal_overlap_rate_completed": positive["matched_80pct_reciprocal_overlap"]["completed"][ "rate" ], "exact_endpoints_rate_intention_to_test": positive["exact_endpoints"]["intention_to_test"]["rate"], "exact_target_set_rate_intention_to_test": positive["exact_target_set"]["intention_to_test"]["rate"], "all_truth_targets_callable_rate": result["positive_truth_target_callability"][ "all_truth_targets_callable" ]["rate"], "negative_samples_intention_to_test": counts["negative_samples_intention_to_test"], "negative_samples_completed": counts["negative_samples_completed"], "negative_sample_qc_no_calls": counts["negative_sample_qc_no_calls"], "negative_primary_referral_rate_completed": negative["any_primary_reflex_trigger"]["completed"]["rate"], "negative_no_referral_rate_completed": negative["no_referral_rate_completed"]["rate"], } def run_experiment( reference_dir: Path, output_path: Path, model_path: Path, *, quick: bool = False, verification_seed_offset: int = 300, ) -> dict: started = time.perf_counter() panel = Panel.load(str(reference_dir)) assay = build_assay(panel) sizes = ( { "fit": 250, "calibration": 400, "development_normal": 300, "development_positive": 500, "locked_same_assay_normal": 400, "locked_same_assay_positive": 700, "stress_normal": 250, "stress_positive": 400, "transport_normal": 300, "transport_positive": 500, } if quick else { "fit": 1300, "calibration": 3000, "development_normal": 1200, "development_positive": 2400, "locked_same_assay_normal": 2000, "locked_same_assay_positive": 4000, "stress_normal": 1200, "stress_positive": 2400, "transport_normal": 1200, "transport_positive": 2400, } ) # The locked model is assay/site/chemistry specific. Fit, calibrate, # develop, and perform same-assay verification at one intended synthetic # site. Low-depth/site/mosaic stress and target-specific assay transport are # separate evaluations and are deliberately reported first in the artifact. deployment_site_shift = 0.55 intended_site_shifts = (deployment_site_shift,) random_streams = { "fit_normals": BASE_SEED + 10, "calibration_normals": BASE_SEED + 11, "development_normals": BASE_SEED + 12, "development_event_design": BASE_SEED + 13, "development_positives": BASE_SEED + 14, "locked_same_assay_normals": BASE_SEED + verification_seed_offset, "locked_same_assay_event_design": BASE_SEED + verification_seed_offset + 1, "locked_same_assay_positives": BASE_SEED + verification_seed_offset + 2, "locked_same_assay_independent_dosage": BASE_SEED + verification_seed_offset + 3, "within_assay_stress_normals": BASE_SEED + verification_seed_offset + 99, "within_assay_stress_event_design": BASE_SEED + verification_seed_offset + 100, "within_assay_stress_positives": BASE_SEED + verification_seed_offset + 101, "within_assay_stress_independent_dosage": BASE_SEED + verification_seed_offset + 102, "transport_assay_construction": BASE_SEED + verification_seed_offset + 200, "transport_normals": BASE_SEED + verification_seed_offset + 201, "transport_event_design": BASE_SEED + verification_seed_offset + 202, "transport_positives": BASE_SEED + verification_seed_offset + 203, "transport_independent_dosage": BASE_SEED + verification_seed_offset + 204, } fit_counts = simulate_normal_site_mixture( assay, sizes["fit"], seed=random_streams["fit_normals"], site_shifts=intended_site_shifts, ) calibration_counts = simulate_normal_site_mixture( assay, sizes["calibration"], seed=random_streams["calibration_normals"], site_shifts=intended_site_shifts, ) fit_sample_ids = tuple( f"synthetic-fit-{index:06d}" for index in range(sizes["fit"]) ) calibration_sample_ids = tuple( f"synthetic-calibration-{index:06d}" for index in range(sizes["calibration"]) ) if ( len(set(fit_sample_ids)) != len(fit_sample_ids) or len(set(calibration_sample_ids)) != len(calibration_sample_ids) or set(fit_sample_ids).intersection(calibration_sample_ids) ): raise RuntimeError("synthetic fit/calibration IDs are not disjoint") development_normal, development_normal_truth = simulate_counts( assay, sizes["development_normal"], seed=random_streams["development_normals"], site_shift=deployment_site_shift, artifact_sample_fraction=0.015, ) development_events = event_design( assay, n=sizes["development_positive"], seed=random_streams["development_event_design"], ) development_positive, development_positive_truth = simulate_counts( assay, sizes["development_positive"], seed=random_streams["development_positives"], events=development_events, site_shift=deployment_site_shift, artifact_sample_fraction=0.0, ) development_counts = np.vstack( (development_normal, development_positive) ) development_truth = ( development_normal_truth + development_positive_truth ) variants: list[dict] = [] models: dict[tuple[int, float], AnchoredDepthModel] = {} total_reportable = sum(target.reportable for target in assay.targets) for n_components in COMPONENT_GRID: for alpha in ALPHA_GRID: model = AnchoredDepthModel.fit( fit_counts, calibration_counts, assay.targets, alpha=alpha, n_components=n_components, min_median_control_depth=100.0, min_target_median_depth=100.0, max_target_sigma=0.13, fit_sample_ids=fit_sample_ids, calibration_sample_ids=calibration_sample_ids, ) result = evaluate( model, development_counts, development_truth, prefix=f"development_k{n_components}_a{alpha}", ) callable_reportable = int( sum( bool(callable_target and target.reportable) for callable_target, target in zip( model.callable_mask, model.targets ) ) ) variants.append( { "n_components": n_components, "alpha": alpha, "threshold": model.family_threshold, "noise_threshold": model.noise_threshold, "callable_reportable_targets": callable_reportable, "result": result, "selection_utility": list( _development_utility( result, callable_reportable, total_reportable ) ), } ) models[(n_components, alpha)] = model # Utility is evaluated only on the development cohort. Ties deliberately # select the smaller model. selected = max( variants, key=lambda row: ( tuple(row["selection_utility"]), -int(row["n_components"]), -float(row["alpha"]), ), ) chosen_components = int(selected["n_components"]) chosen_alpha = float(selected["alpha"]) model = models[(chosen_components, chosen_alpha)] locked_normal, locked_normal_truth = simulate_counts( assay, sizes["locked_same_assay_normal"], seed=random_streams["locked_same_assay_normals"], site_shift=deployment_site_shift, artifact_sample_fraction=0.015, ) locked_events = event_design( assay, n=sizes["locked_same_assay_positive"], seed=random_streams["locked_same_assay_event_design"], ) locked_positive, locked_positive_truth = simulate_counts( assay, sizes["locked_same_assay_positive"], seed=random_streams["locked_same_assay_positives"], events=locked_events, site_shift=deployment_site_shift, artifact_sample_fraction=0.0, ) locked_truth = locked_normal_truth + locked_positive_truth locked_independent_dosage = simulate_orthogonal_log2( assay, locked_truth, seed=random_streams["locked_same_assay_independent_dosage"], ) locked_same_assay = evaluate( model, np.vstack((locked_normal, locked_positive)), locked_truth, prefix="locked_same_assay_verification", orthogonal_log2=locked_independent_dosage, ) stress_normal, stress_normal_truth = simulate_counts( assay, sizes["stress_normal"], seed=random_streams["within_assay_stress_normals"], site_shift=0.85, median_depth=280.0, low_depth_fraction=0.30, artifact_sample_fraction=0.015, ) stress_events = event_design( assay, n=sizes["stress_positive"], seed=random_streams["within_assay_stress_event_design"], include_mosaic=True, ) stress_positive, stress_positive_truth = simulate_counts( assay, sizes["stress_positive"], seed=random_streams["within_assay_stress_positives"], events=stress_events, site_shift=0.85, median_depth=280.0, low_depth_fraction=0.30, artifact_sample_fraction=0.0, ) stress_truth = stress_normal_truth + stress_positive_truth stress_independent_dosage = simulate_orthogonal_log2( assay, stress_truth, seed=random_streams["within_assay_stress_independent_dosage"], ) stress = evaluate( model, np.vstack((stress_normal, stress_positive)), stress_truth, prefix="within_assay_stress", orthogonal_log2=stress_independent_dosage, ) transport_assay, transport_shift = build_transport_shift_assay( assay, seed=random_streams["transport_assay_construction"] ) if tuple(target.target_id for target in transport_assay.targets) != tuple( target.target_id for target in model.targets ): raise RuntimeError("transport assay target IDs are not aligned") transport_normal, transport_normal_truth = simulate_counts( transport_assay, sizes["transport_normal"], seed=random_streams["transport_normals"], site_shift=deployment_site_shift, artifact_sample_fraction=0.015, ) transport_events = event_design( transport_assay, n=sizes["transport_positive"], seed=random_streams["transport_event_design"], ) transport_positive, transport_positive_truth = simulate_counts( transport_assay, sizes["transport_positive"], seed=random_streams["transport_positives"], events=transport_events, site_shift=deployment_site_shift, artifact_sample_fraction=0.0, ) transport_truth = ( transport_normal_truth + transport_positive_truth ) transport_independent_dosage = simulate_orthogonal_log2( transport_assay, transport_truth, seed=random_streams["transport_independent_dosage"], ) transport = evaluate( model, np.vstack((transport_normal, transport_positive)), transport_truth, prefix="transport_shift", orthogonal_log2=transport_independent_dosage, ) model.save(model_path) reference_files = { name: { "path": str((reference_dir / name).resolve()), "sha256": sha256_file(reference_dir / name), } for name in ("panel.json", "panel.fa") } script_path = Path(__file__).resolve() callable_reportable_target_ids = [ target.target_id for target, callable_target in zip( model.targets, model.callable_mask ) if target.reportable and callable_target ] payload = { "schema": SCHEMA, "artifact_type": "synthetic software verification", "generated_utc": datetime.now(timezone.utc).isoformat().replace( "+00:00", "Z" ), "runtime_seconds": round(time.perf_counter() - started, 3), "claims_boundary": [ "Every endpoint is conditional synthetic software verification, " "not analytical or clinical performance.", "Conditional Wilson intervals describe finite simulated draws only " "and exclude simulator, assay, and transport uncertainty.", "Full event truth is never clipped to the caller callable mask; " "truth-target callability is reported separately.", "The locked same-assay stream was not used to select model " "components or thresholds.", "The independent dosage channel tests workflow logic only; it is " "not evidence about MLPA, aCGH, ddPCR, or clinical confirmation.", "Stress and transport-shift findings must be reviewed before the " "more favorable same-generator results.", ], "stress_and_transport_findings_review_first": { "reason": ( "distribution shift is the principal limitation of a " "same-generator simulation" ), "transport_shift": _headline_findings(transport), "within_assay_low_depth_site_mosaic_stress": _headline_findings(stress), }, "design": { "base_seed": BASE_SEED, "quick": quick, "sizes": sizes, "n_targets": len(assay.targets), "n_control_targets": sum( not target.reportable for target in assay.targets ), "n_reportable_targets": sum( target.reportable for target in assay.targets ), "unstable_reportable_targets": [ assay.targets[index].target_id for index in assay.unstable_reportable ], "component_grid": list(COMPONENT_GRID), "family_alpha_grid": list(ALPHA_GRID), "selected_family_alpha": chosen_alpha, "selection_rule": ( "engineering-only: maximize development gates (completed " "80%-reciprocal-overlap rate>=0.98, completed negative " "referral rate<=0.025, >=90% targets callable, >=85% positive " "completion), then callable fraction, completion, overlap " "rate, fewer negative referrals, and fewer components" ), "selection_gates_are_clinical_acceptance_criteria": False, "deployment_site_shift": deployment_site_shift, "locked_same_assay_site_shift": deployment_site_shift, "verification_seed_offset": verification_seed_offset, "intended_site_shifts_in_fit_and_calibration": list( intended_site_shifts ), "transport_shift": transport_shift, "normal_artifact_fraction": 0.015, "independent_dosage_channel_workflow_logic_only": { "independent_log2_sigma": 0.055, "independent_artifact_fraction": 0.002, "deletion_threshold": -0.25, "duplication_threshold": 0.20, "minimum_fraction_of_called_targets": 0.5, "scope": ( "synthetic workflow verification only; not a claim about " "MLPA or aCGH clinical performance" ), }, }, "transport_shift_synthetic_software_verification": transport, "within_assay_stress_synthetic_software_verification": stress, "locked_same_assay_synthetic_software_verification": locked_same_assay, "development_model_selection": { "variants": variants, "selected_components": chosen_components, "selected_alpha": chosen_alpha, }, "selected_model": { "path": str(model_path.resolve()), "sha256": sha256_file(model_path), "schema": model.to_json()["schema"], "family_threshold": model.family_threshold, "noise_threshold": model.noise_threshold, "callable_reportable_targets": { "n": len(callable_reportable_target_ids), "target_ids": callable_reportable_target_ids, }, "n_fit_synthetic_normals": model.n_fit, "n_calibration_synthetic_normals": model.n_calibration, }, "provenance": { "source_script": { "path": str(script_path), "sha256": sha256_file(script_path), }, "reference_files": reference_files, "numpy_random_generator": "default_rng/PCG64", "random_streams": random_streams, "fit_and_calibration_cohorts": { "unique_and_disjoint_sample_ids": True, "fit_sample_id_pattern": "synthetic-fit-{index:06d}", "calibration_sample_id_pattern": "synthetic-calibration-{index:06d}", "fit_sample_ids_sha256": model.fit_sample_ids_sha256, "calibration_sample_ids_sha256": model.calibration_sample_ids_sha256, "fit_count_rows_sha256": model.fit_count_rows_sha256, "calibration_count_rows_sha256": model.calibration_count_rows_sha256, }, "truth_hashes": { "development_full_truth_sha256": sha256_json(development_truth), "locked_same_assay_full_truth_sha256": sha256_json(locked_truth), "within_assay_stress_full_truth_sha256": sha256_json(stress_truth), "transport_shift_full_truth_sha256": sha256_json(transport_truth), }, }, "environment": { "python": platform.python_version(), "numpy": np.__version__, "sklearn": sklearn.__version__, }, } atomic_write_json(str(output_path), payload) return payload def main() -> int: parser = argparse.ArgumentParser() parser.add_argument( "--reference-dir", type=Path, default=Path("reference_grch37") ) parser.add_argument( "--output", type=Path, default=Path("analysis/clinical_panel_simulation.json"), ) parser.add_argument( "--model", type=Path, default=Path("analysis/anchored_clinical_depth_model.json"), ) parser.add_argument("--quick", action="store_true") parser.add_argument( "--verification-seed-offset", type=int, default=300 ) args = parser.parse_args() payload = run_experiment( args.reference_dir, args.output, args.model, quick=args.quick, verification_seed_offset=args.verification_seed_offset, ) print( json.dumps( { "artifact_type": payload["artifact_type"], "selected_model": { "n_components": payload[ "development_model_selection" ]["selected_components"], "family_alpha": payload[ "development_model_selection" ]["selected_alpha"], }, "stress_and_transport_findings_review_first": payload["stress_and_transport_findings_review_first"], "locked_same_assay_synthetic_software_verification": _headline_findings( payload[ "locked_same_assay_synthetic_software_verification" ] ), "runtime_seconds": payload["runtime_seconds"], }, indent=2, ) ) return 0 if __name__ == "__main__": raise SystemExit(main())