"""Clinically oriented, control-anchored read-depth screening. This module implements the part of an exon-CNV workflow that decides which specimens need orthogonal dosage confirmation. It is intentionally a *screen*, not an autonomous clinical sign-out: * normalization is anchored on targets outside the reportable genes; * latent-factor coefficients are estimated from those controls, so a BRCA event cannot normalize itself away; * target dispersion, callability, weights, and QC limits are learned only from a fit cohort; an untouched calibration cohort sets one joint DEL/DUP, all-window threshold; * low-depth/noisy specimens and unstable targets are no-calls; * every emitted candidate is labelled ``REFLEX_POSITIVE``. The architecture follows the common clinical pattern used by hereditary-cancer panels: broad-panel normalization/PCA, exon-bin segmentation, explicit QC, and confirmation of NGS CNV flags by MLPA, aCGH, or another independently designed dosage assay. Internal calibration or simulation must never be described as clinical validation. """ from __future__ import annotations import hashlib import json import math from dataclasses import asdict, dataclass from numbers import Integral, Real from pathlib import Path from typing import Sequence import numpy as np from sklearn.covariance import OAS from .provenance import atomic_write_json SCHEMA = "lgasieve.anchored-clinical-depth-screen.v2" LEGACY_SCHEMA = "lgasieve.anchored-clinical-depth-screen.v1" _EPS = 1e-9 CLAIMS_BOUNDARY = ( "Research/development reflex screen only. A candidate requires " "confirmation by an independently designed dosage assay. No reflex flag " "does not establish absence of an LGA and must not be used for clinical " "exclusion." ) def _robust_sigma(values: np.ndarray, axis: int | None = None) -> np.ndarray: values = np.asarray(values, dtype=float) centre = np.median(values, axis=axis, keepdims=True) mad = np.median(np.abs(values - centre), axis=axis) return 1.4826 * mad def _hash_strings(values: Sequence[str]) -> str: encoded = json.dumps( list(values), ensure_ascii=False, separators=(",", ":") ).encode("utf-8") return hashlib.sha256(encoded).hexdigest() def _row_fingerprints(matrix: np.ndarray) -> tuple[str, ...]: canonical = np.asarray(matrix, dtype=" tuple[str, ...] | None: if sample_ids is None: return None if isinstance(sample_ids, (str, bytes)): raise ValueError(f"{name} sample IDs must be a sequence of IDs") values = tuple(str(value).strip() for value in sample_ids) if len(values) != expected: raise ValueError(f"{name} sample IDs do not match count rows") if any(not value for value in values): raise ValueError(f"{name} sample IDs must be non-empty") if len(set(values)) != len(values): raise ValueError(f"{name} sample IDs must be unique") return values def _window_index_family( targets: Sequence["DepthTarget"], callable_mask: np.ndarray, max_window_targets: int, ) -> tuple[tuple[str, tuple[int, ...]], ...]: """Return every window in uninterrupted callable reportable runs.""" by_gene: dict[str, list[int]] = {} for index, target in enumerate(targets): by_gene.setdefault(target.gene, []).append(index) windows: list[tuple[str, tuple[int, ...]]] = [] for gene, indices in by_gene.items(): ordered = sorted(indices, key=lambda index: targets[index].order) runs: list[tuple[int, ...]] = [] run: list[int] = [] for index in ordered: if targets[index].reportable and callable_mask[index]: run.append(index) elif run: runs.append(tuple(run)) run = [] if run: runs.append(tuple(run)) for callable_run in runs: width_limit = min(len(callable_run), int(max_window_targets)) for width in range(1, width_limit + 1): for start in range(len(callable_run) - width + 1): windows.append( (gene, callable_run[start : start + width]) ) return tuple(windows) def _valid_sha256(value, *, optional: bool, name: str) -> str | None: if value is None and optional: return None if not isinstance(value, str): raise ValueError(f"{name} must be a lowercase SHA-256 digest") if len(value) != 64 or any( character not in "0123456789abcdef" for character in value ): raise ValueError(f"{name} must be a lowercase SHA-256 digest") return value @dataclass(frozen=True) class DepthTarget: """One analysis bin. ``reportable`` targets are scanned for CNVs. All remaining targets form the invariant normalization domain. Multiple bins may share one exon. """ target_id: str gene: str exon: str order: int reportable: bool @dataclass(frozen=True) class DepthCall: gene: str kind: str first_exon: str last_exon: str target_ids: tuple[str, ...] n_targets: int mean_log2: float statistic: float threshold: float disposition: str = "REFLEX_POSITIVE" @dataclass(frozen=True) class DepthScreenResult: sample: str status: str calls: tuple[DepthCall, ...] no_call_targets: tuple[str, ...] median_control_depth: float control_noise: float reasons: tuple[str, ...] claims_boundary: str = CLAIMS_BOUNDARY def to_json(self) -> dict: return { "sample": self.sample, "status": self.status, "calls": [asdict(call) for call in self.calls], "no_call_targets": list(self.no_call_targets), "median_control_depth": self.median_control_depth, "control_noise": self.control_noise, "reasons": list(self.reasons), "claims_boundary": self.claims_boundary, } @dataclass(frozen=True) class _Window: gene: str indices: tuple[int, ...] target_ids: tuple[str, ...] first_exon: str last_exon: str weights: tuple[float, ...] precision_sum: float class AnchoredDepthModel: """Broad-panel PCA/read-depth model with independent null calibration.""" def __init__(self) -> None: self.targets: tuple[DepthTarget, ...] = () self.control_mask = np.zeros(0, dtype=bool) self.centre = np.zeros(0, dtype=float) self.components = np.zeros((0, 0), dtype=float) self.residual_centre = np.zeros(0, dtype=float) self.residual_sigma = np.zeros(0, dtype=float) self.callable_mask = np.zeros(0, dtype=bool) self.windows: tuple[_Window, ...] = () self.family_threshold = math.inf self.noise_threshold = math.inf self.min_median_control_depth = 100.0 self.min_target_median_depth = 100.0 self.max_target_sigma = 0.22 self.max_window_targets = 60 self.alpha = 0.02 self.n_fit = 0 self.n_calibration = 0 self.n_components = 0 self.requested_n_components = 0 self.calibration_family_max = np.zeros(0, dtype=float) self.fit_sample_ids_sha256: str | None = None self.calibration_sample_ids_sha256: str | None = None self.fit_count_rows_sha256 = "" self.calibration_count_rows_sha256 = "" @staticmethod def _validate_counts( counts: np.ndarray, n_targets: int, name: str, ) -> np.ndarray: matrix = np.asarray(counts, dtype=float) if matrix.ndim != 2 or matrix.shape[1] != n_targets: raise ValueError( f"{name} must be a two-dimensional matrix with " f"{n_targets} target columns" ) if matrix.shape[0] < 1 or not np.all(np.isfinite(matrix)): raise ValueError(f"{name} must contain finite sample rows") if np.any(matrix < 0): raise ValueError(f"{name} cannot contain negative counts") return matrix @staticmethod def _library_log2(counts: np.ndarray, control_mask: np.ndarray) -> np.ndarray: values = np.asarray(counts, dtype=float) one_row = values.ndim == 1 values = np.atleast_2d(values) scale = np.median(values[:, control_mask] + 0.5, axis=1) logged = np.log2(values + 0.5) - np.log2(np.maximum(scale, _EPS))[:, None] return logged[0] if one_row else logged def _residual_one(self, counts: np.ndarray) -> np.ndarray: logged = self._library_log2(counts, self.control_mask) delta = logged - self.centre if self.n_components: design = self.components[:, self.control_mask].T response = delta[self.control_mask] weights = np.ones(response.size, dtype=float) coef = np.zeros(self.n_components, dtype=float) for _ in range(4): weighted_design = design * weights[:, None] gram = design.T @ weighted_design gram += 1e-6 * np.eye(self.n_components) coef = np.linalg.solve(gram, weighted_design.T @ response) error = response - design @ coef scale = max(float(_robust_sigma(error)), 0.03) weights = np.minimum(1.0, 2.5 * scale / np.maximum(np.abs(error), _EPS)) delta = delta - coef @ self.components return delta def residuals(self, counts: np.ndarray) -> np.ndarray: matrix = self._validate_counts(counts, len(self.targets), "counts") return np.vstack([self._residual_one(row) for row in matrix]) @staticmethod def _conformal_threshold(values: np.ndarray, alpha: float) -> float: ordered = np.sort(np.asarray(values, dtype=float)) # Smallest observed order statistic whose conformal upper-tail # probability is no greater than alpha. If calibration is too small # for the requested alpha, the maximum is deliberately returned. rank = int(math.ceil((ordered.size + 1) * (1.0 - alpha))) - 1 rank = min(max(rank, 0), ordered.size - 1) return float(ordered[rank]) @classmethod def fit( cls, fit_counts: np.ndarray, calibration_counts: np.ndarray, targets: Sequence[DepthTarget], *, alpha: float = 0.02, n_components: int = 8, min_median_control_depth: float = 100.0, min_target_median_depth: float = 100.0, max_target_sigma: float = 0.22, max_window_targets: int = 60, fit_sample_ids: Sequence[str] | None = None, calibration_sample_ids: Sequence[str] | None = None, ) -> "AnchoredDepthModel": """Fit all nuisance structure, then calibrate one frozen search family. ``fit_counts`` and ``calibration_counts`` must be disjoint normal cohorts. Centre, dispersion, target callability, window membership, covariance weights, and sample-QC limits are learned exclusively from ``fit_counts``. The untouched calibration cohort only supplies one empirical maximum across every DEL/DUP window. """ target_rows = tuple(targets) if len(target_rows) < 3: raise ValueError("at least three targets are required") if not all(isinstance(target, DepthTarget) for target in target_rows): raise ValueError("targets must be DepthTarget objects") if any( not isinstance(target.target_id, str) or not target.target_id.strip() or not isinstance(target.gene, str) or not target.gene.strip() or not isinstance(target.exon, str) or not target.exon.strip() for target in target_rows ): raise ValueError("target IDs, genes, and exons must be non-empty") if any( isinstance(target.order, bool) or not isinstance(target.order, Integral) for target in target_rows ): raise ValueError("target order must be an integer") if any( not isinstance(target.reportable, (bool, np.bool_)) for target in target_rows ): raise ValueError("target reportable values must be booleans") if len({target.target_id for target in target_rows}) != len( target_rows ): raise ValueError("target IDs must be unique") gene_orders = [(target.gene, target.order) for target in target_rows] if len(set(gene_orders)) != len(gene_orders): raise ValueError("target order must be unique within each gene") control_mask = np.asarray( [not target.reportable for target in target_rows], dtype=bool ) reportable_mask = ~control_mask if control_mask.sum() < 20: raise ValueError( "at least 20 non-reportable control targets are required; " "a BRCA-only matrix has no invariant normalization anchor" ) if reportable_mask.sum() < 1: raise ValueError("at least one reportable target is required") if ( isinstance(alpha, bool) or not isinstance(alpha, Real) or not np.isfinite(alpha) or not 0.0 < float(alpha) < 0.5 ): raise ValueError("alpha must be between zero and 0.5") alpha = float(alpha) if ( isinstance(n_components, bool) or not isinstance(n_components, Integral) or int(n_components) < 0 ): raise ValueError("n_components must be a non-negative integer") n_components = int(n_components) if ( isinstance(min_median_control_depth, bool) or not isinstance(min_median_control_depth, Real) or not np.isfinite(min_median_control_depth) or min_median_control_depth <= 0 ): raise ValueError( "min_median_control_depth must be finite and positive" ) if ( isinstance(min_target_median_depth, bool) or not isinstance(min_target_median_depth, Real) or not np.isfinite(min_target_median_depth) or min_target_median_depth < 0 ): raise ValueError( "min_target_median_depth must be finite and non-negative" ) if ( isinstance(max_target_sigma, bool) or not isinstance(max_target_sigma, Real) or not np.isfinite(max_target_sigma) or max_target_sigma <= 0 ): raise ValueError("max_target_sigma must be finite and positive") if ( isinstance(max_window_targets, bool) or not isinstance(max_window_targets, Integral) or int(max_window_targets) < 1 ): raise ValueError( "max_window_targets must be a positive integer" ) fit_matrix = cls._validate_counts( fit_counts, len(target_rows), "fit_counts" ) calibration_matrix = cls._validate_counts( calibration_counts, len(target_rows), "calibration_counts" ) if fit_matrix.shape[0] < 50: raise ValueError("at least 50 fit normals are required") minimum_calibration = math.ceil(1.0 / alpha) - 1 if calibration_matrix.shape[0] < minimum_calibration: raise ValueError( "calibration cohort cannot attain the requested empirical " f"family alpha; need at least {minimum_calibration} libraries" ) fit_ids = _normalise_sample_ids( fit_sample_ids, fit_matrix.shape[0], "fit" ) calibration_ids = _normalise_sample_ids( calibration_sample_ids, calibration_matrix.shape[0], "calibration", ) if ( fit_ids is not None and calibration_ids is not None and set(fit_ids).intersection(calibration_ids) ): raise ValueError( "fit and calibration sample IDs must be disjoint" ) fit_fingerprints = _row_fingerprints(fit_matrix) calibration_fingerprints = _row_fingerprints( calibration_matrix ) if set(fit_fingerprints).intersection( calibration_fingerprints ): raise ValueError( "fit and calibration cohorts contain identical count rows" ) out = cls() out.targets = target_rows out.control_mask = control_mask out.alpha = float(alpha) out.min_median_control_depth = float(min_median_control_depth) out.min_target_median_depth = float(min_target_median_depth) out.max_target_sigma = float(max_target_sigma) out.max_window_targets = int(max_window_targets) out.n_fit = int(fit_matrix.shape[0]) out.n_calibration = int(calibration_matrix.shape[0]) out.requested_n_components = int(n_components) out.fit_sample_ids_sha256 = ( _hash_strings(fit_ids) if fit_ids is not None else None ) out.calibration_sample_ids_sha256 = ( _hash_strings(calibration_ids) if calibration_ids is not None else None ) out.fit_count_rows_sha256 = _hash_strings( sorted(fit_fingerprints) ) out.calibration_count_rows_sha256 = _hash_strings( sorted(calibration_fingerprints) ) logged = cls._library_log2(fit_matrix, control_mask) out.centre = np.median(logged, axis=0) centred = logged - out.centre max_components = min( int(n_components), fit_matrix.shape[0] - 2, int(control_mask.sum()) - 1, len(target_rows) - 1, ) if max_components: _, _, vt = np.linalg.svd(centred, full_matrices=False) out.components = vt[:max_components] else: out.components = np.zeros((0, len(target_rows)), dtype=float) out.n_components = int(max_components) # Everything below that defines a test statistic is learned only from # the fit cohort. Calibration libraries must never fit themselves. fit_residual = out.residuals(fit_matrix) out.residual_centre = np.median(fit_residual, axis=0) centred_fit_residual = fit_residual - out.residual_centre out.residual_sigma = np.maximum( _robust_sigma(centred_fit_residual, axis=0), 0.03 ) target_depth = np.median(fit_matrix, axis=0) out.callable_mask = ( reportable_mask & (target_depth >= float(min_target_median_depth)) & (out.residual_sigma <= float(max_target_sigma)) ) # Sample-level control scatter is a QC metric, not a multiplicative # correction to event scores. The 99.5th-percentile rule is learned # only on fit normals so calibration remains untouched. control_noise = _robust_sigma( centred_fit_residual[:, control_mask], axis=1 ) out.noise_threshold = float( np.quantile(control_noise, 0.995, method="higher") ) windows: list[_Window] = [] window_family = _window_index_family( target_rows, out.callable_mask, out.max_window_targets ) for gene, indices_tuple in window_family: width = len(indices_tuple) sub = centred_fit_residual[:, indices_tuple] covariance = OAS(store_precision=False).fit(sub).covariance_ covariance = np.atleast_2d(covariance) precision_one = np.linalg.solve( covariance + 1e-6 * np.eye(width), np.ones(width, dtype=float), ) precision_sum = float(np.sum(precision_one)) selected = [target_rows[index] for index in indices_tuple] windows.append( _Window( gene=gene, indices=indices_tuple, target_ids=tuple(row.target_id for row in selected), first_exon=selected[0].exon, last_exon=selected[-1].exon, weights=tuple( float(value) for value in precision_one ), precision_sum=precision_sum, ) ) if not windows: raise ValueError("no contiguous callable reportable window remains") out.windows = tuple(windows) calibration_residual = ( out.residuals(calibration_matrix) - out.residual_centre ) out.calibration_family_max = out._family_statistics( calibration_residual ) out.family_threshold = out._conformal_threshold( out.calibration_family_max, out.alpha ) return out @staticmethod def _directional_statistic( projection: np.ndarray, precision_sum: float, mean_shift: float, ) -> np.ndarray: # Generalized least-squares matched-filter z. The magnitude is learned # from the data rather than fixed at exactly two or three copies, which # preserves sensitivity to mosaic/impure material. ``mean_shift`` is # used only for the predeclared direction. direction = -1.0 if mean_shift < 0 else 1.0 z = direction * np.asarray(projection, dtype=float) / math.sqrt( max(float(precision_sum), _EPS) ) return np.maximum(z, 0.0) def _family_statistics(self, centred_residual: np.ndarray) -> np.ndarray: residual = np.atleast_2d(np.asarray(centred_residual, dtype=float)) family_max = np.zeros(residual.shape[0], dtype=float) for window in self.windows: weights = np.asarray(window.weights, dtype=float) projection = residual[:, window.indices] @ weights deletion = self._directional_statistic( projection, window.precision_sum, -1.0 ) duplication = self._directional_statistic( projection, window.precision_sum, math.log2(1.5) ) family_max = np.maximum( family_max, np.maximum(deletion, duplication) ) return family_max def screen(self, counts: np.ndarray, sample: str = "sample") -> DepthScreenResult: values = np.asarray(counts, dtype=float) if values.ndim != 1 or values.size != len(self.targets): raise ValueError( f"counts must be one vector with {len(self.targets)} values" ) if not np.all(np.isfinite(values)) or np.any(values < 0): raise ValueError("counts must be finite and nonnegative") median_control_depth = float(np.median(values[self.control_mask])) raw_residual = self._residual_one(values) residual = raw_residual - self.residual_centre control_noise = float(_robust_sigma(residual[self.control_mask])) reasons: list[str] = [] if median_control_depth < self.min_median_control_depth: reasons.append( f"median control depth {median_control_depth:.1f} is below " f"{self.min_median_control_depth:.1f}" ) if control_noise > self.noise_threshold: reasons.append( f"control log2 noise {control_noise:.3f} exceeds " f"{self.noise_threshold:.3f}" ) no_call_targets = tuple( target.target_id for index, target in enumerate(self.targets) if target.reportable and not self.callable_mask[index] ) if reasons: return DepthScreenResult( sample=sample, status="NO_CALL_QC", calls=(), no_call_targets=no_call_targets, median_control_depth=round(median_control_depth, 3), control_noise=round(control_noise, 4), reasons=tuple(reasons), ) candidates: list[tuple[float, _Window, str, float, float]] = [] for window in self.windows: weights = np.asarray(window.weights, dtype=float) projection = float(residual[list(window.indices)] @ weights) for kind, shift in (("DEL", -1.0), ("DUP", math.log2(1.5))): statistic = float( self._directional_statistic( np.asarray([projection]), window.precision_sum, shift )[0] ) threshold = self.family_threshold if statistic > threshold: mean_log2 = float(np.mean(residual[list(window.indices)])) candidates.append( (statistic, window, kind, mean_log2, threshold) ) # Keep the strongest non-overlapping explanation. This prevents one # event from producing a ladder of nested reflex requests. selected: list[tuple[float, _Window, str, float, float]] = [] occupied: set[int] = set() for row in sorted(candidates, key=lambda item: item[0], reverse=True): _, window, _, _, _ = row indices = set(window.indices) if occupied & indices: continue selected.append(row) occupied.update(indices) calls = tuple( DepthCall( gene=window.gene, kind=kind, first_exon=window.first_exon, last_exon=window.last_exon, target_ids=window.target_ids, n_targets=len(window.indices), mean_log2=round(mean_log2, 4), statistic=round(statistic, 4), threshold=round(threshold, 4), ) for statistic, window, kind, mean_log2, threshold in selected ) status = "REFLEX_POSITIVE" if calls else ( "NO_CANDIDATE_PARTIAL_COVERAGE" if no_call_targets else "NO_REFLEX_FLAG" ) return DepthScreenResult( sample=sample, status=status, calls=calls, no_call_targets=no_call_targets, median_control_depth=round(median_control_depth, 3), control_noise=round(control_noise, 4), reasons=(), ) def to_json(self) -> dict: return { "schema": SCHEMA, "targets": [asdict(target) for target in self.targets], "control_mask": self.control_mask.astype(int).tolist(), "centre": self.centre.tolist(), "components": self.components.tolist(), "residual_centre": self.residual_centre.tolist(), "residual_sigma": self.residual_sigma.tolist(), "callable_mask": self.callable_mask.astype(int).tolist(), "windows": [ { "gene": window.gene, "indices": list(window.indices), "target_ids": list(window.target_ids), "first_exon": window.first_exon, "last_exon": window.last_exon, "weights": list(window.weights), "precision_sum": window.precision_sum, } for window in self.windows ], "family_threshold": self.family_threshold, "noise_threshold": self.noise_threshold, "alpha": self.alpha, "n_fit": self.n_fit, "n_calibration": self.n_calibration, "n_components": self.n_components, "calibration_family_max": self.calibration_family_max.tolist(), "fit_parameters": { "requested_n_components": self.requested_n_components, "min_median_control_depth": self.min_median_control_depth, "min_target_median_depth": self.min_target_median_depth, "max_target_sigma": self.max_target_sigma, "max_window_targets": self.max_window_targets, }, "cohort_provenance": { "fit_sample_ids_sha256": self.fit_sample_ids_sha256, "calibration_sample_ids_sha256": self.calibration_sample_ids_sha256, "fit_count_rows_sha256": self.fit_count_rows_sha256, "calibration_count_rows_sha256": self.calibration_count_rows_sha256, }, "claims_boundary": CLAIMS_BOUNDARY, } @classmethod def from_json(cls, payload: dict) -> "AnchoredDepthModel": if not isinstance(payload, dict): raise ValueError("anchored depth model must be a JSON object") schema = payload.get("schema") if schema == LEGACY_SCHEMA: raise ValueError( "anchored depth model schema v1 is unsafe and cannot be " "migrated silently; refit the model to schema v2" ) if schema != SCHEMA: raise ValueError("unsupported anchored depth model schema") required = { "targets", "control_mask", "centre", "components", "residual_centre", "residual_sigma", "callable_mask", "windows", "family_threshold", "noise_threshold", "alpha", "n_fit", "n_calibration", "n_components", "calibration_family_max", "fit_parameters", "cohort_provenance", "claims_boundary", } observed_root_fields = set(payload) - {"schema"} missing = sorted(required - observed_root_fields) if missing: raise ValueError( "anchored depth model is missing field(s): " + ", ".join(missing) ) unexpected = sorted(observed_root_fields - required) if unexpected: raise ValueError( "anchored depth model has unexpected field(s): " + ", ".join(unexpected) ) raw_targets = payload["targets"] if not isinstance(raw_targets, list) or not raw_targets: raise ValueError("model targets must be a non-empty list") targets: list[DepthTarget] = [] for index, row in enumerate(raw_targets): if not isinstance(row, dict): raise ValueError(f"model target {index} is not an object") target_required = { "target_id", "gene", "exon", "order", "reportable" } if set(row) != target_required: raise ValueError(f"model target {index} is incomplete") target_id = row["target_id"] gene = row["gene"] exon = row["exon"] order = row["order"] reportable = row["reportable"] if ( not isinstance(target_id, str) or not target_id.strip() or not isinstance(gene, str) or not gene.strip() or not isinstance(exon, str) or not exon.strip() ): raise ValueError( "model target IDs, genes, and exons must be non-empty" ) if isinstance(order, bool) or not isinstance(order, int): raise ValueError("model target order must be an integer") if not isinstance(reportable, bool): raise ValueError( "model target reportable values must be booleans" ) targets.append( DepthTarget( target_id=target_id.strip(), gene=gene.strip(), exon=exon.strip(), order=order, reportable=reportable, ) ) target_rows = tuple(targets) expected = len(target_rows) if len({row.target_id for row in target_rows}) != expected: raise ValueError("model target IDs must be unique") if len({(row.gene, row.order) for row in target_rows}) != expected: raise ValueError( "model target order must be unique within each gene" ) def strict_integer(name: str, minimum: int) -> int: value = payload[name] if isinstance(value, bool) or not isinstance(value, int): raise ValueError(f"{name} must be an integer") if value < minimum: raise ValueError(f"{name} must be at least {minimum}") return int(value) def finite_number(value, name: str) -> float: if isinstance(value, bool) or not isinstance(value, (int, float)): raise ValueError(f"{name} must be numeric") parsed = float(value) if not np.isfinite(parsed): raise ValueError(f"{name} must be finite") return parsed def contains_non_numeric(value) -> bool: if isinstance(value, list): return any(contains_non_numeric(item) for item in value) return ( isinstance(value, bool) or not isinstance(value, (int, float)) ) def numeric_array( name: str, shape: tuple[int, ...], *, nonnegative: bool = False, positive: bool = False, ) -> np.ndarray: raw_value = payload[name] if contains_non_numeric(raw_value): raise ValueError( f"{name} must contain only JSON numeric values" ) try: value = np.asarray(raw_value, dtype=float) except (TypeError, ValueError) as exc: raise ValueError(f"{name} must be numeric") from exc if value.shape != shape: raise ValueError( f"{name} has shape {value.shape}; expected {shape}" ) if not np.all(np.isfinite(value)): raise ValueError(f"{name} must contain finite values") if nonnegative and np.any(value < 0): raise ValueError(f"{name} cannot contain negative values") if positive and np.any(value <= 0): raise ValueError(f"{name} must contain positive values") return value def boolean_vector(name: str) -> np.ndarray: raw = payload[name] if not isinstance(raw, list) or len(raw) != expected: raise ValueError( f"{name} must be a {expected}-element list" ) if any( isinstance(value, bool) is False and ( not isinstance(value, int) or value not in (0, 1) ) for value in raw ): raise ValueError(f"{name} must contain only booleans") return np.asarray(raw, dtype=bool) out = cls() out.targets = target_rows out.control_mask = boolean_vector("control_mask") expected_control = np.asarray( [not target.reportable for target in target_rows], dtype=bool, ) if not np.array_equal(out.control_mask, expected_control): raise ValueError( "model control mask differs from target reportability" ) if int(out.control_mask.sum()) < 20: raise ValueError("model contains fewer than 20 control targets") if int((~out.control_mask).sum()) < 1: raise ValueError("model contains no reportable targets") out.centre = numeric_array("centre", (expected,)) out.residual_centre = numeric_array( "residual_centre", (expected,) ) out.residual_sigma = numeric_array( "residual_sigma", (expected,), positive=True ) out.callable_mask = boolean_vector("callable_mask") if np.any(out.callable_mask & out.control_mask): raise ValueError( "model marks normalization controls as reportable-callable" ) if not np.any(out.callable_mask): raise ValueError("model has no callable reportable targets") out.alpha = finite_number(payload["alpha"], "alpha") if not 0 < out.alpha < 0.5: raise ValueError("model alpha must be finite and in (0, 0.5)") out.n_fit = strict_integer("n_fit", 50) out.n_calibration = strict_integer("n_calibration", 1) minimum_calibration = math.ceil(1.0 / out.alpha) - 1 if out.n_calibration < minimum_calibration: raise ValueError( "model calibration cohort cannot attain its family alpha" ) out.n_components = strict_integer("n_components", 0) maximum_components = min( out.n_fit - 2, int(out.control_mask.sum()) - 1, expected - 1, ) if out.n_components > maximum_components: raise ValueError( "model contains more components than its fit design permits" ) raw_components = payload["components"] if out.n_components == 0 and raw_components == []: out.components = np.zeros((0, expected), dtype=float) else: out.components = numeric_array( "components", (out.n_components, expected) ) if out.n_components and not np.allclose( out.components @ out.components.T, np.eye(out.n_components), rtol=1e-8, atol=1e-8, ): raise ValueError("model components must be orthonormal") fit_parameters = payload["fit_parameters"] if not isinstance(fit_parameters, dict): raise ValueError("fit_parameters must be an object") fit_required = { "requested_n_components", "min_median_control_depth", "min_target_median_depth", "max_target_sigma", "max_window_targets", } if set(fit_parameters) != fit_required: raise ValueError( "fit_parameters must contain exactly the schema-v2 fields" ) requested = fit_parameters["requested_n_components"] max_width = fit_parameters["max_window_targets"] if ( isinstance(requested, bool) or not isinstance(requested, int) or requested < 0 or isinstance(max_width, bool) or not isinstance(max_width, int) or max_width < 1 ): raise ValueError( "fit component and window parameters are invalid" ) out.requested_n_components = requested out.max_window_targets = max_width if out.n_components > out.requested_n_components: raise ValueError( "fitted components exceed requested components" ) scalar_parameters = ( "min_median_control_depth", "min_target_median_depth", "max_target_sigma", ) parsed_scalars: dict[str, float] = {} for name in scalar_parameters: value = finite_number( fit_parameters[name], f"fit parameter {name}" ) parsed_scalars[name] = value if ( parsed_scalars["min_median_control_depth"] <= 0 or parsed_scalars["min_target_median_depth"] < 0 or parsed_scalars["max_target_sigma"] <= 0 ): raise ValueError("fit depth/sigma parameters are invalid") out.min_median_control_depth = parsed_scalars[ "min_median_control_depth" ] out.min_target_median_depth = parsed_scalars[ "min_target_median_depth" ] out.max_target_sigma = parsed_scalars["max_target_sigma"] if np.any(out.residual_sigma < 0.03): raise ValueError( "model residual sigma is below the fitted variance floor" ) if np.any( out.callable_mask & (out.residual_sigma > out.max_target_sigma) ): raise ValueError( "model callable mask contradicts max_target_sigma" ) out.noise_threshold = finite_number( payload["noise_threshold"], "noise_threshold" ) out.family_threshold = finite_number( payload["family_threshold"], "family_threshold" ) if ( out.noise_threshold < 0 or out.family_threshold < 0 ): raise ValueError("model thresholds must be finite and non-negative") out.calibration_family_max = numeric_array( "calibration_family_max", (out.n_calibration,), nonnegative=True, ) expected_threshold = cls._conformal_threshold( out.calibration_family_max, out.alpha ) if not np.isclose( out.family_threshold, expected_threshold, rtol=1e-12, atol=1e-12, ): raise ValueError( "model family threshold differs from calibration maxima" ) expected_family = _window_index_family( out.targets, out.callable_mask, out.max_window_targets, ) raw_windows = payload["windows"] if not isinstance(raw_windows, list) or not raw_windows: raise ValueError("model windows must be a non-empty list") windows: list[_Window] = [] observed_family: list[tuple[str, tuple[int, ...]]] = [] for window_index, row in enumerate(raw_windows): if not isinstance(row, dict): raise ValueError( f"model window {window_index} is not an object" ) window_required = { "gene", "indices", "target_ids", "first_exon", "last_exon", "weights", "precision_sum", } if set(row) != window_required: raise ValueError( f"model window {window_index} is incomplete" ) raw_gene = row["gene"] gene = raw_gene.strip() if isinstance(raw_gene, str) else "" raw_indices = row["indices"] if ( not gene or not isinstance(raw_indices, list) or not raw_indices or any( isinstance(value, bool) or not isinstance(value, int) for value in raw_indices ) ): raise ValueError( f"model window {window_index} has invalid indices" ) indices = tuple(int(value) for value in raw_indices) if ( len(set(indices)) != len(indices) or any(value < 0 or value >= expected for value in indices) ): raise ValueError( f"model window {window_index} indices are invalid" ) selected = [out.targets[value] for value in indices] raw_target_ids = row["target_ids"] if ( not isinstance(raw_target_ids, list) or any( not isinstance(value, str) or not value for value in raw_target_ids ) ): raise ValueError( f"model window {window_index} target IDs are invalid" ) target_ids = tuple(raw_target_ids) if target_ids != tuple(value.target_id for value in selected): raise ValueError( f"model window {window_index} target IDs differ" ) first_exon = row["first_exon"] last_exon = row["last_exon"] if ( any(value.gene != gene for value in selected) or not np.all(out.callable_mask[list(indices)]) or not isinstance(first_exon, str) or not isinstance(last_exon, str) or first_exon != selected[0].exon or last_exon != selected[-1].exon ): raise ValueError( f"model window {window_index} target contract differs" ) raw_weights = row["weights"] if ( not isinstance(raw_weights, list) or contains_non_numeric(raw_weights) ): raise ValueError( f"model window {window_index} weights are invalid" ) try: weights = tuple(float(value) for value in raw_weights) except (TypeError, ValueError) as exc: raise ValueError( f"model window {window_index} weights are invalid" ) from exc precision_sum = finite_number( row["precision_sum"], f"model window {window_index} precision_sum", ) if ( len(weights) != len(indices) or not np.all(np.isfinite(weights)) or not np.isfinite(precision_sum) or precision_sum <= 0 or not np.isclose( sum(weights), precision_sum, rtol=1e-6, atol=1e-8, ) ): raise ValueError( f"model window {window_index} precision is invalid" ) windows.append( _Window( gene=gene, indices=indices, target_ids=target_ids, first_exon=selected[0].exon, last_exon=selected[-1].exon, weights=weights, precision_sum=precision_sum, ) ) observed_family.append((gene, indices)) if tuple(observed_family) != expected_family: raise ValueError( "model windows differ from the frozen callable family" ) out.windows = tuple(windows) provenance = payload["cohort_provenance"] if not isinstance(provenance, dict): raise ValueError("cohort_provenance must be an object") provenance_required = { "fit_sample_ids_sha256", "calibration_sample_ids_sha256", "fit_count_rows_sha256", "calibration_count_rows_sha256", } if set(provenance) != provenance_required: raise ValueError( "cohort_provenance must contain exactly the schema-v2 fields" ) out.fit_sample_ids_sha256 = _valid_sha256( provenance["fit_sample_ids_sha256"], optional=True, name="fit_sample_ids_sha256", ) out.calibration_sample_ids_sha256 = _valid_sha256( provenance["calibration_sample_ids_sha256"], optional=True, name="calibration_sample_ids_sha256", ) out.fit_count_rows_sha256 = _valid_sha256( provenance["fit_count_rows_sha256"], optional=False, name="fit_count_rows_sha256", ) out.calibration_count_rows_sha256 = _valid_sha256( provenance["calibration_count_rows_sha256"], optional=False, name="calibration_count_rows_sha256", ) if ( out.fit_sample_ids_sha256 is not None and out.fit_sample_ids_sha256 == out.calibration_sample_ids_sha256 ): raise ValueError( "fit and calibration sample-ID hashes must be disjoint" ) if ( out.fit_count_rows_sha256 == out.calibration_count_rows_sha256 ): raise ValueError( "fit and calibration count-row hashes must be disjoint" ) if payload["claims_boundary"] != CLAIMS_BOUNDARY: raise ValueError("model claims boundary is missing or altered") return out def save(self, path: str | Path) -> None: destination = Path(path) atomic_write_json(str(destination), self.to_json(), indent=2) @classmethod def load(cls, path: str | Path) -> "AnchoredDepthModel": return cls.from_json(json.loads(Path(path).read_text(encoding="utf-8")))