diff --git a/CHANGELOG.md b/CHANGELOG.md index a111086..88b88ff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,36 @@ # Changelog +## 1.4.4 - 2026-10-01 + +### Fixed +- Per-run CSV reuse now requires `interfaces.csv.meta.json`, matching scoring options, inputs, software/calibration, and CSV checksum. The default checks file size, mtime, ctime and file identity without reading input payloads; `--cache_validation content` hashes every input for strict validation. Metadata can miss same-size edits on filesystems with coarse or preserved timestamps. Old caches and incomplete scoring runs are recomputed; CSVs and manifests are published with atomic file replacement. +- Native AlphaFold-Multimer rankings with only `order` and `iptm+ptm` are accepted. Individual confidences and missing PAE JSON are recovered from result pickles, including gzip/xz; unavailable individual scores remain missing. +- Parser backend and selected structure are written to `backend` and `structure_file`. Calibration and report labels use the explicit backend, preventing Boltz filenames from selecting an AF2 ladder. Legacy CSV inference is retained. +- Corrupt JSON reports the actual file, PAE dimensions and values are checked before scoring, and unmappable Boltz PAE is rejected instead of trimmed. Failed models log tracebacks and cannot leave partially scored model rows or a reusable completed cache. +- AF2 structure choice is deterministic: exact preferred state, other state, model name, model directory, then ranked output; CIF precedes PDB within each category. `--af2_structure relaxed` is the default; `unrelaxed` reverses the first two choices. Legacy fallback matching respects model-name boundaries and sorts candidates. This can change scores when both coordinate sets exist: the old glob order was unspecified, and the frozen calibration predates per-row structure provenance, so its relaxed/unrelaxed selection cannot be reconstructed reliably. The existing calibration is retained, not represented as a uniformly relaxed or unrelaxed reference. +- PAE PNG requests and rendering failures no longer invalidate CSV scores. Missing, stale or failed plots can be regenerated on a cache hit without rerunning interface scoring or biophysical calculations. +- A missing AF3 full-confidence file now raises a file-not-found error even when the summary contains chain-pair minima. Existing minima-only files retain a distinct insufficient-PAE error; combined exports with actual full PAE remain supported. +- Filtered residues retain a full Bio.PDB hierarchy, so equal residue numbers on different chains have distinct identities. The existing pairwise pDockQ2 formula and coefficients are unchanged, with an IPSAE compatibility regression. +- Per-interface and aggregate ranking tables paginate without dropping rows, with repeated headers and sequential page numbers. Reports show the actual metascore feature count instead of implying that every row has all eleven features. +- Report imports preserve the caller's Matplotlib backend and figures; rendering restores plotting settings even on errors. Invalid `--aggregate_report` usage is rejected before scoring starts. +- CSVs record scoring cutoffs, a content-derived calibration ID, and the contributing metascore features/count. Reports disclose recalculation versus stored scores and flag custom or unknown cutoffs; custom settings retain the existing formulas and reference distributions, with percentiles labelled exploratory. The unsupported interface-area entry was removed from the percentile registry; raw area remains available. + +### CSV additions + +| Column | Meaning | +| --- | --- | +| `backend` | Prediction source: `af2`, `af3`, or `boltz2`. | +| `structure_file` | Selected structure path relative to the run directory. | +| `contact_thresh` | Residue contact cutoff in Å. | +| `pae_filter` | Maximum mean interface PAE accepted for a CSV row, in Å. | +| `ipsae_pae_cutoff` | PAE cutoff used by ipSAE, in Å. | +| `metascore_calibration` | Content-derived identifier of the frozen calibration. | +| `metascore_calibration_status` | Whether scoring cutoffs match the calibration defaults: `default`, `custom`, or `unknown`. | +| `metascore_features` | Semicolon-separated names of the features contributing to this row's metascore. | +| `metascore_feature_count` | Number of contributing features. | + +Custom parsers should accept `parse_run(directory, *, options=ParseOptions())`, ignoring unused options; `cache_options` identifies which parser options affect scores. + ## 1.4.3 - 2026-09-28 ### Fixed diff --git a/README.md b/README.md index cc20b7a..bbb3c2c 100644 --- a/README.md +++ b/README.md @@ -237,9 +237,8 @@ AlphaJudge writes `interfaces.csv` with one row per interface (and includes the - **interface_hb, interface_sb, interface_sc, interface_area, interface_solv_en** Exact header is asserted in tests to be consistent across AF2 and AF3 runs. -When a normal or summary run finds a cached per-run CSV from an older release that lacks current -required fields, AlphaJudge recomputes that run automatically; `--force_recompute` remains -available when values, rather than the schema, must be refreshed. +Cached results are reused when inputs and scoring settings match; use `--force_recompute` to rescore. +See the [changelog](CHANGELOG.md) for release details. --- diff --git a/pyproject.toml b/pyproject.toml index c6e039f..fe3df75 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "alphajudge" -version = "1.4.3" +version = "1.4.4" description = "Evaluate AlphaFold-predicted protein complexes using confidence metrics and interface biophysics." readme = { file = "README.md", content-type = "text/markdown" } requires-python = ">=3.10" diff --git a/src/alphajudge/cache.py b/src/alphajudge/cache.py new file mode 100644 index 0000000..0622065 --- /dev/null +++ b/src/alphajudge/cache.py @@ -0,0 +1,139 @@ +"""Versioned provenance for per-run CSV reuse. + +Input metadata is checked by default; content validation hashes compressed +inputs as stored. Package source (including calibration) and CSV output are +always hashed. Missing/old manifests are cache misses. +""" +from __future__ import annotations + +from contextlib import contextmanager +from dataclasses import asdict, dataclass +import hashlib +from importlib.metadata import PackageNotFoundError, version +import json +from pathlib import Path +import platform +import tempfile + +from .parsers import ParseOptions + +SCHEMA_VERSION = 2 +_INPUT_SUFFIXES = (".json", ".json.gz", ".json.xz", ".pkl", ".pkl.gz", ".pkl.xz", + ".npz", ".cif", ".pdb") + + +@dataclass(frozen=True) +class ScoringSettings: + contact_thresh: float = 8.0 + pae_filter: float = 100.0 + models_to_analyse: str = "best" + ipsae_pae_cutoff: float = 10.0 + skip_biophysical_scores: bool = False + parser_options: ParseOptions | None = None + + +def file_digest(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 manifest_path(csv_path: Path) -> Path: + return csv_path.with_name(csv_path.name + ".meta.json") + + +def file_fingerprint(path: Path) -> dict[str, int]: + """Cheap change detection, including file replacements. + + ctime supplements mtime, but coarse/preserved timestamps can still hide + same-size rewrites. Use content validation when metadata is insufficient. + """ + stat = path.stat() + return {"size": stat.st_size, "mtime_ns": stat.st_mtime_ns, "ctime_ns": stat.st_ctime_ns, + "device": stat.st_dev, "inode": stat.st_ino} + + +def software_identity() -> dict: + root = Path(__file__).parent + dependencies = {} + for package in ("alphajudge", "numpy", "scipy", "biopython", "matplotlib"): + try: + dependencies[package] = version(package) + except PackageNotFoundError: + dependencies[package] = "uninstalled" + return { + "python": platform.python_version(), "packages": dependencies, + "source": {str(p.relative_to(root)): file_digest(p) for p in sorted(root.rglob("*.py"))}, + } + + +def request_identity(directory: Path, csv_path: Path, *, + settings: ScoringSettings = ScoringSettings(), + cache_validation: str = "stat") -> dict: + if cache_validation not in {"stat", "content"}: + raise ValueError("cache_validation must be 'stat' or 'content'") + inputs = {} + for path in sorted(directory.rglob("*")): + if not path.is_file() or path == csv_path or path.name.endswith(".meta.json"): + continue + if path.name.endswith(_INPUT_SUFFIXES) or path.name.endswith("ranking_scores.csv"): + inputs[str(path.relative_to(directory))] = ( + file_digest(path) if cache_validation == "content" else file_fingerprint(path) + ) + return {"schema": SCHEMA_VERSION, "directory": str(directory.resolve()), + "cache_validation": cache_validation, + "settings": asdict(settings), "inputs": inputs, + "software": software_identity()} + + +@contextmanager +def atomic_text(path: Path): + """Publish a complete file with a same-directory atomic rename.""" + path.parent.mkdir(parents=True, exist_ok=True) + with tempfile.NamedTemporaryFile(mode="w", newline="", dir=path.parent, + prefix=f".{path.name}.", suffix=".tmp", delete=False) as handle: + temporary = Path(handle.name) + try: + yield handle + handle.close() + temporary.replace(path) + finally: + temporary.unlink(missing_ok=True) + + +def write_manifest(csv_path: Path, request: dict, *, csv_sha256: str, complete: bool, + backend: str, models: list[str], structure_files: dict[str, str], + pae_pngs: dict | None = None) -> None: + manifest = {"request": request, "complete": complete, + "csv_sha256": csv_sha256, "backend": backend, + "models": models, "structure_files": structure_files, "pae_pngs": pae_pngs or {}} + _write_manifest_data(csv_path, manifest) + + +def _write_manifest_data(csv_path: Path, manifest: dict) -> None: + with atomic_text(manifest_path(csv_path)) as handle: + json.dump(manifest, handle, indent=2, sort_keys=True) + handle.write("\n") + + +def read_manifest(csv_path: Path) -> dict: + return json.loads(manifest_path(csv_path).read_text()) + + +def update_pae_pngs(csv_path: Path, request: dict, pngs: dict) -> None: + """Update optional plots without certifying new CSV bytes or scoring options.""" + manifest = read_manifest(csv_path) + if manifest.get("request") == request: + manifest["pae_pngs"] = pngs + _write_manifest_data(csv_path, manifest) + + +def matches(csv_path: Path, request: dict) -> bool: + try: + manifest = read_manifest(csv_path) + return (manifest.get("complete") is True and manifest.get("request") == request + and manifest.get("csv_sha256") == file_digest(csv_path)) + except (OSError, ValueError, AttributeError): + return False diff --git a/src/alphajudge/cli.py b/src/alphajudge/cli.py index 124fd1f..43eae6e 100644 --- a/src/alphajudge/cli.py +++ b/src/alphajudge/cli.py @@ -14,6 +14,10 @@ def main() -> None: p.add_argument("--pae_filter", type=float, default=100.0) p.add_argument("--ipsae_pae_cutoff", type=float, default=10.0) p.add_argument("--models_to_analyse", choices=["best","all"], default="best") + p.add_argument("--af2_structure", choices=["relaxed", "unrelaxed"], default="relaxed", + help="Preferred AF2 coordinates when both states exist (default: relaxed).") + p.add_argument("--cache_validation", choices=["stat", "content"], default="stat", + help="Validate input metadata (fast default) or hash full input content.") p.add_argument("-r","--recursive", action="store_true", help="Recursively search for runs under given PATHS") p.add_argument("-o","--summary", help="Write aggregated CSV across runs to this path") p.add_argument( @@ -67,6 +71,8 @@ def main() -> None: logging.basicConfig(level=logging.INFO, format="%(levelname)s:%(message)s") if not args.paths: p.error("Provide PATHS") + if args.aggregate_report and not args.summary: + p.error("--aggregate_report requires --summary") write_per_run_report = args.report if write_per_run_report is None: @@ -86,11 +92,11 @@ def main() -> None: skip_pae_png=args.skip_pae_png, skip_biophysical_scores=args.skip_biophysical_scores, write_per_run_report=write_per_run_report, + cache_validation=args.cache_validation, + af2_structure=args.af2_structure, ) if args.aggregate_report: - if not args.summary: - p.error("--aggregate_report requires --summary") if summary_path is None: p.error( f"--aggregate_report requested but no summary was written to " diff --git a/src/alphajudge/complex.py b/src/alphajudge/complex.py index deeff7d..4b4f3b9 100644 --- a/src/alphajudge/complex.py +++ b/src/alphajudge/complex.py @@ -5,11 +5,11 @@ from functools import cached_property from typing import Any -from Bio.PDB import Chain +from Bio.PDB import Chain, Model, Structure -from .confidence import Confidence +from .confidence import Confidence, validate_pae from .docking_scores import MPDOCKQ -from .geometry import is_pae_token_residue +from .geometry import chain_boundaries, is_pae_token_residue from .interface import Interface logger = logging.getLogger(__name__) @@ -35,7 +35,13 @@ def __init__( self.pae_filter = float(pae_filter) self.ipsae_pae_cutoff = None if ipsae_pae_cutoff is None else float(ipsae_pae_cutoff) - self._res_index_map, self._chain_indices_by_id, self._chains = self._build_maps() + # Retain the Bio.PDB hierarchy so equally numbered residues on distinct + # chains retain distinct identities after filtering. + model = next(structure.get_models()) + self._scored_structure = Structure.Structure(structure.id) + scored_model = Model.Model(model.id) + self._scored_structure.add(scored_model) + self._res_index_map, self._chain_indices_by_id, self._chains = self._build_maps(model, scored_model) self.interfaces: list[Interface] = [] for i in range(len(self._chains)): @@ -44,8 +50,7 @@ def __init__( if iface.num_intf_residues > 0: self.interfaces.append(iface) - def _build_maps(self) -> tuple[dict[tuple[str, Any], int], dict[str, list[int]], list[Chain.Chain]]: - model = next(self.structure.get_models()) + def _build_maps(self, model, scored_model) -> tuple[dict[tuple[str, Any], int], dict[str, list[int]], list[Chain.Chain]]: chains = list(model.get_chains()) res_index_map: dict[tuple[str, Any], int] = {} @@ -59,6 +64,7 @@ def _build_maps(self) -> tuple[dict[tuple[str, Any], int], dict[str, list[int]], continue new_chain = Chain.Chain(chain.id) + scored_model.add(new_chain) for residue in kept: new_chain.add(residue.copy()) filtered_chains.append(new_chain) @@ -70,12 +76,7 @@ def _build_maps(self) -> tuple[dict[tuple[str, Any], int], dict[str, list[int]], idx += 1 chain_indices_by_id[new_chain.id] = idxs - pae_n = len(self.conf.pae_matrix) - if idx != pae_n: - logger.warning( - f"token residues counted = {idx}, but PAE is {pae_n}x{pae_n}. " - f"Indexing may be misaligned for this structure." - ) + validate_pae(self.conf.pae_matrix, idx) return res_index_map, chain_indices_by_id, filtered_chains @@ -86,8 +87,7 @@ def num_chains(self) -> int: @property def chain_boundaries(self) -> list[float]: """PAE-matrix positions between consecutive chains (heatmap separators).""" - ends = sorted(max(idxs) + 0.5 for idxs in self._chain_indices_by_id.values() if idxs) - return ends[:-1] + return chain_boundaries(self._chain_indices_by_id) @cached_property def average_interface_plddt(self) -> float: diff --git a/src/alphajudge/confidence.py b/src/alphajudge/confidence.py index c6e5c7e..c6b9704 100644 --- a/src/alphajudge/confidence.py +++ b/src/alphajudge/confidence.py @@ -19,6 +19,17 @@ IPTM_SCOPE_GLOBAL = "global" +def validate_pae(pae: np.ndarray, residue_count: int, *, source: str = "PAE") -> None: + """Reject incomplete or invalid residue matrices before interface scoring.""" + if pae.shape != (residue_count, residue_count): + raise ValueError( + f"{source}: PAE shape {pae.shape} does not match {residue_count} scored residues " + f"(expected {residue_count}x{residue_count}); cannot align PAE to structure" + ) + if not np.all(np.isfinite(pae)) or np.any(pae < 0): + raise ValueError(f"{source}: PAE values must be finite and non-negative") + + @dataclass(frozen=True) class Confidence: pae_matrix: np.ndarray # Keep as numpy array for memory efficiency (7-10x reduction) diff --git a/src/alphajudge/geometry.py b/src/alphajudge/geometry.py index eececa7..3498d35 100644 --- a/src/alphajudge/geometry.py +++ b/src/alphajudge/geometry.py @@ -1,5 +1,7 @@ from __future__ import annotations +from collections.abc import Mapping, Sequence + from Bio.PDB.Polypeptide import is_aa @@ -21,6 +23,11 @@ } +def chain_boundaries(chain_indices: Mapping[str, Sequence[int]]) -> list[float]: + """PAE-matrix positions between nonempty scored chains.""" + return sorted(max(indices) + .5 for indices in chain_indices.values() if indices)[:-1] + + def representative_atom(res): """ Representative atom for distance / bfactor use. diff --git a/src/alphajudge/interface.py b/src/alphajudge/interface.py index a427d4f..8562c65 100644 --- a/src/alphajudge/interface.py +++ b/src/alphajudge/interface.py @@ -345,10 +345,7 @@ def _avg_plddt_union(self) -> float: def _avg_pae_over_pairs(self) -> float: vals = [] for i, j in self._pair_indices(): - try: - vals.extend((float(self._pae[i, j]), float(self._pae[j, i]))) - except IndexError: # PAE smaller than the scored residues - continue + vals.extend((float(self._pae[i, j]), float(self._pae[j, i]))) return sum(vals) / len(vals) if vals else float("nan") def _frac(self, names: set[str]) -> float: diff --git a/src/alphajudge/meta_score.py b/src/alphajudge/meta_score.py index fece51a..f2d0778 100644 --- a/src/alphajudge/meta_score.py +++ b/src/alphajudge/meta_score.py @@ -1,5 +1,7 @@ from __future__ import annotations +import hashlib +import json import math from bisect import bisect_right from collections.abc import Mapping @@ -31,7 +33,6 @@ "pDockQ/mpDockQ": 1.0, "interface_sc": 1.0, "interface_hb": 1.0, - "interface_area": 1.0, "interface_solv_en": -1.0, "interface_contact_prob_top10_mean": 1.0, "interface_ccc": 1.0, @@ -52,8 +53,8 @@ ) # Frozen deciles from the AlphaJudge interacting reference set. Calibrated on -# POSITIVE (interacting) pairs ONLY: 3,878 AF2/AF3 positive rows out of the -# 7,756-row balanced table. The database-negative re-pairings are deliberately +# POSITIVE (interacting) pairs ONLY: 12,163 AF2/AF3 positive rows from the full +# non-downsampled v3 table. The database-negative re-pairings are deliberately # excluded so a new prediction is ranked against the distribution of real # interfaces, not against a 50% non-interacting decoy population. Regenerate # manually with @@ -548,6 +549,20 @@ } +# Content-addressed identity travels with standalone CSVs, independently of +# the run cache. Include the feature/scope policy as well as the frozen data. +CALIBRATION_ID = "aj-positive-v3-" + hashlib.sha256(json.dumps({ + "policy": "available-mean_scope-v1_explicit-backend-v1", + "features": META_SCORE_FEATURES, + "directions": FEATURE_DIRECTIONS, + "levels": CALIBRATION_LEVELS, + "pooled": BENCHMARK_QUANTILES, + "backends": BENCHMARK_QUANTILES_BY_BACKEND, +}, sort_keys=True, separators=(",", ":")).encode()).hexdigest()[:12] + +_REFERENCE_PARAMETERS = {"contact_thresh": 8., "pae_filter": 100., "ipsae_pae_cutoff": 10.} + + def _safe_float(value: Any) -> float: try: parsed = float(value) @@ -556,6 +571,21 @@ def _safe_float(value: Any) -> float: return parsed if math.isfinite(parsed) else float("nan") +def calibration_parameter_status(row: Mapping[str, Any]) -> str: + """Flag custom or unknown settings; this does not refit the calibration. + + A custom PAE row filter changes selection rather than the per-row scores, + but still changes the population being compared with the reference set. + """ + values = {key: _safe_float(row.get(key)) for key in _REFERENCE_PARAMETERS} + if any(math.isfinite(value) and value != _REFERENCE_PARAMETERS[key] + for key, value in values.items()): + return "custom" + if any(not math.isfinite(value) for value in values.values()): + return "unknown" + return "default" + + def calibrated_feature_percentile( feature: str, value: Any, backend: str | None = None ) -> float | None: @@ -593,6 +623,10 @@ def calibrated_feature_percentile( def infer_backend(row: Mapping[str, Any]) -> str | None: """Best-effort backend tag for a scored row ("af2"/"af3", else None).""" + backend = str(row.get("backend") or "").strip().lower() + if backend: + # An explicit uncalibrated backend (e.g. Boltz) uses the pooled ladder. + return backend if backend in {"af2", "af3"} else None src = str(row.get("interface_contact_prob_source") or "") if src.startswith("af2"): return "af2" @@ -620,6 +654,19 @@ def feature_is_comparable(row: Mapping[str, Any], feature: str) -> bool: return feature != "iptm" or row.get("iptm_scope") == IPTM_SCOPE_CHAIN_PAIR +def meta_score_components(row: Mapping[str, Any]) -> dict[str, float]: + """The available, scope-compatible percentile contributions for this row.""" + backend = infer_backend(row) + percentiles = {} + for feature in META_SCORE_FEATURES: + if not feature_is_comparable(row, feature): + continue + percentile = calibrated_feature_percentile(feature, row.get(feature), backend) + if percentile is not None: + percentiles[feature] = percentile + return percentiles + + def interface_meta_score(row: Mapping[str, Any]) -> float: """ Transparent rank-style interface metascore. @@ -629,14 +676,7 @@ def interface_meta_score(row: Mapping[str, Any]) -> float: non-finite and scope-incompatible inputs are ignored. The final score is the mean of the remaining percentiles; no AF3x-specific calibration is implied. """ - backend = infer_backend(row) - percentiles = [] - for feature in META_SCORE_FEATURES: - if not feature_is_comparable(row, feature): - continue - percentile = calibrated_feature_percentile(feature, row.get(feature), backend) - if percentile is not None: - percentiles.append(percentile) + percentiles = meta_score_components(row) if not percentiles: return float("nan") - return float(sum(percentiles) / len(percentiles)) + return float(sum(percentiles.values()) / len(percentiles)) diff --git a/src/alphajudge/parsers/__init__.py b/src/alphajudge/parsers/__init__.py index b2ca61f..0ee7328 100644 --- a/src/alphajudge/parsers/__init__.py +++ b/src/alphajudge/parsers/__init__.py @@ -1,22 +1,33 @@ from __future__ import annotations from abc import ABC, abstractmethod from collections.abc import Callable -from dataclasses import dataclass +from dataclasses import dataclass, field from pathlib import Path from typing import Any import gzip import json import lzma +import re from Bio.PDB import PDBParser, MMCIFParser from ..confidence import Confidence from ..geometry import is_pae_token_residue, representative_atom +@dataclass(frozen=True) +class ParseOptions: + structure_preference: str = "relaxed" + + def __post_init__(self): + if self.structure_preference not in {"relaxed", "unrelaxed"}: + raise ValueError("Structure preference must be 'relaxed' or 'unrelaxed'") + + @dataclass class Run: order: list[str] source: str load_model: Callable[[str], tuple[Any, Confidence]] + structure_files: dict[str, str] = field(default_factory=dict) class BaseParser(ABC): @@ -26,7 +37,11 @@ class BaseParser(ABC): @abstractmethod def detect(d: Path) -> bool: ... @abstractmethod - def parse_run(self, d: Path) -> Run: ... + def parse_run(self, d: Path, *, options: ParseOptions = ParseOptions()) -> Run: ... + + def cache_options(self, options: ParseOptions) -> ParseOptions | None: + """Only options used by this parser belong in the scoring cache.""" + return None # Compression magic numbers (file header bytes). _MAGIC = ((b"\xfd7zXZ\x00", lzma.open), (b"\x1f\x8b", gzip.open)) @@ -47,7 +62,7 @@ def _open_maybe_compressed(cls, p: Path, mode: str = "rt"): return p.open(mode) @classmethod - def _read_json(cls, p: Path) -> dict: + def _read_json(cls, p: Path, *, required: bool = False) -> Any: """Read a JSON file, transparently handling xz/gz compression. AlphaPulldown's ``--storage_mode slim/minimal`` may store large JSON @@ -56,17 +71,18 @@ def _read_json(cls, p: Path) -> dict: scoring is unaffected by the storage mode. Compression is detected from the file's magic bytes, not its extension. """ + target = cls._first_existing( + [p, p.with_name(p.name + ".xz"), p.with_name(p.name + ".gz")] + ) + if target is None: + if required: + raise FileNotFoundError(f"Required JSON file not found: {p} (or .xz/.gz sibling)") + return {} try: - # Plain path absent: try a compressed sibling from slim/minimal. - target = cls._first_existing( - [p, p.with_name(p.name + ".xz"), p.with_name(p.name + ".gz")] - ) - if target is None: - return {} with cls._open_maybe_compressed(target) as fh: return json.load(fh) - except Exception: - return {} + except (OSError, EOFError, ValueError, lzma.LZMAError) as exc: + raise ValueError(f"Could not read JSON {target}: {exc}") from exc @staticmethod def _first_existing(paths) -> Path | None: @@ -79,14 +95,32 @@ def _load_structure(path: str): return parser.get_structure("complex", str(p)) @staticmethod - def _guess_struct(d: Path, model: str) -> str: - if (d / model / "model.cif").exists(): - return str(d / model / "model.cif") - for ext in ("cif", "pdb"): - hits = list(d.glob(f"*{model}*.{ext}")) - if hits: - return str(hits[0]) - raise ValueError(f"struct for model {model} not found") + def _guess_struct(d: Path, model: str, rank: int | None = None, *, options: ParseOptions = ParseOptions()) -> str: + """AF2 precedence: preferred state, other state, exact model, directory, rank. + + Prefer CIF to PDB *within* each category. The preferred state's PDB + wins over the other state's CIF. Decorated names are a sorted fallback. + """ + states = (("relaxed_", "unrelaxed_") if options.structure_preference == "relaxed" + else ("unrelaxed_", "relaxed_")) + candidates = [d / f"{prefix}{model}.{ext}" + for prefix in (*states, "") + for ext in ("cif", "pdb")] + candidates.extend(d / model / f"model.{ext}" for ext in ("cif", "pdb")) + if rank is not None: + candidates.extend(d / f"ranked_{rank}.{ext}" for ext in ("cif", "pdb")) + exact = BaseParser._first_existing(candidates) + if exact is not None: + return str(exact) + pattern = re.compile(r"(?:^|[_\-.])" + re.escape(model) + r"(?:$|[_\-.])") + hits = [p for p in d.iterdir() if p.suffix in (".cif", ".pdb") + and pattern.search(p.stem)] + if hits: + def priority(p): + kind = next((i for i, state in enumerate(states) if p.name.startswith(state)), 2) + return kind, p.suffix != ".cif", p.name + return str(sorted(hits, key=priority)[0]) + raise ValueError(f"structure for model {model} not found in {d}") @staticmethod def _maps(struct): diff --git a/src/alphajudge/parsers/af2.py b/src/alphajudge/parsers/af2.py index 731fd14..2771f87 100644 --- a/src/alphajudge/parsers/af2.py +++ b/src/alphajudge/parsers/af2.py @@ -1,10 +1,11 @@ from __future__ import annotations import logging import pickle +import re from pathlib import Path import numpy as np -from . import BaseParser, Run -from ..confidence import Confidence +from . import BaseParser, ParseOptions, Run +from ..confidence import Confidence, validate_pae from ..contact_probs import ( AF2_DISTOGRAM_CONTACT_CUTOFF, contact_probs_from_distogram, @@ -17,28 +18,53 @@ class AF2Parser(BaseParser): _warned_missing_distogram = False def detect(self, d: Path) -> bool: - return (d / "ranking_debug.json").exists() + return any((d / f"ranking_debug.json{suffix}").exists() for suffix in ("", ".xz", ".gz")) - def parse_run(self, d: Path) -> Run: - rj = self._read_json(d / "ranking_debug.json") + def cache_options(self, options: ParseOptions) -> ParseOptions: + return options + + def parse_run(self, d: Path, *, options: ParseOptions = ParseOptions()) -> Run: + rj = self._read_json(d / "ranking_debug.json", required=True) order = rj["order"] + structure_files = {} def load_model(model: str): - struct = self._load_structure(self._guess_struct(d, model)) + structure_path = Path(self._guess_struct(d, model, order.index(model), options=options)) + struct = self._load_structure(structure_path) + structure_files[model] = str(structure_path.relative_to(d)) chains, rim, _ = self._maps(struct) - # AF2: full residue×residue matrix in pae_{model}.json - pae_payload = self._read_json(d / f"pae_{model}.json") - pae = np.array(pae_payload[0]["predicted_aligned_error"], dtype=float) - max_pae = float(np.nanmax(pae) if pae.size else np.nan) - contact_probs = self._load_contact_probs_from_result_pkl(d, model, pae.shape) + result_path, result = self._load_result(d, model) + pae_path = d / f"pae_{model}.json" + json_path = self._first_existing([pae_path.with_name(pae_path.name + s) for s in ("", ".xz", ".gz")]) + if json_path is not None: + payload = self._read_json(json_path, required=True) + try: + raw_pae = (payload[0] if isinstance(payload, list) else payload)["predicted_aligned_error"] + except (KeyError, IndexError, TypeError) as exc: + raise ValueError(f"{json_path}: missing predicted_aligned_error matrix") from exc + pae_source = str(json_path) + elif "predicted_aligned_error" in result: + raw_pae = result["predicted_aligned_error"] + pae_source = str(result_path) + else: + raise ValueError(f"No full PAE matrix for {model}: expected {pae_path} or predicted_aligned_error in result pickle") + try: + pae = np.asarray(raw_pae, dtype=float) + except (TypeError, ValueError) as exc: + raise ValueError(f"{pae_source}: invalid PAE matrix: {exc}") from exc + validate_pae(pae, len(rim), source=pae_source) + max_pae = float(np.max(pae)) if pae.size else float("nan") + contact_probs = self._contact_probs_from_result(result, result_path, pae.shape) # AF2 rankings - is_multimer = ("iptm+ptm" in rj) and ("iptm" in rj) + is_multimer = "iptm+ptm" in rj or "iptm" in rj or "iptm" in result + def score(name): + value = self._safe_float(rj.get(name, {}).get(model)) + return value if value is not None else self._safe_float(result.get(name)) + if is_multimer: - iptm = self._safe_float(rj["iptm"].get(model)) - ptm = self._safe_float(rj.get("ptm", {}).get(model)) - iptm_ptm = self._safe_float(rj["iptm+ptm"].get(model)) + iptm, ptm, iptm_ptm = score("iptm"), score("ptm"), score("iptm+ptm") # Backfill when PTM is not provided in AF2 multimer JSON if ptm is None and (iptm_ptm is not None) and (iptm is not None): ptm = (iptm_ptm - 0.8 * iptm) / 0.2 @@ -47,7 +73,7 @@ def load_model(model: str): iptm_ptm = 0.8 * iptm + 0.2 * ptm conf = iptm_ptm else: - iptm, ptm = 0.0, self._safe_float(rj["ptm"][model]) + iptm, ptm = 0.0, score("ptm") iptm_ptm = conf = ptm plddt = self._plddt(chains, rim) @@ -62,24 +88,38 @@ def load_model(model: str): else None ), ) - return Run(order=order, source="af2", load_model=load_model) + return Run(order=order, source="af2", load_model=load_model, structure_files=structure_files) @classmethod def _load_contact_probs_from_result_pkl( cls, d: Path, model: str, expected_shape: tuple[int, int] ) -> np.ndarray | None: + result_pkl, payload = cls._load_result(d, model) + return cls._contact_probs_from_result(payload, result_pkl, expected_shape) + + @classmethod + def _load_result(cls, d: Path, model: str) -> tuple[Path | None, dict]: result_pkl = cls._find_result_pkl(d, model) if result_pkl is None: - return None + return None, {} try: with cls._open_maybe_compressed(result_pkl, "rb") as f: payload = pickle.load(f) except Exception as e: logger.warning(f"could not read AF2 result pickle {result_pkl}: {e}") - return None + return result_pkl, {} if not isinstance(payload, dict): + logger.warning("AF2 result pickle %s does not contain a mapping", result_pkl) + return result_pkl, {} + return result_pkl, payload + + @classmethod + def _contact_probs_from_result( + cls, payload: dict, result_pkl: Path | None, expected_shape: tuple[int, int] + ) -> np.ndarray | None: + if result_pkl is None or not payload: return None distogram = payload.get("distogram") if not isinstance(distogram, dict): @@ -131,7 +171,8 @@ def _find_result_pkl(cls, d: Path, model: str) -> Path | None: candidates = [ stem.with_name(stem.name + suffix) for stem in stems for suffix in ("", ".gz", ".xz") ] - candidates.extend(sorted(d.glob(f"result*{model}*.pkl*"))) + pattern = re.compile(r"(?:^|[_\-.])" + re.escape(model) + r"(?:$|[_\-.])") + candidates.extend(p for p in sorted(d.glob("result*.pkl*")) if pattern.search(p.name)) if (d / model).is_dir(): candidates.extend(sorted((d / model).glob("result*.pkl*"))) return cls._first_existing(candidates) diff --git a/src/alphajudge/parsers/af3.py b/src/alphajudge/parsers/af3.py index 8384af8..1a06d96 100644 --- a/src/alphajudge/parsers/af3.py +++ b/src/alphajudge/parsers/af3.py @@ -5,7 +5,7 @@ import logging import re import numpy as np -from . import BaseParser, Run +from . import BaseParser, ParseOptions, Run from ..confidence import SCOPE_INCLUDES_EXCLUDED_TOKENS, SCOPE_SCORED_RESIDUES, Confidence from ..geometry import is_pae_token_residue @@ -42,28 +42,34 @@ class AF3Parser(BaseParser): def detect(d: Path) -> bool: return AF3Parser._ranking_scores_file(d) is not None - def parse_run(self, d: Path) -> Run: + def parse_run(self, d: Path, *, options: ParseOptions = ParseOptions()) -> Run: ranking_file = self._ranking_scores_file(d) if ranking_file is None: raise ValueError(f"AF3 ranking scores file not found in {d}") order, ranking_scores = self._read_csv_order(ranking_file) job_prefix = self._job_prefix_from_ranking_file(ranking_file) + structure_files = {} def load_model(model: str): is_best_model = bool(order and model == order[0]) - struct = self._load_structure( - self._guess_af3_struct(d, model, job_prefix, is_best_model) - ) + structure_path = Path(self._guess_af3_struct(d, model, job_prefix, is_best_model)) + struct = self._load_structure(structure_path) + structure_files[model] = str(structure_path.relative_to(d)) chains, rim, cid = self._maps(struct) - summary = self._read_json( - self._find_af3_json( - d, model, "summary_confidences", job_prefix, is_best_model - ) - ) - matrix = self._read_json( - self._find_af3_json(d, model, "confidences", job_prefix, is_best_model) - ) or summary + summary_path = self._find_af3_json(d, model, "summary_confidences", job_prefix, is_best_model) + summary = self._read_json(summary_path) + if not isinstance(summary, dict): + raise ValueError(f"{summary_path}: expected a JSON object") + matrix_path = self._find_af3_json(d, model, "confidences", job_prefix, is_best_model) + if not matrix_path.exists() and any(summary.get(key) is not None for key in ("pae", "predicted_aligned_error")): + # Some combined exports carry the actual full matrix here. + # Chain-pair minima alone cannot replace a missing PAE file. + matrix, matrix_path = summary, summary_path + else: + matrix = self._read_json(matrix_path, required=True) + if not isinstance(matrix, dict): + raise ValueError(f"{matrix_path}: expected a JSON object carrying full PAE") iptm = self._safe_float(summary.get("iptm")) ptm = self._safe_float(summary.get("ptm")) @@ -90,7 +96,10 @@ def load_model(model: str): logger.warning("chain_pair_iptm dimensions do not match the source chains; skipping pair ipTM.") chain_pair_iptm = None - pae, max_pae, residue_tokens = self._normalize_pae_af3(matrix, chains, cid) + try: + pae, max_pae, residue_tokens = self._normalize_pae_af3(matrix, chains, cid) + except (ValueError, TypeError) as exc: + raise ValueError(f"{matrix_path}: {exc}") from exc contact_probs = self._normalize_contact_probs_af3(matrix, residue_tokens) plddt = self._plddt(chains, rim) global_scope = ( @@ -108,7 +117,7 @@ def load_model(model: str): chain_pair_iptm_chain_ids=summary_chain_ids, global_confidence_scope=global_scope, ) - return Run(order=order, source="af3", load_model=load_model) + return Run(order=order, source="af3", load_model=load_model, structure_files=structure_files) # ---- AF3-specific helpers ---- @staticmethod diff --git a/src/alphajudge/parsers/boltz.py b/src/alphajudge/parsers/boltz.py index b3e09fe..0755df6 100644 --- a/src/alphajudge/parsers/boltz.py +++ b/src/alphajudge/parsers/boltz.py @@ -9,8 +9,8 @@ import numpy as np from Bio.PDB.Polypeptide import is_aa -from . import BaseParser, Run -from ..confidence import Confidence +from . import BaseParser, ParseOptions, Run +from ..confidence import Confidence, validate_pae from ..geometry import is_pae_token_residue logger = logging.getLogger(__name__) @@ -37,7 +37,7 @@ class Boltz2Parser(BaseParser): def detect(d: Path) -> bool: return bool(Boltz2Parser._model_entries(d)) - def parse_run(self, d: Path) -> Run: + def parse_run(self, d: Path, *, options: ParseOptions = ParseOptions()) -> Run: entries = self._model_entries(d) if not entries: raise ValueError(f"Boltz-2 prediction files not found in {d}") @@ -49,7 +49,7 @@ def load_model(model: str): struct = self._load_structure(str(entry.structure_file)) chains, rim, cid = self._maps(struct) - summary = self._read_json(entry.confidence_file) + summary = self._read_json(entry.confidence_file, required=True) token_count, token_index = self._residue_tokens(chains, cid) pae, max_pae = self._load_pae(entry.pae_file, len(rim), token_count, token_index) plddt = ( @@ -87,7 +87,8 @@ def load_model(model: str): chain_pair_iptm_chain_ids=chain_ids, ) - return Run(order=order, source="boltz2", load_model=load_model) + return Run(order=order, source="boltz2", load_model=load_model, + structure_files={e.name: str(e.structure_file.relative_to(d)) for e in entries}) @staticmethod def _model_entries(d: Path) -> list[_BoltzModel]: @@ -124,12 +125,14 @@ def _rank_from_model_name(model_name: str) -> int: return int(match.group(1)) @staticmethod - def _load_npz_array(path: Path, preferred_key: str) -> np.ndarray | None: + def _load_npz_array(path: Path, preferred_key: str, *, required: bool = False) -> np.ndarray | None: try: with np.load(path) as payload: key = preferred_key if preferred_key in payload else payload.files[0] return np.array(payload[key], dtype=float) except Exception as exc: + if required: + raise ValueError(f"Could not read {preferred_key} array {path}: {exc}") from exc logger.warning(f"could not read {path}: {exc}") return None @@ -167,28 +170,19 @@ def _load_pae( if pae_file is None: return pae, 100.0 - matrix = cls._load_npz_array(pae_file, "pae") - if matrix is None or not matrix.size: - return pae, 100.0 + matrix = cls._load_npz_array(pae_file, "pae", required=True) if token_index is not None and matrix.shape == (token_count, token_count): selected = matrix[np.ix_(token_index, token_index)] + validate_pae(selected, total, source=str(pae_file)) return selected, float(np.nanmax(selected)) if selected.size else 100.0 if matrix.shape == pae.shape: + validate_pae(matrix, total, source=str(pae_file)) return matrix, float(np.nanmax(matrix)) - if matrix.ndim == 2 and matrix.shape[0] >= total and matrix.shape[1] >= total: - logger.warning( - f"Cannot map Boltz-2 PAE shape {matrix.shape} to the structure's tokens; " - f"using its leading {total}x{total} block." - ) - trimmed = matrix[:total, :total] - return trimmed, float(np.nanmax(trimmed)) if trimmed.size else 100.0 - - logger.warning( - f"Boltz-2 PAE shape {matrix.shape} != expected {pae.shape}; " - "using default PAE=100 for all residue pairs." + raise ValueError( + f"{pae_file}: Cannot align Boltz-2 PAE shape {matrix.shape} with " + f"{total} scored residues / {token_count} structure tokens" ) - return pae, float(np.nanmax(matrix)) @classmethod def _load_plddt( diff --git a/src/alphajudge/report.py b/src/alphajudge/report.py index c6506a5..dcc20a0 100644 --- a/src/alphajudge/report.py +++ b/src/alphajudge/report.py @@ -18,35 +18,43 @@ import logging import math import statistics +import textwrap from collections import Counter from collections.abc import Mapping, Sequence +from contextlib import contextmanager from datetime import datetime from pathlib import Path from typing import Any import matplotlib -matplotlib.use("Agg") import matplotlib.image as mpimg -import matplotlib.pyplot as plt import numpy as np +from matplotlib.backends.backend_agg import FigureCanvasAgg from matplotlib.backends.backend_pdf import PdfPages from matplotlib.colors import LinearSegmentedColormap +from matplotlib.figure import Figure from matplotlib.patches import Rectangle from matplotlib.ticker import FuncFormatter, MaxNLocator from .confidence import SCOPE_INCLUDES_EXCLUDED_TOKENS from .meta_score import ( + CALIBRATION_ID, META_SCORE_FEATURES, + calibration_parameter_status, calibrated_feature_percentile, feature_is_comparable, infer_backend, interface_meta_score, + meta_score_components, ) logger = logging.getLogger(__name__) _A4 = (8.27, 11.69) +_INTERFACE_ROWS_PER_PAGE = 28 +_AGGREGATE_COVER_ROWS = 10 +_AGGREGATE_ROWS_PER_PAGE = 30 # Percentile graphic: red -> pale centre -> blue. _SLIDER_CMAP = LinearSegmentedColormap.from_list( @@ -87,6 +95,8 @@ _GRADIENT = np.tile(np.linspace(0.0, 1.0, 1024), (2, 1)) _FEATURE_DISPLAY = { + "interface_LIS": "LIS", + "interface_pDockQ2": "pDockQ2", "interface_contact_prob_top10_mean": "Contact probability", "interface_ccc": "Confident contacts", "interface_ipSAE": "Interface ipSAE", @@ -145,8 +155,9 @@ # style + utility helpers # --------------------------------------------------------------------------- -def _setup_rcparams() -> None: - """Use a Computer-Modern-like serif PDF look, close to wwPDB reports.""" +@contextmanager +def _report_style(): + """Apply the report style temporarily, including on failed renders.""" rcparams = { "font.family": "serif", "font.serif": [ @@ -174,7 +185,8 @@ def _setup_rcparams() -> None: "ps.fonttype": 42, "axes.unicode_minus": False, } - matplotlib.rcParams.update(rcparams) + with matplotlib.rc_context(rcparams): + yield def _safe_float(value: Any) -> float | None: @@ -260,6 +272,10 @@ def _shorten_path(path: str, max_len: int = 64) -> str: def _detect_backend(rows: Sequence[Mapping[str, Any]]) -> str: + explicit = {str(r.get("backend") or "").strip().lower() for r in rows} - {""} + if explicit: + names = {"af2": "AlphaFold 2", "af3": "AlphaFold 3", "boltz2": "Boltz-2"} + return ", ".join(names.get(backend, backend) for backend in sorted(explicit)) for r in rows: model = str(r.get("model_used") or "") if "multimer" in model.lower(): @@ -301,11 +317,12 @@ def _decile_label(pct: float | None) -> str: # page primitives # --------------------------------------------------------------------------- -def _new_figure() -> plt.Figure: - return plt.figure(figsize=_A4, facecolor="white") +def _new_figure() -> Figure: + # Unmanaged figures do not touch pyplot's backend or figure registry. + return Figure(figsize=_A4, facecolor="white") -def _text_axes(fig: plt.Figure, rect: tuple[float, float, float, float]): +def _text_axes(fig: Figure, rect: tuple[float, float, float, float]): """Invisible unit-square axes for placing text and shapes in ``transAxes``.""" ax = fig.add_axes(rect) ax.set_xlim(0, 1) @@ -314,7 +331,7 @@ def _text_axes(fig: plt.Figure, rect: tuple[float, float, float, float]): return ax -def _label_axes(fig: plt.Figure, rect: tuple[float, float, float, float], text: str, **kwargs) -> None: +def _label_axes(fig: Figure, rect: tuple[float, float, float, float], text: str, **kwargs) -> None: """A single line of text centred in its own invisible axes.""" ax = _text_axes(fig, rect) ax.text(0.5, 0.5, text, ha="center", va="center", transform=ax.transAxes, **kwargs) @@ -329,12 +346,11 @@ def _kv_rows(ax, pairs: Sequence[tuple[str, str]], *, xs: tuple[float, float, fl ax.text(x, ypos, text, fontsize=fontsize, ha=ha, va="top", transform=ax.transAxes) -def _save_page(pdf: PdfPages, fig: plt.Figure) -> None: +def _save_page(pdf: PdfPages, fig: Figure) -> None: pdf.savefig(fig) - plt.close(fig) -def _add_page_header(fig: plt.Figure, *, page_no: int, entry: str) -> None: +def _add_page_header(fig: Figure, *, page_no: int, entry: str) -> None: """RCSB-style running header (the cover page has none).""" ax = _text_axes(fig, (0.07, 0.952, 0.86, 0.036)) for x, text, ha in ((0.0, f"Page {page_no}", "left"), (0.5, _REPORT_TITLE, "center"), (1.0, entry, "right")): @@ -342,7 +358,7 @@ def _add_page_header(fig: plt.Figure, *, page_no: int, entry: str) -> None: ax.plot([0.0, 1.0], [0.18, 0.18], color=_HEADER_RULE, linewidth=0.6, transform=ax.transAxes) -def _draw_info_box(fig: plt.Figure, *, x: float, y: float, w: float, h: float, lines: Sequence[str]) -> None: +def _draw_info_box(fig: Figure, *, x: float, y: float, w: float, h: float, lines: Sequence[str]) -> None: """Square-corner pink/red cover callout, closer to wwPDB style.""" ax = _text_axes(fig, (x, y, w, h)) ax.add_patch( @@ -358,7 +374,7 @@ def _draw_info_box(fig: plt.Figure, *, x: float, y: float, w: float, h: float, l color="#111111", transform=ax.transAxes) -def _draw_meta_block(fig: plt.Figure, *, x: float, y: float, w: float, h: float, pairs: Sequence[tuple[str, str]]) -> None: +def _draw_meta_block(fig: Figure, *, x: float, y: float, w: float, h: float, pairs: Sequence[tuple[str, str]]) -> None: """Right-aligned label, colon, then value (RCSB style).""" ax = _text_axes(fig, (x, y, w, h)) if pairs: @@ -366,7 +382,7 @@ def _draw_meta_block(fig: plt.Figure, *, x: float, y: float, w: float, h: float, def _draw_section_heading( - fig: plt.Figure, *, x: float, y: float, w: float, h: float, number: str, title: str + fig: Figure, *, x: float, y: float, w: float, h: float, number: str, title: str ) -> None: """Large numbered section heading with RCSB-like spacing.""" ax = _text_axes(fig, (x, y, w, h)) @@ -443,7 +459,7 @@ def _metric_rows_for_slider_panel( def _draw_percentile_legend( - fig: plt.Figure, + fig: Figure, *, x: float, y: float, @@ -462,7 +478,7 @@ def _draw_percentile_legend( def _draw_slider_panel( - fig: plt.Figure, + fig: Figure, *, top: float, height: float, @@ -576,7 +592,7 @@ def _draw_slider_panel( # --------------------------------------------------------------------------- def _draw_fixed_table( - fig: plt.Figure, + fig: Figure, *, x: float, y_top: float, @@ -701,9 +717,39 @@ def _quality_page( ) _draw_slider_panel(fig, top=0.775, height=0.56, row=row, include_overall=True) + _draw_score_provenance(fig, row, y=0.055) _save_page(pdf, fig) +def _draw_score_provenance(fig: Figure, row: Mapping[str, Any], *, y: float) -> None: + """Make score availability, recalculation and cutoff assumptions visible.""" + components = meta_score_components(row) + if components: + lines = [f"Meta recomputed from {len(components)}/{len(META_SCORE_FEATURES)} available features. " + f"CSV score: {_format_raw(_safe_float(row.get('interface_meta_score')))}; " + f"report score: {_format_raw(_row_meta_score(row))}."] + elif _row_meta_score(row) is not None: + lines = ["Stored CSV score used: raw calibrated features unavailable; calibration is unverified."] + else: + lines = ["Meta unavailable: no comparable calibrated features."] + backend = infer_backend(row) or "pooled" + lines.append(f"Report calibration: {CALIBRATION_ID} ({backend}); " + f"CSV calibration: {row.get('metascore_calibration') or 'unknown'}.") + lines.append("Contributing features: " + (", ".join(_FEATURE_DISPLAY[f] for f in components) or "none")) + status = calibration_parameter_status(row) + if status == "unknown": + lines.append("Scoring cutoffs unknown in this CSV; calibration applicability is unknown.") + else: + lines.append("Cutoffs (Å): contacts={}; PAE row filter={}; ipSAE PAE={}.".format( + *(_format_raw(_safe_float(row.get(k))) for k in ("contact_thresh", "pae_filter", "ipsae_pae_cutoff")))) + if status == "custom": + lines.append("Custom cutoffs: reference calibration is unchanged; percentiles are exploratory.") + wrapped = [part for line in lines for part in textwrap.wrap(line, width=108)] + ax = _text_axes(fig, (0.10, y, 0.80, 0.12)) + for i, line in enumerate(wrapped): + ax.text(0., 1. - i * .14, line, fontsize=8, va="top", color="#555", transform=ax.transAxes) + + def _raw_cells(row: Mapping[str, Any], columns: Sequence[str]) -> list[str]: return [_format_raw(_safe_float(row.get(c))) for c in columns] @@ -715,21 +761,9 @@ def _per_interface_page( section_no: str, rows: Sequence[Mapping[str, Any]], page_no: int, -) -> None: - fig = _new_figure() - _add_page_header(fig, page_no=page_no, entry=entry_id) - _draw_section_heading(fig, x=0.07, y=0.91, w=0.86, h=0.03, - number=section_no, title="Per-interface raw scores") - - intro_ax = _text_axes(fig, (0.10, 0.83, 0.80, 0.06)) - for y, line in ( - (1.0, "Each row is one chain pair detected by AlphaJudge."), - (0.55, f"The Meta column is the averaged percentile across the {len(META_SCORE_FEATURES)} " - "metascore features (higher is better)."), - ): - intro_ax.text(0.0, y, line, fontsize=9, ha="left", va="top", transform=intro_ax.transAxes) - - headers = ["Model", "Interface", "Residues", "Meta", "LIS", "ipSAE", "pDockQ2", "ipTM", "PAE", "Sc"] +) -> int: + """Write as many table pages as needed and return the final page number.""" + headers = ["Model", "Interface", "Residues", "Meta", "Nfeat", "LIS", "ipSAE", "pDockQ2", "ipTM", "PAE", "Sc"] raw_columns = ("interface_LIS", "interface_ipSAE", "interface_pDockQ2", "iptm", "average_interface_pae", "interface_sc") body = [ @@ -738,16 +772,31 @@ def _per_interface_page( str(r.get("interface") or ""), str(r.get("interface_num_intf_residues") or ""), _format_raw(_row_meta_score(r)), + str(len(meta_score_components(r))), *_raw_cells(r, raw_columns), ] for r in sorted(rows, key=_meta_sort_key, reverse=True) ] - _draw_fixed_table( - fig, x=0.07, y_top=0.78, w=0.86, headers=headers, rows=body, - col_fracs=[0.18, 0.10, 0.10, 0.08, 0.08, 0.09, 0.10, 0.07, 0.08, 0.12], - row_height=0.024, - ) - _save_page(pdf, fig) + for start in range(0, len(body), _INTERFACE_ROWS_PER_PAGE): + fig = _new_figure() + current_page = page_no + start // _INTERFACE_ROWS_PER_PAGE + _add_page_header(fig, page_no=current_page, entry=entry_id) + title = "Per-interface raw scores" + (" (continued)" if start else "") + _draw_section_heading(fig, x=0.07, y=0.91, w=0.86, h=0.03, + number=section_no, title=title) + intro_ax = _text_axes(fig, (0.10, 0.83, 0.80, 0.06)) + for y, line in ( + (1.0, "Each row is one chain pair detected by AlphaJudge."), + (0.55, f"Meta averages up to {len(META_SCORE_FEATURES)} available feature percentiles; " + "Nfeat gives the count used for each row."), + ): + intro_ax.text(0.0, y, line, fontsize=9, ha="left", va="top", transform=intro_ax.transAxes) + _draw_fixed_table( + fig, x=0.07, y_top=0.78, w=0.86, headers=headers, rows=body[start:start + _INTERFACE_ROWS_PER_PAGE], + col_fracs=[.17, .09, .09, .08, .06, .08, .09, .10, .08, .08, .08], row_height=.024, + ) + _save_page(pdf, fig) + return page_no + max(0, (len(body) - 1) // _INTERFACE_ROWS_PER_PAGE) def _format_residue_tick(value: float, _pos: int | None = None) -> str: @@ -777,6 +826,7 @@ def _pae_vmax(matrix: np.ndarray, max_error: float | None) -> float: return float(math.ceil(observed / 5.0) * 5.0) +@_report_style() def render_pae_png( out_path: str | Path, pae_matrix: Any, @@ -793,8 +843,6 @@ def render_pae_png( look like the in-report graphic) and indirectly by reports that embed the resulting PNG. """ - _setup_rcparams() - try: matrix = pae_matrix if isinstance(pae_matrix, np.ndarray) else np.asarray(pae_matrix, dtype=float) except Exception as e: @@ -813,7 +861,9 @@ def render_pae_png( out_path = Path(out_path) out_path.parent.mkdir(parents=True, exist_ok=True) - fig, ax = plt.subplots(figsize=figsize) + fig = Figure(figsize=figsize) + FigureCanvasAgg(fig) + ax = fig.subplots() im = ax.imshow( matrix, cmap=_PAE_CMAP, @@ -861,7 +911,6 @@ def render_pae_png( fig.tight_layout() fig.savefig(str(out_path), dpi=dpi, bbox_inches="tight") - plt.close(fig) return out_path @@ -947,7 +996,7 @@ def _aggregate_cover_page( scores: Sequence[float], top_rows: Sequence[tuple[str, float, Mapping[str, Any]]], backends: Mapping[str, int], -) -> None: +) -> int: fig = _new_figure() _label_axes(fig, (0.07, 0.87, 0.86, 0.06), _REPORT_TITLE, fontsize=22, fontweight="bold") @@ -1007,12 +1056,24 @@ def _aggregate_cover_page( [str(i), _truncate(name, 34), _format_raw(score), *_raw_cells(row, raw_columns)] for i, (name, score, row) in enumerate(top_rows, start=1) ] + col_fracs = [0.07, 0.34, 0.09, 0.09, 0.10, 0.09, 0.10, 0.12] _draw_fixed_table( - fig, x=0.07, y_top=0.285, w=0.86, headers=headers, rows=body, - col_fracs=[0.07, 0.34, 0.09, 0.09, 0.10, 0.09, 0.10, 0.12], + fig, x=0.07, y_top=0.285, w=0.86, headers=headers, rows=body[:_AGGREGATE_COVER_ROWS], + col_fracs=col_fracs, row_height=0.020, ) _save_page(pdf, fig) + page_no = 1 + for start in range(_AGGREGATE_COVER_ROWS, len(body), _AGGREGATE_ROWS_PER_PAGE): + page_no += 1 + fig = _new_figure() + _add_page_header(fig, page_no=page_no, entry="Aggregate ranking") + _label_axes(fig, (0.07, 0.90, 0.86, 0.03), + f"Top {len(top_rows)} interfaces by meta score (continued)", fontsize=13, fontweight="bold") + _draw_fixed_table(fig, x=.07, y_top=.85, w=.86, headers=headers, rows=body[start:start + _AGGREGATE_ROWS_PER_PAGE], + col_fracs=col_fracs, row_height=.024) + _save_page(pdf, fig) + return page_no def _interface_summary_page( @@ -1047,13 +1108,7 @@ def _interface_summary_page( _draw_slider_panel(fig, top=0.79, height=0.62, row=row, include_overall=True) - note_ax = _text_axes(fig, (0.10, 0.07, 0.80, 0.06)) - note_ax.text( - 0.5, 1.0, - "Black marker shows this interface's percentile rank against the AlphaJudge " - "interacting (positive) benchmark pairs (higher = better).", - ha="center", va="top", fontsize=9, color="#555", transform=note_ax.transAxes, - ) + _draw_score_provenance(fig, row, y=0.025) _save_page(pdf, fig) @@ -1076,6 +1131,7 @@ def _find_pae_png(run_dir: Path, model_used: str) -> Path | None: return next((c for c in candidates if c.is_file()), None) +@_report_style() def generate_per_run_report( run_dir: str | Path, *, @@ -1084,8 +1140,6 @@ def generate_per_run_report( ) -> Path | None: """Build a per-run report.pdf next to ``interfaces.csv``.""" - _setup_rcparams() - run_dir = Path(run_dir) interfaces_csv = run_dir / csv_name if not interfaces_csv.exists(): @@ -1151,7 +1205,7 @@ def generate_per_run_report( quality_section_no = 1 if show_interface_table: page_no += 1 - _per_interface_page(pdf, entry_id=entry_id, section_no="1", rows=rows, page_no=page_no) + page_no = _per_interface_page(pdf, entry_id=entry_id, section_no="1", rows=rows, page_no=page_no) quality_section_no = 2 for i, row in enumerate(interface_rows): @@ -1209,6 +1263,7 @@ def generate_per_run_report( return out_pdf +@_report_style() def generate_aggregate_report( summary_csv: str | Path, *, @@ -1222,8 +1277,6 @@ def generate_aggregate_report( in the merged CSV). A multimer with 15 interfaces contributes 15 points. """ - _setup_rcparams() - summary_csv = Path(summary_csv) if not summary_csv.exists(): logger.warning("summary CSV not found: %s", summary_csv) @@ -1283,7 +1336,7 @@ def generate_aggregate_report( out_pdf = Path(out_pdf) out_pdf.parent.mkdir(parents=True, exist_ok=True) with PdfPages(str(out_pdf)) as pdf: - _aggregate_cover_page( + ranking_pages = _aggregate_cover_page( pdf, summary_csv=summary_csv, n_complexes=len(best_per_complex), @@ -1299,7 +1352,7 @@ def generate_aggregate_report( interface_label=iface or "?", row=r, cohort_position=(rank, len(ranked_per_page)), - page_no=1 + rank, + page_no=ranking_pages + rank, ) for ev_rank, (cname, crow) in enumerate(complex_evidence, start=1): @@ -1313,7 +1366,7 @@ def generate_aggregate_report( row=crow, pae_path=pae_path, model_label=model_label, - page_no=1 + len(ranked_per_page) + ev_rank, + page_no=ranking_pages + len(ranked_per_page) + ev_rank, complex_label=cname, ) diff --git a/src/alphajudge/runner.py b/src/alphajudge/runner.py index 9be8331..234e33f 100644 --- a/src/alphajudge/runner.py +++ b/src/alphajudge/runner.py @@ -6,20 +6,27 @@ import os from concurrent.futures import ProcessPoolExecutor, as_completed from functools import partial +from dataclasses import dataclass -from .parsers import pick_parser +from . import cache +from .parsers import BaseParser, ParseOptions, Run, pick_parser from .complex import Complex +from .geometry import chain_boundaries from .confidence import IPTM_SCOPE_CHAIN_PAIR, IPTM_SCOPE_GLOBAL -from .meta_score import interface_meta_score +from .meta_score import ( + CALIBRATION_ID, calibration_parameter_status, interface_meta_score, meta_score_components, +) from .report import generate_per_run_report, render_pae_png logger = logging.getLogger(__name__) -# Cached per-run CSVs from older releases are safe to reuse only when they -# contain every field introduced by the current output contract. +# Schema validation supplements the provenance manifest. _REQUIRED_CACHE_COLUMNS = frozenset({ "interface_ccc", "interface_expected_contacts", "global_confidence_scope", "iptm_scope", + "backend", "structure_file", + "contact_thresh", "pae_filter", "ipsae_pae_cutoff", "metascore_calibration", + "metascore_calibration_status", "metascore_features", "metascore_feature_count", }) @@ -27,8 +34,64 @@ def _float_or_nan(value) -> float: return float(value) if value is not None else float("nan") +def _write_pae_png(path: Path, confidence, chain_boundaries: list[float]) -> dict | None: + """Rendering is optional output and must never invalidate scored rows.""" + try: + result = render_pae_png(path, confidence.pae_matrix, + chain_boundaries=chain_boundaries, figsize=(8, 8)) + if result is None: + logger.warning("PAE renderer produced no output for %s", path) + return None + fingerprint = cache.file_fingerprint(path) + logger.info("wrote %s", path) + return fingerprint + except Exception: + logger.exception("Could not create PAE heatmap %s", path) + return None + + +def _png_is_current(path: Path, fingerprint: dict | None) -> bool: + try: + return fingerprint == cache.file_fingerprint(path) + except OSError: + return False + + +def _render_model_pae(directory: Path, run: Run, model: str) -> dict | None: + try: + structure, confidence = run.load_model(model) + _, _, indices = BaseParser._maps(structure) + return _write_pae_png(directory / f"pae_{model}.png", confidence, chain_boundaries(indices)) + except Exception: + logger.exception("Could not load PAE for %s in %s", model, directory) + return None + + +def _ensure_pae_pngs(directory: Path, csv_path: Path, request: dict, *, af2_structure: str) -> None: + """Repair optional plots; every failure leaves the scoring cache usable.""" + try: + manifest = cache.read_manifest(csv_path) + if manifest.get("request") != request: + return # A concurrent scoring run replaced the manifest. + pngs = manifest.get("pae_pngs", {}).copy() + missing = [model for model in manifest["models"] + if not _png_is_current(directory / f"pae_{model}.png", pngs.get(model))] + if not missing: + return + run = pick_parser(directory).parse_run(directory, options=ParseOptions(af2_structure)) + for model in missing: + pngs.pop(model, None) # A failed repair must remain eligible for retry. + fingerprint = _render_model_pae(directory, run, model) + if fingerprint is not None: + pngs[model] = fingerprint + cache.update_pae_pngs(csv_path, request, pngs) + except Exception: + logger.exception("Could not refresh PAE heatmaps in %s; cached scores remain usable", directory) + + def _interface_row( - job: str, model: str, iface, confidence, global_score: float, skip_biophysical_scores: bool + job: str, model: str, iface, confidence, global_score: float, skip_biophysical_scores: bool, + *, backend: str = "", structure_file: str = "", ) -> dict: """One output CSV row for a scored interface.""" pd2, _ = iface.pDockQ2() @@ -36,6 +99,11 @@ def _interface_row( row = { "jobs": job, "model_used": model, + "backend": backend, + "structure_file": structure_file, + "contact_thresh": iface.c.contact_thresh, + "pae_filter": iface.c.pae_filter, + "ipsae_pae_cutoff": iface.c.ipsae_pae_cutoff, "interface": iface.label, "global_confidence_scope": confidence.global_confidence_scope, "iptm_scope": IPTM_SCOPE_CHAIN_PAIR if pair_iptm is not None else IPTM_SCOPE_GLOBAL, @@ -73,23 +141,27 @@ def _interface_row( "interface_solv_en": iface.int_solv_en, }) row["interface_meta_score"] = interface_meta_score(row) + components = meta_score_components(row) + row["metascore_calibration"] = CALIBRATION_ID + row["metascore_calibration_status"] = calibration_parameter_status(row) + row["metascore_features"] = ";".join(components) + row["metascore_feature_count"] = len(components) return row -def process( - directory: str, - contact_thresh: float, - pae_filter: float, - models_to_analyse: str, - ipsae_pae_cutoff: float = 10.0, - *, - per_run_csv_name: str = "interfaces.csv", - skip_pae_png: bool = False, - skip_biophysical_scores: bool = False, -) -> Path: - d = Path(directory) - parser = pick_parser(d) - run = parser.parse_run(d) +@dataclass +class _ScoredRun: + rows: list[dict] + models: list[str] + complete: bool + pae_pngs: dict + + +def _score_run(d: Path, run: Run, settings: cache.ScoringSettings, *, skip_pae_png: bool) -> _ScoredRun: + contact_thresh, pae_filter = settings.contact_thresh, settings.pae_filter + models_to_analyse = settings.models_to_analyse + ipsae_pae_cutoff = settings.ipsae_pae_cutoff + skip_biophysical_scores = settings.skip_biophysical_scores models = run.order[:1] if models_to_analyse == "best" else run.order job = d.resolve().name @@ -98,11 +170,12 @@ def process( models_processed = 0 total_interfaces = 0 dropped_by_pae = 0 + complete = True + pae_pngs = {} for m in models: try: structure, confidence = run.load_model(m) comp = Complex(structure, confidence, contact_thresh, pae_filter, ipsae_pae_cutoff) - models_processed += 1 total_interfaces += len(comp.interfaces) global_score = ( @@ -111,35 +184,30 @@ def process( else (comp.interfaces[0].pDockQ if comp.interfaces else float("nan")) ) + model_rows = [] for iface in comp.interfaces: if iface.num_intf_residues == 0: continue if iface.average_interface_pae > pae_filter: dropped_by_pae += 1 continue - rows.append( - _interface_row(job, m, iface, confidence, global_score, skip_biophysical_scores) + model_rows.append( + _interface_row(job, m, iface, confidence, global_score, skip_biophysical_scores, + backend=run.source, structure_file=run.structure_files.get(m, "")) ) if not skip_pae_png: - pae_png = d / f"pae_{m}.png" - try: - render_pae_png( - pae_png, - confidence.pae_matrix, - chain_boundaries=comp.chain_boundaries, - figsize=(8, 8), - ) - logger.info(f"wrote {pae_png}") - except Exception as e: - logger.error(f"Could not create PAE heatmap {pae_png}: {e}") - - logger.info(f"processed model: {m} via {parser.name}") - except Exception as e: - logger.error(f"error processing model {m}: {e}") + fingerprint = _write_pae_png(d / f"pae_{m}.png", confidence, comp.chain_boundaries) + if fingerprint is not None: + pae_pngs[m] = fingerprint + + rows.extend(model_rows) + models_processed += 1 + logger.info("processed model: %s via %s", m, run.source) + except Exception: + complete = False + logger.exception("error processing model %s in %s", m, d) - out = d / per_run_csv_name - out.parent.mkdir(parents=True, exist_ok=True) # per_run_csv_name may hold a subdir if not rows: # Explain *why* the CSV is empty instead of writing a silent zero-byte file # (see https://github.com/KosinskiLab/AlphaJudge/issues/17). The common case @@ -160,14 +228,59 @@ def process( ) else: reason = "all detected interfaces had zero interface residues" - logger.warning(f"no interface rows for {job}: {reason}; writing empty {out}") + logger.warning("no interface rows for %s: %s", job, reason) + + return _ScoredRun(rows, list(models), complete, pae_pngs) - with out.open("w", newline="") as f: - if rows: - w = csv.DictWriter(f, fieldnames=list(rows[0])) + +def _publish_scores(d: Path, out: Path, run: Run, result: _ScoredRun, + request: dict, settings: cache.ScoringSettings, cache_validation: str) -> None: + """Publish rows and provenance without certifying another writer's CSV.""" + with cache.atomic_text(out) as f: + if result.rows: + w = csv.DictWriter(f, fieldnames=list(result.rows[0])) w.writeheader() - w.writerows(rows) - logger.info(f"wrote {out}") + w.writerows(result.rows) + f.flush() + # Hash our own temporary file, before publishing it. A concurrent writer + # must not get its CSV certified with this run's settings. + csv_sha256 = cache.file_digest(Path(f.name)) + # Input files may be written by a predictor while scoring is in progress. + # Such a run must be recomputed on the next invocation. + unchanged = request == cache.request_identity(d, out, cache_validation=cache_validation, settings=settings) + cache.write_manifest(out, request, csv_sha256=csv_sha256, complete=result.complete and unchanged, + backend=run.source, models=result.models, + structure_files={m: run.structure_files[m] for m in result.models if m in run.structure_files}, + pae_pngs=result.pae_pngs) + logger.info("wrote %s", out) + + +def process( + directory: str, + contact_thresh: float, + pae_filter: float, + models_to_analyse: str, + ipsae_pae_cutoff: float = 10.0, + *, + per_run_csv_name: str = "interfaces.csv", + skip_pae_png: bool = False, + skip_biophysical_scores: bool = False, + cache_validation: str = "stat", + af2_structure: str = "relaxed", +) -> Path: + d = Path(directory) + out = d / per_run_csv_name + cache.manifest_path(out).unlink(missing_ok=True) + parser = pick_parser(d) + options = ParseOptions(af2_structure) + settings = cache.ScoringSettings(contact_thresh=contact_thresh, pae_filter=pae_filter, + models_to_analyse=models_to_analyse, ipsae_pae_cutoff=ipsae_pae_cutoff, + skip_biophysical_scores=skip_biophysical_scores, + parser_options=parser.cache_options(options)) + request = cache.request_identity(d, out, cache_validation=cache_validation, settings=settings) + run = parser.parse_run(d, options=options) + result = _score_run(d, run, settings, skip_pae_png=skip_pae_png) + _publish_scores(d, out, run, result, request, settings, cache_validation) return out @@ -194,18 +307,20 @@ def _read_csv_rows(path: Path) -> list[dict]: return list(csv.DictReader(f)) -def _read_reusable_csv(path: Path) -> list[dict] | None: - """Return cached rows only when the file satisfies the current schema.""" +def _read_reusable_csv(path: Path, request: dict) -> list[dict] | None: + """Return cached rows only when schema, inputs, options and software match.""" + if not cache.matches(path, request): + logger.info("existing %s has absent or stale provenance; recomputing", path) + return None rows = _read_csv_rows(path) if not rows: - logger.info(f"existing {path} is empty; recomputing") + logger.info("existing %s is empty; recomputing", path) return None missing = sorted(_REQUIRED_CACHE_COLUMNS - set(rows[0])) if missing: logger.info( - f"existing {path} is missing required column(s) " - f"{', '.join(missing)}; recomputing" + "existing %s is missing required column(s) %s; recomputing", path, ", ".join(missing) ) return None return rows @@ -223,6 +338,8 @@ def _process_one_run( skip_pae_png: bool, skip_biophysical_scores: bool, write_per_run_report: bool = False, + cache_validation: str = "stat", + af2_structure: str = "relaxed", ) -> tuple[str, list[dict]]: """ Worker: process a single run dir (or reuse its per-run CSV) and return @@ -231,18 +348,23 @@ def _process_one_run( d = Path(d_str) csv_path = d / per_run_csv_name - # Reuse a precomputed CSV only if its header satisfies the current output - # contract. This prevents an upgrade from silently aggregating stale rows. rows = None if csv_path.exists() and not force_recompute: try: - rows = _read_reusable_csv(csv_path) + settings = cache.ScoringSettings( + contact_thresh=contact_thresh, pae_filter=pae_filter, + models_to_analyse=models_to_analyse, ipsae_pae_cutoff=ipsae_pae_cutoff, + skip_biophysical_scores=skip_biophysical_scores, + parser_options=pick_parser(d).cache_options(ParseOptions(af2_structure)), + ) + request = cache.request_identity(d, csv_path, settings=settings, cache_validation=cache_validation) + rows = _read_reusable_csv(csv_path, request) except Exception as e: - logger.warning(f"could not reuse {csv_path}; recomputing: {e}") + logger.warning("could not reuse %s; recomputing: %s", csv_path, e) if rows is not None and summary_csv is not None: - logger.info(f"reused existing {csv_path} for aggregation") + logger.info("reused existing %s for aggregation", csv_path) elif rows is not None: - logger.info(f"reused existing {csv_path}; skipping recompute") + logger.info("reused existing %s; skipping recompute", csv_path) if rows is None: csv_path = process( @@ -254,13 +376,17 @@ def _process_one_run( per_run_csv_name=per_run_csv_name, skip_pae_png=skip_pae_png, skip_biophysical_scores=skip_biophysical_scores, + cache_validation=cache_validation, + af2_structure=af2_structure, ) + elif not skip_pae_png: + _ensure_pae_pngs(d, csv_path, request, af2_structure=af2_structure) if write_per_run_report: try: generate_per_run_report(d, csv_name=per_run_csv_name) except Exception as e: # pragma: no cover - defensive - logger.warning(f"per-run report failed for {d}: {e}") + logger.warning("per-run report failed for %s: %s", d, e) if summary_csv is None: return d_str, [] @@ -268,7 +394,7 @@ def _process_one_run( try: rows = _read_csv_rows(csv_path) except Exception as e: - logger.error(f"failed reading {csv_path} for aggregation: {e}") + logger.error("failed reading %s for aggregation: %s", csv_path, e) return d_str, [] # Absolute source_dir lets the aggregate report locate per-run side files # (PAE PNGs, etc.) from the summary CSV. @@ -300,6 +426,8 @@ def process_many( skip_pae_png: bool = False, skip_biophysical_scores: bool = False, write_per_run_report: bool = False, + cache_validation: str = "stat", + af2_structure: str = "relaxed", ) -> Path | None: """ Process one or more directories. Optionally recurse into nested directories @@ -320,14 +448,14 @@ def process_many( for p in paths: rp = Path(p).resolve() if not rp.exists(): - logger.warning(f"path does not exist: {rp}") + logger.warning("path does not exist: %s", rp) elif recursive and rp.is_dir(): run_dirs.update(dict.fromkeys(_discover_run_dirs(rp))) elif _is_run_dir(rp): run_dirs[rp] = None else: logger.warning( - f"no supported run detected at {rp} (use --recursive to search within)" + "no supported run detected at %s (use --recursive to search within)", rp ) if not run_dirs: @@ -350,23 +478,25 @@ def process_many( skip_pae_png=skip_pae_png, skip_biophysical_scores=skip_biophysical_scores, write_per_run_report=write_per_run_report, + cache_validation=cache_validation, + af2_structure=af2_structure, ) aggregated_rows: list[dict] = [] - logger.info(f"Processing {len(run_dirs)} runs with {cores} cores") + logger.info("Processing %s runs with %s cores", len(run_dirs), cores) if cores == 1: for d in run_dirs: try: aggregated_rows.extend(worker(str(d))[1]) - except Exception as e: - logger.error(f"failed processing {d}: {e}") + except Exception: + logger.exception("failed processing %s", d) else: # Important: logging from multiple processes can interleave; acceptable. with ProcessPoolExecutor(max_workers=cores) as ex: for fut in as_completed([ex.submit(worker, str(d)) for d in run_dirs]): try: aggregated_rows.extend(fut.result()[1]) - except Exception as e: - logger.error(f"worker failed: {e}") + except Exception: + logger.exception("worker failed") if not summary_csv: return None @@ -386,6 +516,6 @@ def process_many( w.writerows(aggregated_rows) logger.info( - f"wrote summary {summary_path} ({len(aggregated_rows)} rows from {len(run_dirs)} runs)" + "wrote summary %s (%s rows from %s runs)", summary_path, len(aggregated_rows), len(run_dirs) ) return summary_path diff --git a/test/structure_helpers.py b/test/structure_helpers.py new file mode 100644 index 0000000..dc0a8f6 --- /dev/null +++ b/test/structure_helpers.py @@ -0,0 +1,26 @@ +"""Small coordinate fixtures shared by scoring regressions.""" +from collections.abc import Mapping + +import numpy as np +from Bio.PDB import Atom, Chain, Model, Residue, Structure + + +def make_structure(*, chains="AB", residues=2, separation=5., plddt=90.): + structure = Structure.Structure("test") + model = Model.Model(0) + structure.add(model) + serial = 1 + for chain_no, chain_id in enumerate(chains): + chain = Chain.Chain(chain_id) + model.add(chain) + bfactor = plddt[chain_id] if isinstance(plddt, Mapping) else plddt + for res_no in range(1, residues + 1): + residue = Residue.Residue((" ", res_no, " "), "ALA", " ") + chain.add(residue) + for name, dx in (("CA", 0.), ("CB", .5)): + residue.add(Atom.Atom( + name, np.array([res_no * 3.8 + dx, chain_no * separation, 0.]), + bfactor, 1., " ", name, serial, element="C", + )) + serial += 1 + return structure diff --git a/test/test_af3x.py b/test/test_af3x.py index 0baf76e..4f64b40 100644 --- a/test/test_af3x.py +++ b/test/test_af3x.py @@ -98,11 +98,32 @@ def test_af3x_summary_minima_cannot_replace_missing_full_pae(tmp_path, caplog): summary = json.loads((model_dir / "summary_confidences.json").read_text()) summary["chain_pair_pae_min"] = [[1.] * 3] * 3 (model_dir / "summary_confidences.json").write_text(json.dumps(summary)) - with pytest.raises(ValueError, match="residue-level PAE"): + with pytest.raises(FileNotFoundError, match="confidences.json"): AF3Parser().parse_run(run_dir).load_model("seed-1_sample-0") output = process(str(run_dir), 8., 100., "best", skip_pae_png=True, skip_biophysical_scores=True) assert output.read_text() == "" - assert "residue-level PAE" in caplog.text + assert "Required JSON file not found" in caplog.text + + +def test_af3_combined_summary_with_full_pae_remains_supported(tmp_path): + run_dir, _, _ = _write_run(tmp_path) + model_dir = run_dir / "seed-1_sample-0" + matrix = json.loads((model_dir / "confidences.json").read_text()) + summary_path = model_dir / "summary_confidences.json" + summary = json.loads(summary_path.read_text()) + summary.update(matrix) + summary_path.write_text(json.dumps(summary)) + (model_dir / "confidences.json").unlink() + _, confidence = AF3Parser().parse_run(run_dir).load_model("seed-1_sample-0") + assert confidence.pae_matrix.shape == (4, 4) + + +def test_af3_present_minima_only_file_has_distinct_error(tmp_path): + run_dir, _, _ = _write_run(tmp_path) + path = run_dir / "seed-1_sample-0/confidences.json" + path.write_text(json.dumps({"chain_pair_pae_min": [[1., 1.], [1., 1.]]})) + with pytest.raises(ValueError, match="confidences.json.*only chain-pair minima"): + AF3Parser().parse_run(run_dir).load_model("seed-1_sample-0") @pytest.mark.parametrize("field", ["pae", "predicted_aligned_error"]) @@ -261,27 +282,26 @@ def test_plain_af3_global_confidence_still_contributes_to_metascore(): def test_cached_scores_without_scope_metadata_are_recomputed(tmp_path, monkeypatch): from alphajudge import runner - stale = tmp_path / "interfaces.csv" + run_dir, _, _ = _write_run(tmp_path) + stale = run_dir / "interfaces.csv" stale.write_text("jobs,interface_ccc,interface_expected_contacts\nold,1,2\n") calls = [] + original_process = runner.process def recompute(directory, *args, **kwargs): calls.append(directory) - stale.write_text( - "jobs,interface_ccc,interface_expected_contacts,global_confidence_scope,iptm_scope\n" - "new,1,2,includes_excluded_tokens,chain_pair\n" - ) - return stale + return original_process(directory, *args, **kwargs) monkeypatch.setattr(runner, "process", recompute) for _ in range(2): _, rows = runner._process_one_run( - str(tmp_path), contact_thresh=8., pae_filter=100., models_to_analyse="best", + str(run_dir), contact_thresh=8., pae_filter=100., models_to_analyse="best", summary_csv="summary.csv", ipsae_pae_cutoff=10., force_recompute=False, per_run_csv_name="interfaces.csv", skip_pae_png=True, skip_biophysical_scores=True, ) - assert rows[0]["jobs"] == "new" - assert calls == [str(tmp_path)] + assert rows[0]["jobs"] == "mixed" + assert rows[0]["global_confidence_scope"] == "includes_excluded_tokens" + assert calls == [str(run_dir)] @pytest.mark.parametrize("ligand_first", [False, True]) diff --git a/test/test_cli.py b/test/test_cli.py new file mode 100644 index 0000000..fd913b2 --- /dev/null +++ b/test/test_cli.py @@ -0,0 +1,26 @@ +import sys + +import pytest + +from alphajudge import cli + + +def test_aggregate_requires_summary_before_scoring(monkeypatch, capsys): + calls = [] + monkeypatch.setattr(cli, "process_many", lambda *a, **kw: calls.append(kw)) + monkeypatch.setattr(sys, "argv", ["alphajudge", "run", "--aggregate_report", "out.pdf"]) + with pytest.raises(SystemExit) as error: + cli.main() + assert error.value.code == 2 + assert "--aggregate_report requires --summary" in capsys.readouterr().err + assert calls == [] + + +def test_cache_and_structure_options_reach_scoring(monkeypatch): + calls = [] + monkeypatch.setattr(cli, "process_many", lambda *a, **kw: calls.append(kw)) + monkeypatch.setattr(sys, "argv", ["alphajudge", "run", "--no-report", + "--cache_validation", "content", "--af2_structure", "unrelaxed"]) + cli.main() + assert calls[0]["cache_validation"] == "content" + assert calls[0]["af2_structure"] == "unrelaxed" diff --git a/test/test_meta_score.py b/test/test_meta_score.py index e47809f..450e42f 100644 --- a/test/test_meta_score.py +++ b/test/test_meta_score.py @@ -2,6 +2,8 @@ import math +import pytest + from alphajudge.meta_score import ( META_SCORE_FEATURES, calibrated_feature_percentile, @@ -56,3 +58,32 @@ def test_inverted_features_score_in_the_expected_direction() -> None: def test_feature_percentiles_clamp_to_unit_interval() -> None: assert calibrated_feature_percentile("interface_LIS", -1.0) == 0.0 assert calibrated_feature_percentile("interface_LIS", 10.0) == 1.0 + + +def test_metascore_components_exclude_missing_and_incompatible_features(): + from alphajudge.meta_score import meta_score_components + + row = {**_complete_row(), "global_confidence_scope": "includes_excluded_tokens", + "iptm_scope": "global", "interface_hb": float("nan")} + parts = meta_score_components(row) + assert set(parts) == set(_complete_row()) - {"interface_area", "iptm", "confidence_score"} + assert interface_meta_score(row) == pytest.approx(sum(parts.values()) / len(parts)) + + +@pytest.mark.parametrize("parameters,status", [ + ({}, "unknown"), + ({"contact_thresh": 8., "pae_filter": 100., "ipsae_pae_cutoff": 10.}, "default"), + ({"contact_thresh": "8", "pae_filter": "100", "ipsae_pae_cutoff": "10"}, "default"), + ({"contact_thresh": 8., "pae_filter": 100., "ipsae_pae_cutoff": 15.}, "custom"), + ({"contact_thresh": 8., "pae_filter": 100., "ipsae_pae_cutoff": "nan"}, "unknown"), +]) +def test_calibration_parameter_status(parameters, status): + from alphajudge.meta_score import calibration_parameter_status + + assert calibration_parameter_status(parameters) == status + + +def test_percentile_registry_only_lists_supported_features(): + from alphajudge.meta_score import BENCHMARK_QUANTILES, FEATURE_DIRECTIONS + + assert set(FEATURE_DIRECTIONS) == set(BENCHMARK_QUANTILES) diff --git a/test/test_pdockq2_compatibility.py b/test/test_pdockq2_compatibility.py new file mode 100644 index 0000000..6b43a78 --- /dev/null +++ b/test/test_pdockq2_compatibility.py @@ -0,0 +1,25 @@ +"""Preserve the pairwise pDockQ2 variant used by AlphaJudge and IPSAE.""" +import numpy as np +import pytest +from structure_helpers import make_structure + +from alphajudge.complex import Complex +from alphajudge.confidence import Confidence +from alphajudge.parsers import BaseParser + + +def test_pdockq2_matches_pinned_ipsae_asymmetric_trimer(): + # Reference: DunbrackLab/IPSAE ipsae.py at + # 6174cf9e71cb1bd660cc805856a18c4871a6dec3, pDockQ2 block. + # Expected values were captured by running that script on this structure + # and PAE, with shared pLDDT inputs. They are the maximum of its two + # directional scores for each pair, not a chain-versus-rest score. + structure = make_structure(chains="ABC", residues=3, plddt={"A": 95., "B": 30., "C": 80.}) + pae = np.full((9, 9), 2.) + pae[3:6, :3] = 18. + pae[6:, 3:6] = 12. + chains, residue_map, _ = BaseParser._maps(structure) + confidence = Confidence(pae, 31., .7, .6, .68, .68, BaseParser._plddt(chains, residue_map)) + comp = Complex(structure, confidence, 8., 100., 10.) + scores = {iface.label: iface.pDockQ2()[0] for iface in comp.interfaces} + assert scores == pytest.approx({"A_B": .18333777506612065, "B_C": .11509754028543784}, abs=1e-12) diff --git a/test/test_report.py b/test/test_report.py index 0f651a8..1ba4e88 100644 --- a/test/test_report.py +++ b/test/test_report.py @@ -199,7 +199,7 @@ def test_per_run_report_text_counts_features_and_numbers_appendices( _write_csv(tmp_path / "interfaces.csv", rows) assert generate_per_run_report(tmp_path) is not None - assert any(f"across the {len(META_SCORE_FEATURES)} metascore features" in t for t in drawn_text) + assert any(f"up to {len(META_SCORE_FEATURES)} available feature percentiles" in t for t in drawn_text) assert "A.1" in drawn_text assert "Appendix – model model_2_multimer_v3_pred_0" in drawn_text @@ -218,3 +218,123 @@ def test_aggregate_cover_median_averages_the_middle_pair(tmp_path: Path, drawn_t scores = sorted(rep._row_meta_score(r) for r in rows) median = (scores[1] + scores[2]) / 2 assert f"median = {median:.3f}" in drawn_text + + +@pytest.fixture +def captured_pages(monkeypatch): + """Capture actual PDF page artists and their physical positions before saving.""" + from matplotlib.backends.backend_pdf import PdfPages + from matplotlib.text import Text + + pages = [] + original = PdfPages.savefig + + def record(pdf, figure, **kwargs): + texts = [] + for artist in figure.findobj(Text): + point = artist.get_transform().transform(artist.get_position()) + x, y = figure.transFigure.inverted().transform(point) + texts.append((artist.get_text(), x, y)) + pages.append(texts) + return original(pdf, figure, **kwargs) + + monkeypatch.setattr(PdfPages, "savefig", record) + return pages + + +def test_interface_table_paginates_without_losing_rows(tmp_path, captured_pages): + from alphajudge import report as rep + from matplotlib.backends.backend_pdf import PdfPages + + rows = [{**_BASE_ROW, "model_used": f"R{i:02}"} for i in range(36)] + with PdfPages(tmp_path / "table.pdf") as pdf: + rep._per_interface_page(pdf, entry_id="test", section_no="1", rows=rows, page_no=2) + assert len(captured_pages) == 2 + assert all(any(text == "Model" for text, _, _ in page) for page in captured_pages) + table_rows = [(text, y) for page in captured_pages for text, _, y in page if text.startswith("R") and text[1:].isdigit()] + assert sorted(text for text, _ in table_rows) == [f"R{i:02}" for i in range(36)] + assert all(.05 < y < .85 for _, y in table_rows) + assert [text for page in captured_pages for text, _, _ in page if text.startswith("Page ")] == ["Page 2", "Page 3"] + + +def test_large_aggregate_ranking_paginates(tmp_path, captured_pages): + rows = [{**_BASE_ROW, "jobs": f"JOB_{i:02}"} for i in range(42)] + summary = tmp_path / "summary.csv" + _write_csv(summary, rows) + generate_aggregate_report(summary, out_pdf=tmp_path / "aggregate.pdf", top_n=42, max_complexes=1) + # Cover with ten entries, two continuation tables, one interface, one complex. + assert len(captured_pages) == 5 + table_entries = [(text, y) for page in captured_pages[:3] for text, _, y in page if text.startswith("JOB_")] + assert len(table_entries) == 42 + assert all(.05 < y < .9 for _, y in table_entries) + assert [text for page in captured_pages for text, _, _ in page if text.startswith("Page ")] == [f"Page {i}" for i in range(2, 6)] + + +def test_import_preserves_matplotlib_backend_and_figures(): + import subprocess + import sys + + code = """ +import matplotlib +matplotlib.use('svg') +import matplotlib.pyplot as plt +fig = plt.figure() +import alphajudge.report +assert matplotlib.get_backend().lower() == 'svg' +assert plt.fignum_exists(fig.number) +""" + subprocess.run([sys.executable, "-c", code], check=True, capture_output=True, text=True) + + +@pytest.mark.parametrize("operation", ["png", "per_run", "aggregate", "error"]) +def test_render_restores_matplotlib_settings(tmp_path, operation, monkeypatch): + import matplotlib + import matplotlib.pyplot as plt + from matplotlib.figure import Figure + from alphajudge import report as rep + + _write_csv(tmp_path / "interfaces.csv", [dict(_BASE_ROW)]) + with matplotlib.rc_context({"font.family": ["sans-serif"], "font.size": 17, "pdf.fonttype": 3}): + figure = plt.figure() + before = dict(matplotlib.rcParams) + try: + if operation == "png": + rep.render_pae_png(tmp_path / "pae.png", [[1., 2.], [3., 4.]]) + elif operation == "per_run": + rep.generate_per_run_report(tmp_path) + elif operation == "aggregate": + rep.generate_aggregate_report(tmp_path / "interfaces.csv", out_pdf=tmp_path / "a.pdf") + else: + def fail(*args, **kwargs): + raise OSError("write failed") + monkeypatch.setattr(Figure, "savefig", fail) + with pytest.raises(OSError, match="write failed"): + rep.render_pae_png(tmp_path / "pae.png", [[1.]]) + assert dict(matplotlib.rcParams) == before + assert plt.get_fignums() == [figure.number] + finally: + plt.close(figure) + + +def test_report_explains_recalibration_and_custom_cutoffs(tmp_path, drawn_text): + row = {**_BASE_ROW, "contact_thresh": "12", "pae_filter": "100", "ipsae_pae_cutoff": "10", + "metascore_calibration": "historical-fit", "interface_meta_score": "0.99"} + _write_csv(tmp_path / "interfaces.csv", [row]) + generate_per_run_report(tmp_path) + from alphajudge.meta_score import CALIBRATION_ID + text = "\n".join(drawn_text) + assert CALIBRATION_ID in text + assert "recomputed" in text + assert "10/11" in text + assert "CSV score: 0.99" in text + assert "historical-fit" in text + assert "Custom cutoffs" in text + assert "exploratory" in text + + +def test_legacy_report_marks_unknown_parameters_and_stored_score(tmp_path, drawn_text): + _write_csv(tmp_path / "interfaces.csv", [{"jobs": "old", "interface": "A_B", "interface_meta_score": "0.6"}]) + generate_per_run_report(tmp_path) + text = "\n".join(drawn_text) + assert "Stored CSV score" in text + assert "cutoffs unknown" in text diff --git a/test/test_run_integrity.py b/test/test_run_integrity.py new file mode 100644 index 0000000..0a80d64 --- /dev/null +++ b/test/test_run_integrity.py @@ -0,0 +1,531 @@ +"""Regression tests for input provenance, cache validity and failed-model isolation.""" +from __future__ import annotations + +import gzip +import json +import lzma +import os +import pickle +import shutil +from dataclasses import replace +from pathlib import Path + +import numpy as np +import pytest +from Bio.PDB import MMCIFIO, PDBIO +from structure_helpers import make_structure + +from alphajudge import cache, meta_score, report, runner +from alphajudge.complex import Complex +from alphajudge.parsers import BaseParser, ParseOptions +from alphajudge.parsers.af2 import AF2Parser +from alphajudge.parsers.af3 import AF3Parser +from alphajudge.parsers.boltz import Boltz2Parser + + +def write_structure(path: Path, separation: float = 5., plddt: float = 90., chains: str = "AB") -> None: + structure = make_structure(chains=chains, separation=separation, plddt=plddt) + writer = PDBIO() + writer.set_structure(structure) + writer.save(str(path)) + + +@pytest.fixture +def record_process(monkeypatch): + def start(): + calls = [] + original = runner.process + def record(*args, **kwargs): + calls.append(1) + return original(*args, **kwargs) + monkeypatch.setattr(runner, "process", record) + return calls + return start + + +@pytest.fixture +def af2_run(tmp_path): + models = ["model_1_multimer_v3_pred_0", "model_2_multimer_v3_pred_0"] + ranking = {"order": models, "iptm+ptm": dict.fromkeys(models, .68), + "iptm": dict.fromkeys(models, .7), "ptm": dict.fromkeys(models, .6)} + (tmp_path / "ranking_debug.json").write_text(json.dumps(ranking)) + for index, model in enumerate(models): + write_structure(tmp_path / f"unrelaxed_{model}.pdb", separation=7. if index == 0 else 5.) + (tmp_path / f"pae_{model}.json").write_text(json.dumps([ + {"predicted_aligned_error": np.full((4, 4), 2.).tolist()} + ])) + return tmp_path, models + + +def score(directory, **changes): + options = dict(contact_thresh=8., pae_filter=100., models_to_analyse="best", + summary_csv="summary.csv", ipsae_pae_cutoff=10., force_recompute=False, + per_run_csv_name="interfaces.csv", skip_pae_png=True, + skip_biophysical_scores=True) + options.update(changes) + return runner._process_one_run(str(directory), **options)[1] + + +def test_filtered_residues_retain_chain_identity(af2_run): + directory, models = af2_run + structure, confidence = AF2Parser().parse_run(directory).load_model(models[0]) + comp = Complex(structure, confidence, 8., 100., 10.) + first, second = (chain.child_list[0] for chain in comp._chains) + assert first.id == second.id # Same residue number, different chains. + assert first != second + assert first not in list(comp._chains[1]) + residues = [res for chain in comp._chains for res in chain] + for a in residues: + for b in residues: + assert (a == b) == (a in {b}) + assert len(set(residues)) == 4 + assert first.get_parent().get_parent().get_parent() is not None + # Filtering must not reparent the input structure. + assert structure[0]["A"].child_list[0] is not first + assert structure[0]["A"].get_parent() is structure[0] + + +@pytest.mark.parametrize("contact_thresh,expected_status", [(8., "default"), (9., "custom")]) +def test_csv_records_score_provenance(af2_run, contact_thresh, expected_status): + directory, _ = af2_run + row = score(directory, contact_thresh=contact_thresh)[0] + assert float(row["contact_thresh"]) == contact_thresh + assert float(row["pae_filter"]) == 100. + assert float(row["ipsae_pae_cutoff"]) == 10. + assert row["metascore_calibration"] == meta_score.CALIBRATION_ID + assert row["metascore_calibration_status"] == expected_status + features = row["metascore_features"].split(";") + assert int(row["metascore_feature_count"]) == len(features) == 7 + assert "interface_hb" not in features + assert float(row["interface_meta_score"]) == pytest.approx( + sum(meta_score.calibrated_feature_percentile(f, row[f], "af2") for f in features) / len(features) + ) + + +def test_cache_changed_selection_and_threshold(af2_run): + directory, models = af2_run + assert {r["model_used"] for r in score(directory)} == {models[0]} + assert {r["model_used"] for r in score(directory, contact_thresh=6., models_to_analyse="all")} == {models[1]} + + +def test_parser_options_select_structure_and_only_affect_relevant_caches(af2_run): + directory, models = af2_run + write_structure(directory / f"relaxed_{models[0]}.pdb", plddt=40.) + options = ParseOptions(structure_preference="unrelaxed") + run = AF2Parser().parse_run(directory, options=options) + _, confidence = run.load_model(models[0]) + assert confidence.plddt_residue == [90.] * 4 + assert AF2Parser().cache_options(options) == options + assert AF3Parser().cache_options(options) is None + assert Boltz2Parser().cache_options(options) is None + + +def test_invalid_cache_or_parser_option_is_rejected(af2_run): + directory, _ = af2_run + with pytest.raises(TypeError, match="pae_filtter"): + cache.request_identity(directory, directory / "interfaces.csv", pae_filtter=10.) + with pytest.raises(TypeError, match="pae_filtter"): + cache.ScoringSettings(pae_filtter=10.) + with pytest.raises(ValueError, match="Structure preference"): + ParseOptions(structure_preference="arbitrary") + + +@pytest.mark.parametrize("change", [ + {"pae_filter": 1.}, {"ipsae_pae_cutoff": 1.}, + {"skip_biophysical_scores": False}, +]) +def test_cache_invalidates_scoring_options(af2_run, record_process, change): + directory, _ = af2_run + score(directory) + calls = record_process() + score(directory) + assert calls == [] + score(directory, **change) + assert calls == [1] + + +@pytest.mark.parametrize("cache_validation", ["stat", "content"]) +def test_cache_checks_input_and_csv_contents(af2_run, cache_validation): + directory, models = af2_run + first = score(directory, cache_validation=cache_validation)[0] + structure_path = directory / f"unrelaxed_{models[0]}.pdb" + before = structure_path.stat() + write_structure(structure_path, separation=7., plddt=40.) + # Content mode must detect a same-size edit even with a preserved mtime. + # Stat mode needs an observable metadata change; the shared filesystem can + # quantize both mtime and ctime to whole seconds. + mtime = before.st_mtime_ns if cache_validation == "content" else before.st_mtime_ns + 2_000_000_000 + os.utime(structure_path, ns=(before.st_atime_ns, mtime)) + assert structure_path.stat().st_size == before.st_size + second = score(directory, cache_validation=cache_validation)[0] + assert float(first["interface_average_plddt"]) == 90. + assert float(second["interface_average_plddt"]) == 40. + path = directory / "interfaces.csv" + path.write_text(path.read_text().replace(models[0], "edited_model")) + assert score(directory, cache_validation=cache_validation)[0]["model_used"] == models[0] + + +def test_warm_cache_does_not_read_input_payloads(af2_run, monkeypatch): + directory, _ = af2_run + original = score(directory) + digest = cache.file_digest + def guarded_digest(path): + assert path.parent != directory or path.name == "interfaces.csv", f"Re-read input: {path}" + return digest(path) + monkeypatch.setattr(cache, "file_digest", guarded_digest) + assert score(directory) == original + + +def test_content_validation_hashes_inputs_on_warm_reuse(af2_run, monkeypatch): + directory, _ = af2_run + score(directory, cache_validation="content") + calls = [] + digest = cache.file_digest + def recording_digest(path): + calls.append(path) + return digest(path) + monkeypatch.setattr(cache, "file_digest", recording_digest) + score(directory, cache_validation="content") + assert directory / "ranking_debug.json" in calls + assert any(path.suffix == ".pdb" for path in calls) + + +def test_content_validation_detects_edits_hidden_from_metadata(af2_run, monkeypatch): + directory, models = af2_run + csv_path = directory / "interfaces.csv" + stat_request = cache.request_identity(directory, csv_path) + content_request = cache.request_identity(directory, csv_path, cache_validation="content") + original = cache.file_fingerprint + def frozen_input_metadata(path): + key = str(path.relative_to(directory)) + return stat_request["inputs"].get(key) or original(path) + monkeypatch.setattr(cache, "file_fingerprint", frozen_input_metadata) + write_structure(directory / f"unrelaxed_{models[0]}.pdb", plddt=40.) + assert cache.request_identity(directory, csv_path) == stat_request + assert cache.request_identity(directory, csv_path, cache_validation="content") != content_request + + +def test_changing_validation_mode_revalidates_scores(af2_run, record_process): + directory, _ = af2_run + score(directory) + calls = record_process() + score(directory, cache_validation="content") + score(directory, cache_validation="content") + assert calls == [1] + + +def test_cache_detects_added_removed_and_replaced_inputs(af2_run): + directory, models = af2_run + csv_path = directory / "interfaces.csv" + before = cache.request_identity(directory, csv_path) + extra = directory / "extra.json" + extra.write_text("{}") + assert cache.request_identity(directory, csv_path) != before + extra.unlink() + assert cache.request_identity(directory, csv_path) == before + structure = directory / f"unrelaxed_{models[0]}.pdb" + stat = structure.stat() + replacement = structure.with_suffix(".tmp") + replacement.write_bytes(structure.read_bytes()) + os.utime(replacement, ns=(stat.st_atime_ns, stat.st_mtime_ns)) + replacement.replace(structure) + assert cache.request_identity(directory, csv_path) != before + + +def test_png_toggle_generates_missing_plot_without_rescoring(af2_run, monkeypatch): + directory, models = af2_run + original = score(directory) + csv_path = directory / "interfaces.csv" + csv_stat = csv_path.stat() + def no_rescore(*args, **kwargs): + pytest.fail("PNG request recomputed scores") + monkeypatch.setattr(runner, "process", no_rescore) + monkeypatch.setattr(runner, "Complex", no_rescore) + assert score(directory, skip_pae_png=False) == original + assert (directory / f"pae_{models[0]}.png").is_file() + assert csv_path.stat().st_mtime_ns == csv_stat.st_mtime_ns + # A second request should not even reload the model when the plot is current. + monkeypatch.setattr(AF2Parser, "parse_run", no_rescore) + assert score(directory, skip_pae_png=False) == original + assert score(directory, skip_pae_png=True) == original + + +@pytest.mark.parametrize("failure", ["exception", "no_output"]) +def test_failed_png_is_retried_without_invalidating_scores(af2_run, monkeypatch, failure): + directory, _ = af2_run + attempts = [] + def broken_plot(*args, **kwargs): + attempts.append(1) + if failure == "exception": + raise OSError("PNG destination unavailable") + return None + monkeypatch.setattr(runner, "render_pae_png", broken_plot) + original = score(directory, skip_pae_png=False) + manifest = json.loads(cache.manifest_path(directory / "interfaces.csv").read_text()) + assert manifest["complete"] is True + monkeypatch.setattr(runner, "process", lambda *a, **kw: pytest.fail("Rescored after PNG failure")) + assert score(directory, skip_pae_png=False) == original + assert attempts == [1, 1] + + +def test_skipped_png_is_refreshed_after_input_change(af2_run, monkeypatch): + directory, models = af2_run + score(directory, skip_pae_png=False, cache_validation="content") + (directory / f"pae_{models[0]}.json").write_text(json.dumps([ + {"predicted_aligned_error": np.full((4, 4), 7.).tolist()} + ])) + score(directory, skip_pae_png=True, cache_validation="content") + plotted = [] + original = runner.render_pae_png + def recording_plot(path, matrix, **kwargs): + plotted.append(float(matrix[0, 0])) + return original(path, matrix, **kwargs) + monkeypatch.setattr(runner, "render_pae_png", recording_plot) + monkeypatch.setattr(runner, "process", lambda *a, **kw: pytest.fail("Rescored for stale PNG")) + score(directory, skip_pae_png=False, cache_validation="content") + assert plotted == [7.] + + +def test_explicit_unrelaxed_preference_invalidates_cache(af2_run): + directory, models = af2_run + write_structure(directory / f"relaxed_{models[0]}.pdb", plddt=40.) + relaxed = score(directory)[0] + unrelaxed = score(directory, af2_structure="unrelaxed")[0] + assert float(relaxed["interface_average_plddt"]) == 40. + assert float(unrelaxed["interface_average_plddt"]) == 90. + assert unrelaxed["structure_file"] == f"unrelaxed_{models[0]}.pdb" + assert score(directory, af2_structure="unrelaxed")[0] == unrelaxed + + +def test_unrelaxed_preference_falls_back_to_available_relaxed_file(af2_run): + directory, models = af2_run + model = models[0] + (directory / f"unrelaxed_{model}.pdb").rename(directory / f"relaxed_{model}.pdb") + row = score(directory, af2_structure="unrelaxed")[0] + assert row["structure_file"] == f"relaxed_{model}.pdb" + + +def test_copied_cache_does_not_keep_the_old_job_name(af2_run): + directory, _ = af2_run + original_job = score(directory)[0]["jobs"] + relocated = directory.parent / (directory.name + "_relocated") + shutil.copytree(directory, relocated) + row = score(relocated)[0] + assert row["jobs"] == relocated.name != original_job + + +@pytest.mark.parametrize("invalidate", ["software", "missing_manifest", "corrupt_manifest", "force"]) +def test_cache_recomputes_without_current_provenance(af2_run, monkeypatch, record_process, invalidate): + directory, _ = af2_run + filename = "scores/custom.csv" + score(directory, per_run_csv_name=filename) + manifest = cache.manifest_path(directory / filename) + if invalidate == "software": + monkeypatch.setattr(cache, "software_identity", lambda: {"changed_calibration": True}) + elif invalidate == "missing_manifest": + manifest.unlink() + elif invalidate == "corrupt_manifest": + manifest.write_text("{BROKEN") + calls = record_process() + score(directory, per_run_csv_name=filename, force_recompute=invalidate == "force") + assert calls == [1] + score(directory, per_run_csv_name=filename) + assert calls == [1] + + +def test_concurrent_csv_replacement_cannot_receive_wrong_provenance(af2_run, monkeypatch): + directory, models = af2_run + original = cache.write_manifest + def replace_before_manifest(csv_path, *args, **kwargs): + csv_path.write_text(csv_path.read_text().replace(models[0], "other_writer")) + return original(csv_path, *args, **kwargs) + monkeypatch.setattr(cache, "write_manifest", replace_before_manifest) + score(directory) + path = directory / "interfaces.csv" + manifest = json.loads(cache.manifest_path(path).read_text()) + assert not cache.matches(path, manifest["request"]) + + +@pytest.mark.parametrize("suffix,compress", [("", lambda x: x), (".gz", gzip.compress), (".xz", lzma.compress)]) +def test_native_multimer_ranking_and_pickle_pae(af2_run, suffix, compress): + directory, models = af2_run + model = models[0] + (directory / "ranking_debug.json").write_text(json.dumps({"order": [model], "iptm+ptm": {model: .68}})) + (directory / f"pae_{model}.json").unlink() + (directory / f"result_{model}.pkl{suffix}").write_bytes(compress(pickle.dumps({ + "iptm": .7, "ptm": .6, "predicted_aligned_error": np.full((4, 4), 2.) + }))) + _, confidence = AF2Parser().parse_run(directory).load_model(model) + assert (confidence.iptm, confidence.ptm, confidence.confidence_score) == (.7, .6, .68) + assert confidence.pae_matrix.shape == (4, 4) + assert len(score(directory)) == 1 + + +def test_native_multimer_missing_individual_scores_remain_missing(af2_run): + directory, models = af2_run + model = models[0] + (directory / "ranking_debug.json").write_text(json.dumps({"order": [model], "iptm+ptm": {model: .68}})) + _, confidence = AF2Parser().parse_run(directory).load_model(model) + assert confidence.iptm is None and confidence.ptm is None + assert confidence.confidence_score == .68 + + +def test_boltz_backend_does_not_depend_on_model_name(tmp_path): + rows = [] + for name in ("model_foo_model_0", "myinput_model_0"): + directory = tmp_path / name + directory.mkdir() + write_structure(directory / f"{name}.pdb") + (directory / f"confidence_{name}.json").write_text(json.dumps({"iptm": .7, "ptm": .6, "confidence_score": .68})) + np.savez(directory / f"pae_{name}.npz", pae=np.full((4, 4), 2.)) + row = score(directory)[0] + assert row["backend"] == "boltz2" + assert meta_score.infer_backend(row) is None # pooled ladder for Boltz + assert report._detect_backend([row]) == "Boltz-2" + rows.append(row) + assert rows[0]["interface_meta_score"] == rows[1]["interface_meta_score"] + + +@pytest.mark.parametrize("suffix,compress", [("", lambda x: x), (".gz", gzip.compress), (".xz", lzma.compress)]) +def test_corrupt_json_reports_actual_file(tmp_path, suffix, compress): + path = tmp_path / f"confidences.json{suffix}" + path.write_bytes(compress(b"{BROKEN")) + with pytest.raises(ValueError, match="confidences.json"): + BaseParser._read_json(tmp_path / "confidences.json") + assert BaseParser._read_json(tmp_path / "optional.json") == {} + + +@pytest.mark.parametrize("suffix,compress", [(".gz", gzip.compress), (".xz", lzma.compress)]) +def test_corrupt_compressed_stream_reports_actual_file(tmp_path, suffix, compress): + path = tmp_path / f"confidences.json{suffix}" + data = bytearray(compress(b'{"pae": [[0.0]]}')) + data[-5] ^= 0xff # damage the compression trailer, not the JSON + path.write_bytes(data) + with pytest.raises(ValueError, match="confidences.json"): + BaseParser._read_json(tmp_path / "confidences.json") + + +@pytest.mark.parametrize("payload", [None, "{BROKEN", "[]", "{}"]) +def test_af3_required_confidence_errors_keep_file_context(tmp_path, payload): + directory = tmp_path / "seed-1_sample-0" + directory.mkdir() + write_structure(directory / "model.pdb") + (tmp_path / "ranking_scores.csv").write_text("seed,sample,ranking_score\n1,0,0.8\n") + (directory / "summary_confidences.json").write_text(json.dumps({"iptm": .7, "ptm": .6})) + if payload is not None: + (directory / "confidences.json").write_text(payload) + run = AF3Parser().parse_run(tmp_path) + with pytest.raises((FileNotFoundError, ValueError), match="confidences.json") as caught: + run.load_model(run.order[0]) + assert "only chain-pair minima" not in str(caught.value) + + +def test_af2_corrupt_pae_does_not_silently_fall_back_to_pickle(af2_run): + directory, models = af2_run + model = models[0] + (directory / f"pae_{model}.json").write_text("{BROKEN") + (directory / f"result_{model}.pkl").write_bytes(pickle.dumps({"predicted_aligned_error": np.zeros((4, 4))})) + with pytest.raises(ValueError, match=f"pae_{model}.json"): + AF2Parser().parse_run(directory).load_model(model) + + +@pytest.mark.parametrize("pae", [np.zeros((3, 3)), np.zeros((4, 3)), np.zeros(4), + np.full((4, 4), np.nan), np.full((4, 4), -1.)]) +def test_complex_rejects_invalid_pae_before_scoring(af2_run, pae): + directory, models = af2_run + structure, confidence = AF2Parser().parse_run(directory).load_model(models[0]) + with pytest.raises(ValueError, match="PAE"): + Complex(structure, replace(confidence, pae_matrix=pae), 8., 100.) + + +@pytest.mark.parametrize("pae", [np.zeros((5, 5)), np.zeros((3, 3)), np.full((4, 4), np.nan)]) +def test_boltz_rejects_unmappable_pae(tmp_path, pae): + path = tmp_path / "pae.npz" + np.savez(path, pae=pae) + with pytest.raises(ValueError, match="pae.npz"): + Boltz2Parser._load_pae(path, 4, 4, np.arange(4)) + + +def test_boltz_corrupt_pae_is_distinct_from_missing_optional_pae(tmp_path): + path = tmp_path / "pae.npz" + path.write_bytes(b"BROKEN") + with pytest.raises(ValueError, match="pae.npz"): + Boltz2Parser._load_pae(path, 4, 4, np.arange(4)) + matrix, _ = Boltz2Parser._load_pae(None, 4, 4, np.arange(4)) + assert np.all(matrix == 100.) + + +def test_af2_structure_choice_and_provenance(af2_run): + directory, models = af2_run + model = models[0] + score(directory) + relaxed = directory / f"relaxed_{model}.pdb" + write_structure(relaxed, plddt=40.) + assert Path(BaseParser._guess_struct(directory, model)) == relaxed + row = score(directory)[0] + assert float(row["interface_average_plddt"]) == 40. + assert row["structure_file"] == relaxed.name + + +def test_af2_structure_lookup_does_not_match_model_10(tmp_path): + write_structure(tmp_path / "unrelaxed_model_10.pdb") + with pytest.raises(ValueError, match="model_1"): + BaseParser._guess_struct(tmp_path, "model_1") + + +def test_af2_relaxation_precedence_across_formats_and_rank_fallback(af2_run): + directory, models = af2_run + model = models[0] + unrelaxed = directory / f"unrelaxed_{model}.pdb" + writer = MMCIFIO() + writer.set_structure(BaseParser._load_structure(unrelaxed)) + unrelaxed_cif = unrelaxed.with_suffix(".cif") + writer.save(str(unrelaxed_cif)) + relaxed = directory / f"relaxed_{model}.pdb" + write_structure(relaxed, plddt=40.) + assert Path(BaseParser._guess_struct(directory, model)) == relaxed + unrelaxed.unlink() + unrelaxed_cif.unlink() + relaxed.rename(directory / "ranked_0.pdb") + row = score(directory)[0] + assert row["structure_file"] == "ranked_0.pdb" + assert float(row["interface_average_plddt"]) == 40. + + +def test_legacy_backend_inference_still_works(): + assert meta_score.infer_backend({"model_used": "model_1_multimer_v3_pred_0"}) == "af2" + assert meta_score.infer_backend({"model_used": "seed-1_sample-0"}) == "af3" + assert meta_score.infer_backend({"backend": "af3", "model_used": "model_misleading"}) == "af3" + assert report._detect_backend([{"model_used": "model_1_multimer_v3_pred_0"}]) == "AlphaFold 2" + + +def test_failed_model_is_retried_and_logs_traceback(af2_run, caplog, record_process): + directory, models = af2_run + broken = directory / f"pae_{models[1]}.json" + broken.write_text("{BROKEN") + rows = score(directory, models_to_analyse="all") + assert {r["model_used"] for r in rows} == {models[0]} + assert any(r.exc_info for r in caplog.records) + assert str(broken) in caplog.text + calls = record_process() + score(directory, models_to_analyse="all") + assert calls == [1] + + +def test_failure_after_one_interface_does_not_leave_partial_model_rows(af2_run, monkeypatch): + directory, models = af2_run + write_structure(directory / f"unrelaxed_{models[0]}.pdb", chains="ABC") + (directory / f"pae_{models[0]}.json").write_text(json.dumps([ + {"predicted_aligned_error": np.full((6, 6), 2.).tolist()} + ])) + original = runner._interface_row + calls = [] + def fail_second_interface(*args, **kwargs): + calls.append(1) + if len(calls) == 2: + raise ValueError("synthetic second-interface failure") + return original(*args, **kwargs) + monkeypatch.setattr(runner, "_interface_row", fail_second_interface) + rows = score(directory, models_to_analyse="all") + assert len(calls) == 3 + assert {row["model_used"] for row in rows} == {models[1]}