#!/usr/bin/env python3 """Build a bounded VCC-style metric diagnostic from aggregate perturbation deltas. This adapter intentionally does not claim an official Virtual Cell Challenge execution. It implements the pseudobulk-delta formulas that can be supported by the frozen NMD-VCell prediction bundle and keeps cell-level/DE metrics locked. """ from __future__ import annotations import argparse import csv import hashlib import json import math from pathlib import Path from typing import Any, Callable import numpy as np BUILD_ID = "VCC-METRICS-20260729-02" BOOTSTRAP_SEED = 20260729 BOOTSTRAP_REPLICATES = 10_000 SCALE_GRID = (0.0, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0) def sha256_file(path: Path) -> str: digest = hashlib.sha256() with path.open("rb") as handle: for chunk in iter(lambda: handle.read(1024 * 1024), b""): digest.update(chunk) return digest.hexdigest() def finite_float(value: float | np.floating[Any] | None) -> float | None: if value is None: return None parsed = float(value) return parsed if math.isfinite(parsed) else None def rankdata_average(values: np.ndarray) -> np.ndarray: """Return one-based average ranks with deterministic stable tie handling.""" values = np.asarray(values, dtype=np.float64) order = np.argsort(values, kind="mergesort") sorted_values = values[order] ranks = np.empty(values.size, dtype=np.float64) start = 0 while start < values.size: end = start + 1 while end < values.size and sorted_values[end] == sorted_values[start]: end += 1 ranks[order[start:end]] = (start + 1 + end) / 2.0 start = end return ranks def pearson(x: np.ndarray, y: np.ndarray) -> float | None: x = np.asarray(x, dtype=np.float64) y = np.asarray(y, dtype=np.float64) if x.size < 2 or y.size != x.size: return None x_centered = x - x.mean() y_centered = y - y.mean() denominator = float( np.sqrt(np.dot(x_centered, x_centered) * np.dot(y_centered, y_centered)) ) if denominator == 0.0: return None return finite_float(np.dot(x_centered, y_centered) / denominator) def spearman(x: np.ndarray, y: np.ndarray) -> float | None: return pearson(rankdata_average(x), rankdata_average(y)) def cosine(x: np.ndarray, y: np.ndarray) -> float | None: x = np.asarray(x, dtype=np.float64) y = np.asarray(y, dtype=np.float64) denominator = float(np.linalg.norm(x) * np.linalg.norm(y)) if denominator == 0.0: return None return finite_float(np.dot(x, y) / denominator) def bootstrap_mean( values: np.ndarray, *, seed: int = BOOTSTRAP_SEED, replicates: int = BOOTSTRAP_REPLICATES, ) -> tuple[float | None, float | None]: values = np.asarray(values, dtype=np.float64) values = values[np.isfinite(values)] if values.size == 0: return None, None rng = np.random.default_rng(seed) indices = rng.integers(0, values.size, size=(replicates, values.size)) boot = values[indices].mean(axis=1) low, high = np.quantile(boot, [0.025, 0.975]) return float(low), float(high) def summarize_target_values( values: list[float | None], *, direction: str, seed_offset: int, ) -> dict[str, Any]: finite = np.asarray( [value for value in values if value is not None and math.isfinite(value)], dtype=np.float64, ) if finite.size == 0: return { "estimate": None, "median": None, "ci95": [None, None], "n_targets": 0, "direction": direction, } low, high = bootstrap_mean(finite, seed=BOOTSTRAP_SEED + seed_offset) return { "estimate": float(finite.mean()), "median": float(np.median(finite)), "ci95": [low, high], "n_targets": int(finite.size), "direction": direction, } def pds_l1( prediction: np.ndarray, actual: np.ndarray, perturbations: np.ndarray, response_genes: np.ndarray, ) -> tuple[float, list[float], list[int]]: """Compute cell-eval-aligned L1 discrimination on aggregate deltas. The perturbed target gene is excluded when it is present in the response coordinate set. Rank 0 is the best match and score is 1 - rank / N. """ n_targets = prediction.shape[0] target_scores: list[float] = [] target_ranks: list[int] = [] for target_index, perturbation in enumerate(perturbations): include = response_genes != perturbation distances = np.abs( actual[:, include] - prediction[target_index, include] ).sum(axis=1) sorted_indices = np.argsort(distances, kind="mergesort") rank = int(np.flatnonzero(sorted_indices == target_index)[0]) target_ranks.append(rank + 1) target_scores.append(1.0 - rank / n_targets) return float(np.mean(target_scores)), target_scores, target_ranks def norm_match_predictions( prediction: np.ndarray, actual: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: """Scale each prediction to its paired observed L2 norm. This deliberately uses hidden truth and is therefore a sensitivity diagnostic, never a competition-eligible score. """ prediction_norms = np.linalg.norm(prediction, axis=1) actual_norms = np.linalg.norm(actual, axis=1) factors = np.zeros_like(prediction_norms, dtype=np.float64) nonzero = prediction_norms > 0 factors[nonzero] = actual_norms[nonzero] / prediction_norms[nonzero] return prediction * factors[:, None], factors def target_metrics( prediction: np.ndarray, actual: np.ndarray, train_mean: np.ndarray, perturbations: np.ndarray, response_genes: np.ndarray, ) -> tuple[list[dict[str, Any]], dict[str, Any]]: pds_score, pds_scores, pds_ranks = pds_l1( prediction, actual, perturbations, response_genes ) records: list[dict[str, Any]] = [] for index, gene in enumerate(perturbations): pred = prediction[index] truth = actual[index] records.append( { "gene": str(gene), "rmse_delta": float(np.sqrt(np.mean((pred - truth) ** 2))), "mae_delta": float(np.mean(np.abs(pred - truth))), "delta_pearson": pearson(pred, truth), "delta_spearman": spearman(pred, truth), "raw_cosine": cosine(pred, truth), "residual_cosine_vs_train_mean": cosine( pred - train_mean, truth - train_mean ), "prediction_l2_norm": float(np.linalg.norm(pred)), "actual_l2_norm": float(np.linalg.norm(truth)), "pds_l1_raw": float(pds_scores[index]), "pds_l1_rank": int(pds_ranks[index]), "target_gene_excluded_from_pds": bool( np.any(response_genes == gene) ), } ) summaries: dict[str, Any] = {} metric_specs = ( ("rmse_delta", "lower_is_better"), ("mae_delta", "lower_is_better"), ("delta_pearson", "higher_is_better"), ("delta_spearman", "higher_is_better"), ("raw_cosine", "higher_is_better"), ("residual_cosine_vs_train_mean", "higher_is_better"), ("pds_l1_raw", "higher_is_better"), ) for offset, (metric_id, direction) in enumerate(metric_specs): summaries[metric_id] = summarize_target_values( [record[metric_id] for record in records], direction=direction, seed_offset=offset, ) summaries["pds_l1_raw"]["estimate"] = pds_score summaries["effect_size_spearman"] = { "estimate": spearman( np.linalg.norm(prediction, axis=1), np.linalg.norm(actual, axis=1) ), "median": None, "ci95": [None, None], "n_targets": int(prediction.shape[0]), "direction": "higher_is_better", "definition": "Spearman correlation across perturbations between predicted and observed delta L2 norms.", } return records, summaries def load_target_manifest(path: Path | None) -> dict[str, dict[str, str]]: if path is None: return {} with path.open(newline="", encoding="utf-8") as handle: return {row["gene"]: row for row in csv.DictReader(handle, delimiter="\t")} def load_summary(path: Path | None) -> dict[str, Any]: if path is None: return {} return json.loads(path.read_text(encoding="utf-8")) def validate_bundle(bundle: Any) -> dict[str, np.ndarray]: required = ( "genes", "response_genes", "prediction", "actual", "train_mean", "coordinate_sd", ) missing = [key for key in required if key not in bundle.files] if missing: raise ValueError(f"Prediction bundle missing keys: {', '.join(missing)}") arrays = {key: np.asarray(bundle[key]) for key in required} genes = arrays["genes"] response_genes = arrays["response_genes"] prediction = arrays["prediction"].astype(np.float64) actual = arrays["actual"].astype(np.float64) train_mean = arrays["train_mean"].astype(np.float64) if prediction.shape != actual.shape: raise ValueError("prediction and actual matrices must have identical shapes") if prediction.shape != (genes.size, response_genes.size): raise ValueError("matrix dimensions do not match gene axes") if train_mean.shape != (response_genes.size,): raise ValueError("train_mean must have one value per response coordinate") for label, matrix in ( ("prediction", prediction), ("actual", actual), ("train_mean", train_mean), ): if not np.isfinite(matrix).all(): raise ValueError(f"{label} contains non-finite values") arrays["prediction"] = prediction arrays["actual"] = actual arrays["train_mean"] = train_mean return arrays def diagnostic_object( *, bundle_path: Path, summary_path: Path | None, manifest_path: Path | None, ) -> tuple[dict[str, Any], list[dict[str, Any]]]: with np.load(bundle_path, allow_pickle=False) as bundle: arrays = validate_bundle(bundle) perturbations = arrays["genes"].astype(str) response_genes = arrays["response_genes"].astype(str) prediction = arrays["prediction"] actual = arrays["actual"] train_mean = arrays["train_mean"] source_summary = load_summary(summary_path) manifest = load_target_manifest(manifest_path) method_inputs = { "safe_ridge_transfer": prediction, "training_response_mean": np.repeat( train_mean[None, :], prediction.shape[0], axis=0 ), "zero_change": np.zeros_like(prediction), } target_rows_by_method: dict[str, list[dict[str, Any]]] = {} scoreboard: list[dict[str, Any]] = [] labels = { "safe_ridge_transfer": "Safe ridge transfer", "training_response_mean": "Training-response mean", "zero_change": "Zero change", } for method_index, (method_id, matrix) in enumerate(method_inputs.items()): rows, metrics = target_metrics( matrix, actual, train_mean, perturbations, response_genes ) target_rows_by_method[method_id] = rows scoreboard.append( { "method_id": method_id, "label": labels[method_id], "role": "evaluated_model" if method_id == "safe_ridge_transfer" else "baseline", "metrics": metrics, "interpretation": "Post-pilot diagnostic model; source conclusion remains NO_DIRECTIONAL_TRANSFER_SUPPORT." if method_id == "safe_ridge_transfer" else "Comparator evaluated on the identical 55-target aggregate-delta matrix.", } ) scale_sensitivity: list[dict[str, Any]] = [] for scale in SCALE_GRID: score, _, _ = pds_l1( prediction * scale, actual, perturbations, response_genes ) scale_sensitivity.append( {"prediction_scale": scale, "pds_l1_raw": score} ) norm_matched_prediction, norm_factors = norm_match_predictions( prediction, actual ) norm_score, norm_scores, norm_ranks = pds_l1( norm_matched_prediction, actual, perturbations, response_genes ) norm_low, norm_high = bootstrap_mean( np.asarray(norm_scores), seed=BOOTSTRAP_SEED + 101 ) target_rows: list[dict[str, Any]] = [] for target_index, gene in enumerate(perturbations): metadata = manifest.get(str(gene), {}) combined = { "gene": str(gene), "external_order": int(target_index), "n_cells": int(metadata["n_cells"]) if metadata.get("n_cells") else None, "present_in_hepg2_2160": metadata.get( "present_in_hepg2_2160" ), "model": target_rows_by_method["safe_ridge_transfer"][target_index], "training_response_mean": target_rows_by_method[ "training_response_mean" ][target_index], "zero_change": target_rows_by_method["zero_change"][target_index], "norm_matched_pds_l1": float(norm_scores[target_index]), "norm_matched_pds_rank": int(norm_ranks[target_index]), "truth_derived_norm_factor": float(norm_factors[target_index]), } target_rows.append(combined) source_files: list[dict[str, Any]] = [ { "role": "aggregate_prediction_truth_bundle", "filename": bundle_path.name, "sha256": sha256_file(bundle_path), } ] for role, path in ( ("source_result_summary", summary_path), ("external_target_manifest", manifest_path), ): if path is not None: source_files.append( { "role": role, "filename": path.name, "sha256": sha256_file(path), } ) result = { "diagnostic_schema": "nmd-vcell-vcc-metric-diagnostic/1.0", "feature_build": BUILD_ID, "assessed_at": "2026-07-29", "execution_state": "EXECUTED_AGGREGATE_DELTA_DIAGNOSTIC_NOT_OFFICIAL_VCC", "scientific_context": { "source_stage": source_summary.get("stage", "4.16i"), "source_mode": source_summary.get("mode", "real"), "source_result_label": source_summary.get( "result_label", "NO_DIRECTIONAL_TRANSFER_SUPPORT" ), "post_pilot_diagnostic_only": source_summary.get( "post_pilot_diagnostic_only", True ), "training_context": "HepG2 aggregate perturbation deltas", "external_truth_context": "Frangieh external perturbation response panel", "n_external_targets": int(prediction.shape[0]), "n_response_coordinates": int(prediction.shape[1]), "prediction_unit": "target_level_aggregate_delta", }, "formula_alignment": { "mae_delta": "Mean absolute difference between predicted and observed aggregate perturbation deltas, averaged within target then across targets.", "delta_pearson": "Pearson correlation between predicted and observed perturbation deltas within each target.", "pds_l1_raw": "cell-eval-aligned L1 rank retrieval on aggregate deltas; the perturbed target gene is excluded when present.", "pds_norm_matched": "Prediction vectors are scaled to paired observed L2 norms before L1 retrieval; uses hidden truth and is sensitivity-only.", "upstream": "https://github.com/ArcInstitute/cell-eval/blob/main/src/cell_eval/metrics/_anndata.py", "adapter_difference": "The upstream package consumes predicted and observed AnnData and can compute cell-level and DE metrics. This adapter consumes frozen target-level aggregate deltas only.", }, "uncertainty": { "unit": "external_perturbation_target", "method": "nonparametric target bootstrap", "seed": BOOTSTRAP_SEED, "replicates": BOOTSTRAP_REPLICATES, "note": "Intervals quantify target-sampling variation within this diagnostic panel, not donor, dataset or disease-stage uncertainty.", }, "scoreboard": scoreboard, "pds_scale_sensitivity": { "model_id": "safe_ridge_transfer", "scale_grid": scale_sensitivity, "norm_matched_truth_informed": { "pds_l1": norm_score, "ci95": [norm_low, norm_high], "n_targets": int(prediction.shape[0]), "eligibility": "DIAGNOSTIC_ONLY_TRUTH_DERIVED_SCALING", "median_scaling_factor": float(np.median(norm_factors)), }, "interpretation": "Raw and truth-norm-matched scores are displayed together to reveal scale sensitivity; the norm-matched value is not a fair hidden-test score.", }, "locked_metrics": [ { "metric_id": "des_deg_recovery", "state": "LOCKED_REQUIRES_CELL_LEVEL_OR_DE_TRUTH", }, { "metric_id": "deg_auprc", "state": "LOCKED_REQUIRES_CELL_LEVEL_OR_DE_TRUTH", }, { "metric_id": "single_cell_distribution_distance", "state": "LOCKED_AGGREGATE_PREDICTIONS_ONLY", }, { "metric_id": "cell_state_proportion_recovery", "state": "LOCKED_AGGREGATE_PREDICTIONS_ONLY", }, ], "source_files": source_files, "conclusion": "The multi-metric adapter is operational on the existing 55-target external aggregate panel. It does not overturn the frozen NO_DIRECTIONAL_TRANSFER_SUPPORT conclusion and is not evidence of VCC 2026 task compatibility.", "claim_boundary": "This is a formula-aligned aggregate-delta diagnostic, not an official cell-eval run, not a VCC 2025 rescore, and not a VCC 2026 task execution. No muscle, DMD, treatment or clinical validity follows from these scores.", } return result, target_rows def flatten_target_rows(rows: list[dict[str, Any]]) -> list[dict[str, Any]]: flat_rows: list[dict[str, Any]] = [] for row in rows: flat: dict[str, Any] = { "gene": row["gene"], "external_order": row["external_order"], "n_cells": row["n_cells"], "present_in_hepg2_2160": row["present_in_hepg2_2160"], } for method_id in ("model", "training_response_mean", "zero_change"): for key, value in row[method_id].items(): if key != "gene": flat[f"{method_id}_{key}"] = value flat["norm_matched_pds_l1"] = row["norm_matched_pds_l1"] flat["norm_matched_pds_rank"] = row["norm_matched_pds_rank"] flat["truth_derived_norm_factor"] = row["truth_derived_norm_factor"] flat_rows.append(flat) return flat_rows def write_tsv(path: Path, rows: list[dict[str, Any]]) -> None: path.parent.mkdir(parents=True, exist_ok=True) if not rows: path.write_text("", encoding="utf-8") return with path.open("w", newline="", encoding="utf-8") as handle: writer = csv.DictWriter( handle, fieldnames=list(rows[0].keys()), delimiter="\t", lineterminator="\n", ) writer.writeheader() writer.writerows(rows) def format_metric(metric: dict[str, Any]) -> str: estimate = metric.get("estimate") ci95 = metric.get("ci95", [None, None]) if estimate is None: return "not defined" text = f"{estimate:.6f}" if ci95[0] is not None and ci95[1] is not None: text += f" ({ci95[0]:.6f} to {ci95[1]:.6f})" return text def markdown_report(result: dict[str, Any]) -> str: lines = [ "# VCC-style aggregate-delta metric execution", "", f"Assessed: {result['assessed_at']}", "", "## Scope", "", "This run applies formula-aligned multi-metric diagnostics to the frozen 55-target external Frangieh aggregate-delta panel. It is not an official `cell-eval` AnnData execution and is not a VCC 2026 task result.", "", "## Scoreboard", "", "| method | RMSE delta | MAE delta | PDS L1 | delta Pearson | delta Spearman | effect-size Spearman |", "| --- | ---: | ---: | ---: | ---: | ---: | ---: |", ] for method in result["scoreboard"]: metrics = method["metrics"] lines.append( "| {label} | {rmse} | {mae} | {pds} | {pearson} | {spearman_value} | {effect} |".format( label=method["label"], rmse=format_metric(metrics["rmse_delta"]), mae=format_metric(metrics["mae_delta"]), pds=format_metric(metrics["pds_l1_raw"]), pearson=format_metric(metrics["delta_pearson"]), spearman_value=format_metric(metrics["delta_spearman"]), effect=format_metric(metrics["effect_size_spearman"]), ) ) lines.extend( [ "", "Intervals are 95% target-bootstrap intervals with 10,000 replicates. Undefined correlations for constant baselines remain null rather than being converted to zero.", "", "## PDS scale audit", "", "| prediction scale | raw L1 PDS |", "| ---: | ---: |", ] ) for row in result["pds_scale_sensitivity"]["scale_grid"]: lines.append( f"| {row['prediction_scale']:.2f} | {row['pds_l1_raw']:.6f} |" ) truth_matched = result["pds_scale_sensitivity"][ "norm_matched_truth_informed" ] lines.extend( [ "", f"Truth-norm-matched PDS: {truth_matched['pds_l1']:.6f}. This value uses paired observed norms, is diagnostic only, and is ineligible for hidden-test ranking.", "", "## Locked outputs", "", ] ) for metric in result["locked_metrics"]: lines.append(f"- `{metric['metric_id']}`: {metric['state']}.") lines.extend( [ "", "## Conclusion", "", result["conclusion"], "", "## Claim boundary", "", result["claim_boundary"], "", ] ) return "\n".join(lines) def run_self_test() -> None: perturbations = np.asarray(["A", "B", "C"]) response_genes = np.asarray(["A", "B", "C", "X"]) actual = np.asarray( [ [2.0, 0.0, 0.0, 0.0], [0.0, 2.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0], ] ) perfect_score, _, ranks = pds_l1( actual.copy(), actual, perturbations, response_genes ) if perfect_score != 1.0 or ranks != [1, 1, 1]: raise AssertionError("perfect PDS self-test failed") zero_score, zero_values, _ = pds_l1( np.zeros_like(actual), actual, perturbations, response_genes ) if not (0.0 <= zero_score <= 1.0 and len(zero_values) == 3): raise AssertionError("zero PDS self-test failed") if pearson(actual[0], actual[0]) != 1.0: raise AssertionError("Pearson self-test failed") if spearman(actual[0], actual[0]) != 1.0: raise AssertionError("Spearman self-test failed") if pearson(np.zeros(4), actual[0]) is not None: raise AssertionError("constant-correlation null self-test failed") scaled, factors = norm_match_predictions(actual * 0.25, actual) if not np.allclose(scaled, actual) or not np.allclose(factors, 4.0): raise AssertionError("norm matching self-test failed") print("VCC metric diagnostic self-test: 6/6 passed") def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser() parser.add_argument("--bundle", type=Path) parser.add_argument("--summary", type=Path) parser.add_argument("--target-manifest", type=Path) parser.add_argument("--output-json", type=Path) parser.add_argument("--output-tsv", type=Path) parser.add_argument("--report-md", type=Path) parser.add_argument("--self-test", action="store_true") return parser.parse_args() def main() -> None: args = parse_args() if args.self_test: run_self_test() return required: dict[str, Path | None] = { "--bundle": args.bundle, "--output-json": args.output_json, "--output-tsv": args.output_tsv, "--report-md": args.report_md, } missing = [flag for flag, value in required.items() if value is None] if missing: raise SystemExit(f"Missing required arguments: {', '.join(missing)}") result, target_rows = diagnostic_object( bundle_path=args.bundle, summary_path=args.summary, manifest_path=args.target_manifest, ) args.output_json.parent.mkdir(parents=True, exist_ok=True) args.output_json.write_text( json.dumps(result, indent=2, ensure_ascii=False) + "\n", encoding="utf-8", ) write_tsv(args.output_tsv, flatten_target_rows(target_rows)) args.report_md.parent.mkdir(parents=True, exist_ok=True) args.report_md.write_text(markdown_report(result), encoding="utf-8") print( f"Wrote VCC metric diagnostic for " f"{result['scientific_context']['n_external_targets']} targets" ) if __name__ == "__main__": main()