"""Cross-gene anchored fallback for the supplied BRCA-only count matrix. This is intentionally a constrained screen: BRCA2 targets normalize BRCA1 and BRCA1 targets normalize BRCA2. It preserves gene-wide dosage shifts that are erased by within-gene scaling, but it cannot identify simultaneous changes in both genes and is not a substitute for a broad assay normalization domain. """ from __future__ import annotations import argparse import csv import hashlib import json import math import platform from dataclasses import asdict from datetime import datetime, timezone from pathlib import Path import numpy as np from lgasieve.clinical_depth import AnchoredDepthModel, DepthTarget from lgasieve.provenance import atomic_write_json BASE_SEED = 20260731 COMPONENT_GRID = (0, 2, 4, 6, 8) ALPHA_GRID = (0.05, 0.10, 0.20) def _sha256(path: Path) -> str: digest = hashlib.sha256() with path.open("rb") as handle: for block in iter(lambda: handle.read(1024 * 1024), b""): digest.update(block) return digest.hexdigest() def load_matrix(path: Path) -> tuple[list[dict], tuple[str, ...], np.ndarray]: with path.open(newline="", encoding="utf-8-sig") as handle: reader = csv.DictReader(handle, delimiter="\t") metadata = ( "target_id", "gene", "chrom", "start", "end", "exon", "usable", ) if not set(metadata).issubset(reader.fieldnames or ()): raise ValueError("fragment matrix is missing target metadata") samples = tuple( column for column in reader.fieldnames or () if column not in metadata ) rows = list(reader) counts = np.asarray( [[float(row[sample]) for row in rows] for sample in samples], dtype=float, ) return rows, samples, counts def build_targets( rows: list[dict], reportable_gene: str ) -> tuple[tuple[DepthTarget, ...], np.ndarray, tuple[str, ...]]: usable_indices = np.asarray( [index for index, row in enumerate(rows) if int(row["usable"]) == 1], dtype=int, ) targets = tuple( DepthTarget( target_id=rows[index]["target_id"], gene=rows[index]["gene"], exon=rows[index]["exon"], order=int(rows[index]["exon"]), reportable=rows[index]["gene"] == reportable_gene, ) for index in usable_indices ) excluded = tuple( row["target_id"] for row in rows if row["gene"] == reportable_gene and int(row["usable"]) != 1 ) return targets, usable_indices, excluded def _event_catalog(rows: list[dict]) -> list[dict]: catalog: list[dict] = [] for gene in ("BRCA1", "BRCA2"): indices = sorted( ( index for index, row in enumerate(rows) if row["gene"] == gene ), key=lambda index: int(rows[index]["exon"]), ) for width in (1, 2, 3, 5): for start in range(len(indices) - width + 1): event_indices = indices[start : start + width] for kind, factor in (("DEL", 0.5), ("DUP", 1.5)): catalog.append( { "gene": gene, "kind": kind, "factor": factor, "class": ( "single" if width == 1 else f"width_{width}" ), "matrix_indices": event_indices, "target_ids": [ rows[index]["target_id"] for index in event_indices ], } ) for kind, factor in (("DEL", 0.5), ("DUP", 1.5)): catalog.append( { "gene": gene, "kind": kind, "factor": factor, "class": "whole_gene", "matrix_indices": indices, "target_ids": [ rows[index]["target_id"] for index in indices ], } ) return catalog def _matches(call, event: dict, callable_ids: set[str]) -> bool: if call.gene != event["gene"] or call.kind != event["kind"]: return False expected = set(event["target_ids"]) & callable_ids called = set(call.target_ids) overlap = expected & called return bool( expected and overlap and len(overlap) / len(expected) >= 0.5 and len(overlap) / len(called) >= 0.5 ) def evaluate_variant( models: dict[str, AnchoredDepthModel], matrices: dict[str, np.ndarray], rows: list[dict], development_indices: np.ndarray, *, n_events: int = 2400, ) -> dict: rng = np.random.default_rng(BASE_SEED + 700) catalog = _event_catalog(rows) choices = rng.integers(0, len(catalog), n_events) hosts = rng.choice(development_indices, n_events, replace=True) callable_ids = { gene: { target.target_id for target, callable_target in zip( model.targets, model.callable_mask ) if target.reportable and callable_target } for gene, model in models.items() } presumed_negative_referrals = 0 negative_no_calls = 0 for sample_index in development_indices: calls = [] qc_failed = False for gene, model in models.items(): result = model.screen( matrices[gene][sample_index], sample=f"normal_{sample_index}" ) calls.extend(result.calls) qc_failed = qc_failed or result.status == "NO_CALL_QC" negative_no_calls += int(qc_failed) presumed_negative_referrals += int(bool(calls) and not qc_failed) positive_rows = [] usable_lookup = { gene: { target.target_id: index for index, target in enumerate(model.targets) } for gene, model in models.items() } full_target_lookup = { row["target_id"]: index for index, row in enumerate(rows) } for event_index, (choice, host) in enumerate(zip(choices, hosts)): event = catalog[int(choice)] gene = event["gene"] model = models[gene] vector = matrices[gene][int(host)].copy() for target_id in event["target_ids"]: if target_id in usable_lookup[gene]: vector[usable_lookup[gene][target_id]] *= event["factor"] # Deterministic multiplication is deliberately retained as a # conditional spike-in on the actual library background; the main # full-panel experiment separately includes count sampling. result = model.screen(vector, sample=f"spike_{event_index}") callable_event = bool( set(event["target_ids"]) & callable_ids[gene] ) completed = result.status != "NO_CALL_QC" and callable_event detected = completed and any( _matches(call, event, callable_ids[gene]) for call in result.calls ) positive_rows.append( { "completed": completed, "detected": detected, "kind": event["kind"], "class": event["class"], } ) completed = [row for row in positive_rows if row["completed"]] detected = sum(row["detected"] for row in completed) completed_negatives = len(development_indices) - negative_no_calls metrics = { "positive_spike_ins": len(positive_rows), "positive_completed": len(completed), "positive_no_calls": len(positive_rows) - len(completed), "sensitivity_completed": ( detected / len(completed) if completed else None ), "negative_samples": len(development_indices), "negative_no_calls": negative_no_calls, "presumed_negative_referral_samples": presumed_negative_referrals, "presumed_negative_no_referral_rate_completed": ( 1.0 - presumed_negative_referrals / completed_negatives if completed_negatives else None ), "strata": {}, } for key in sorted( {f"kind:{row['kind']}" for row in positive_rows} | {f"class:{row['class']}" for row in positive_rows} ): field, value = key.split(":", 1) subset = [ row for row in positive_rows if row[field] == value and row["completed"] ] metrics["strata"][key] = { "completed": len(subset), "detected": sum(row["detected"] for row in subset), "sensitivity_completed": ( sum(row["detected"] for row in subset) / len(subset) if subset else None ), } return metrics def _utility( metrics: dict, callable_targets: int, total_usable_reportable: int ) -> tuple: sensitivity = metrics["sensitivity_completed"] or 0.0 no_referral_rate = ( metrics["presumed_negative_no_referral_rate_completed"] or 0.0 ) single = ( metrics["strata"].get("class:single", {}) .get("sensitivity_completed") or 0.0 ) callable_fraction = callable_targets / total_usable_reportable gates = ( int(sensitivity >= 0.95) + int(single >= 0.90) + int(no_referral_rate >= 0.90) + int(callable_fraction >= 0.85) ) return ( gates, sensitivity, single, no_referral_rate, callable_fraction, ) def run( matrix_path: Path, partition_path: Path, evidence_dir: Path, output_path: Path, model_prefix: Path, ) -> dict: rows, samples, counts = load_matrix(matrix_path) partition_payload = json.loads( partition_path.read_text(encoding="utf-8") ) partition = partition_payload["partition"] sample_index = {sample: index for index, sample in enumerate(samples)} fit_indices = np.asarray( [sample_index[sample] for sample in partition["reference"]] ) calibration_indices = np.asarray( [sample_index[sample] for sample in partition["calibration"]] ) development_indices = np.asarray( [sample_index[sample] for sample in partition["development"]] ) fit_sample_ids = tuple(partition["reference"]) calibration_sample_ids = tuple(partition["calibration"]) target_sets = {} matrices = {} usable_by_gene = {} excluded = {} for gene in ("BRCA1", "BRCA2"): targets, usable_indices, excluded_targets = build_targets(rows, gene) target_sets[gene] = targets matrices[gene] = counts[:, usable_indices] usable_by_gene[gene] = usable_indices excluded[gene] = excluded_targets variants = [] candidate_models = {} total_usable = sum(int(row["usable"]) for row in rows) for components in COMPONENT_GRID: for alpha in ALPHA_GRID: models = { gene: AnchoredDepthModel.fit( matrices[gene][fit_indices], matrices[gene][calibration_indices], target_sets[gene], alpha=alpha, n_components=components, min_median_control_depth=30.0, min_target_median_depth=20.0, max_target_sigma=0.60, fit_sample_ids=fit_sample_ids, calibration_sample_ids=calibration_sample_ids, ) for gene in ("BRCA1", "BRCA2") } metrics = evaluate_variant( models, matrices, rows, development_indices ) callable_targets = sum( int(model.callable_mask.sum()) for model in models.values() ) variants.append( { "components": components, "alpha": alpha, "callable_targets": callable_targets, "metrics": metrics, "selection_utility": list( _utility(metrics, callable_targets, total_usable) ), } ) candidate_models[(components, alpha)] = models selected = max( variants, key=lambda row: ( tuple(row["selection_utility"]), -row["components"], -row["alpha"], ), ) models = candidate_models[ (selected["components"], selected["alpha"]) ] model_paths = {} for gene, model in models.items(): path = model_prefix.with_name( f"{model_prefix.name}_{gene.lower()}.json" ) model.save(path) model_paths[gene] = { "path": str(path.resolve()), "sha256": _sha256(path), } def screen_operational( names: tuple[str, ...], count_matrices: dict[str, np.ndarray], ) -> tuple[list[dict], dict[str, int]]: operational_rows = [] for index, sample in enumerate(names): calls = [] qc_reasons = [] for gene, model in models.items(): result = model.screen( count_matrices[gene][index], sample=sample ) calls.extend(asdict(call) for call in result.calls) if result.status == "NO_CALL_QC": qc_reasons.extend( f"{gene}: {reason}" for reason in result.reasons ) status = ( "NO_CALL_QC" if qc_reasons else "REFLEX_POSITIVE" if calls else "NO_CANDIDATE_PARTIAL_COVERAGE" ) operational_rows.append( { "sample": sample, "status": status, "calls": calls, "excluded_targets": sorted( excluded["BRCA1"] + excluded["BRCA2"] ), "reasons": qc_reasons, } ) status_counts_local = {} for result in operational_rows: status = result["status"] status_counts_local[status] = ( status_counts_local.get(status, 0) + 1 ) return operational_rows, status_counts_local operational, status_counts = screen_operational(samples, matrices) expected_target_ids = [row["target_id"] for row in rows] evidence_names = [] evidence_full_counts = [] evidence_median_depth = {} for evidence_path in sorted(evidence_dir.glob("*.evidence.json")): evidence = json.loads( evidence_path.read_text(encoding="utf-8") ) observed = { target_id: count for target_id, count in zip( evidence["target_ids"], evidence["fragments"] ) } if set(observed) != set(expected_target_ids): raise ValueError( f"evidence target IDs differ for {evidence_path.name}" ) evidence_names.append(str(evidence["sample"])) evidence_median_depth[str(evidence["sample"])] = float( evidence["qc"]["median_target_depth"] ) evidence_full_counts.append( [ float(observed[target_id]) for target_id in expected_target_ids ] ) evidence_matrix = np.asarray(evidence_full_counts, dtype=float) evidence_matrices = { gene: evidence_matrix[:, usable_by_gene[gene]] for gene in ("BRCA1", "BRCA2") } evidence_operational, evidence_status_counts = screen_operational( tuple(evidence_names), evidence_matrices ) for result in evidence_operational: median_depth = evidence_median_depth[result["sample"]] result["source_median_target_depth"] = median_depth if median_depth < 30.0: result["status"] = "NO_CALL_QC" result["calls"] = [] result["reasons"].append( f"source median target depth {median_depth:.1f} is below 30.0" ) evidence_status_counts = {} for result in evidence_operational: status = result["status"] evidence_status_counts[status] = ( evidence_status_counts.get(status, 0) + 1 ) payload = { "schema": "lgasieve.cross-gene-reflex-reanalysis.v2", "created_utc": datetime.now(timezone.utc) .isoformat() .replace("+00:00", "Z"), "claims_boundary": [ "BRCA2 anchors BRCA1 and BRCA1 anchors BRCA2; this is a " "two-gene fallback, not the broad-panel clinical model.", "Development libraries are presumed, not verified, negatives.", "Conditional spike-ins are internal software tests, not clinical " "sensitivity.", "Every candidate is REFLEX_POSITIVE and requires MLPA, aCGH, " "ddPCR, or another independently designed dosage assay.", "A simultaneous BRCA1/BRCA2 dosage shift remains non-identifiable.", ], "inputs": { "matrix": { "path": str(matrix_path.resolve()), "sha256": _sha256(matrix_path), }, "partition": { "path": str(partition_path.resolve()), "sha256": _sha256(partition_path), }, "n_samples": len(samples), "n_targets": len(rows), "n_usable_targets": total_usable, "excluded_targets": excluded, "evidence_directory": str(evidence_dir.resolve()), "n_evidence_libraries": len(evidence_names), }, "development_variants": variants, "selected": selected, "models": model_paths, "operational_status_counts": status_counts, "operational_results": operational, "evidence_operational_status_counts": evidence_status_counts, "evidence_operational_results": evidence_operational, "environment": { "python": platform.python_version(), "numpy": np.__version__, }, } atomic_write_json(str(output_path), payload) return payload def main() -> int: parser = argparse.ArgumentParser() parser.add_argument( "--matrix", type=Path, default=Path( "work/exomedepth_development_v2/fragment_counts.tsv" ), ) parser.add_argument( "--partition", type=Path, default=Path("work/exomedepth_development_v2/partition.json"), ) parser.add_argument( "--evidence-dir", type=Path, default=Path("validation_corrected/evidence"), ) parser.add_argument( "--output", type=Path, default=Path("analysis/cross_gene_reflex_reanalysis.json"), ) parser.add_argument( "--model-prefix", type=Path, default=Path("analysis/cross_gene_reflex_model"), ) args = parser.parse_args() payload = run( args.matrix, args.partition, args.evidence_dir, args.output, args.model_prefix, ) print( json.dumps( { "selected": payload["selected"], "operational_status_counts": payload[ "operational_status_counts" ], "evidence_operational_status_counts": payload[ "evidence_operational_status_counts" ], }, indent=2, ) ) return 0 if __name__ == "__main__": raise SystemExit(main())